[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
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);
}
}
More information about the gromacs.org_gmx-developers
mailing list