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

[Bernhard Reuter]

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
Et[m].n
Et[m].a[0]
Et[m].a[1]
Et[m].a[2]
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,
Bernhard

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] = 
cos(Et[m].a[0]*(t-t0))*exp(-sqr(t-t0)/(2.0*sqr(Et[m].a[2])));
             }
             else
             {
                 Ext[m] = cos(Et[m].a[0]*t);
             }
         }
         else
         {
             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];
             }
         }
         else
         {
             Ext[m] = 0;
         }
     }
     if (fp != NULL)
     {
         fprintf(fp, "%10g  %10g  %10g  %10g #FIELD\n", t,
                 Ext[XX]/FIELDFAC, Ext[YY]/FIELDFAC, Ext[ZZ]/FIELDFAC);
     }
}