[gmx-developers] calc_avcsix is N**2!

[Vishal Vaidyanathan]

Hi,
  I noticed that calc_avcsix does an unnecessary double loop over natoms.
This usually results in a small initialization delay before simulation
begins, but for large systems (~300k atoms), this is very significant
(several hours).

  It's relatively trivial to speed it up. The replacement function is
below in case anybody is interested.

Thanks,
Vishal

static double calc_avcsix_fast( FILE *log, real *nbfp, int ntypes,
            int natoms, int type[],  bool bBHAM ) {

    int * typecounts;
    int i, j;
    float c6,c6sum, c6av;

    typecounts = ( int * ) calloc ( ntypes, sizeof( int ) );


    for ( i = 0; i < natoms; i++ ){
        typecounts[ type[i] ]++;
    }

    //Now we have to do a double loop over ntypes
    c6sum = 0;
    for ( i = 0; i < ntypes; i++ ){
        for ( j = 0; j < ntypes; j++ ) {
            if ( bBHAM )
                c6 = BHAMC( nbfp, ntypes, i, j );
            else
                c6 = C6( nbfp, ntypes, i, j );

            c6sum += typecounts[ i ] * typecounts[ j ] * c6;
        }
    }

    c6av = c6sum / ( natoms * natoms );
    free( typecounts );
    return c6av;
}