[gmx-developers] question about variables in calc_f_el() in sim_util.c

Bernhard Reuter b.reuter at uni-kassel.de
Wed Feb 19 11:12:01 CET 2014

Dear all,

I have a question concerning the variables passed to and used in 
calc_f_el() in sim_util.c:
There are two arrays containing the spatial - Ex[ ] - and temporal - Et[ 
] - parameters of the applied electric field.
How to define the values of Ex[ ] is described in the users manual - as 
I understand it:
Ex[ ] is a 3x3-array containing for every dimension (x,y,z) three 
entries like described in the manual.
In the source code the first entry of every column (m=0,1,2=x,y,z) would 
be Ex[m].n, the second Ex[m].a[0] and the third is not used.
In the case of Et[ ] its not that easy, its not described in the manual 
and in the source code one finds four (!) entries for every column of 
the array
which would be a 3x4-array.
This would be very confusing since Et[] (as Ex[]) are declared as type 
t_cosines which is defined as

typedef struct {
     int  n;            /* Number of terms                        */
     real *a;           /* Coeffients (V / nm )                   */
     real *phi;         /* Phase angles                                  */
} t_cosines;

in line 47 of inputrec.h.
Deducing from this definition a column of Et[ ] should only contain 
three entries (Et[m].n Et[m].a Et[m].phi) but obviously it has four. 
Worse Et[m].phi is not even used but therefore a now seems to be an 
array with a[0], a[1] and a[2]. How can that be?
I'm totaly confused by this - could anyone please enlighten me!

Thank you very much and best regards,

source code of calc_f_el from line 329 in sim_util.c:

static void calc_f_el(FILE *fp, int  start, int homenr,
                       real charge[], rvec x[], rvec f[],
                       t_cosines Ex[], t_cosines Et[], double t)
     rvec Ext;
     real t0;
     int  i, m;

     for (m = 0; (m < DIM); m++)
         if (Et[m].n > 0)
             if (Et[m].n == 3)
                 t0     = Et[m].a[1];
                 Ext[m] = 
                 Ext[m] = cos(Et[m].a[0]*t);
             Ext[m] = 1.0;
         if (Ex[m].n > 0)
             /* Convert the field strength from V/nm to MD-units */
             Ext[m] *= Ex[m].a[0]*FIELDFAC;
             for (i = start; (i < start+homenr); i++)
                 f[i][m] += charge[i]*Ext[m];
             Ext[m] = 0;
     if (fp != NULL)
         fprintf(fp, "%10g  %10g  %10g  %10g #FIELD\n", t,
                 Ext[XX]/FIELDFAC, Ext[YY]/FIELDFAC, Ext[ZZ]/FIELDFAC);

More information about the gromacs.org_gmx-developers mailing list