[gmx-users] External Electric field in Gromacs

[Tom]

Dear Gromacs Developers or Users

Althoug the menu of verison 4.5.3 claims that the time-dependent of E
field (E_xt; E_yt; E_zt) has not been implemented yet, as I check the
source code of sim_util.c for the part of calc_f_el, i found the cos
of time function:
----------------------------
     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;
        }

--------------------------
What is meaning of exp part?
How to use the functionality of this subroutine by adding the
parameters in mdp file?

Thanks,

Tom

PS:
The details of the source code was pasted in my previouse email as the
following:


---------- Forwarded message ----------
From: Tom <dnaafm at gmail.com>
Date: Thu, 7 Jul 2011 11:22:15 -0500
Subject: Question about Electrostatic field in Gromacs
To: gmx-users at gromacs.org

 Dear Gromacs developers or users:

I have a question about the implementation of the electrostatic field in
Gromacs (gromacs-4.5.3)

On the menu of this version, it looks that the functionality of
time-dependent has not been implemented yet.
--------
E_x ; E_y ; E_z:
If you want to use an electric field in a direction, enter 3 numbers after
the appropriate E_*,
the first number: the number of cosines, only 1 is implemented (with
frequency 0) so enter
1, the second number: the strength of the electric field in V nm-1, the
third number: the
phase of the cosine, you can enter any number here since a cosine of
frequency zero has no
phase.
E_xt; E_yt; E_zt:
not implemented yet
---------
But in the source code of sim_util.c for the part of calc_f_el :
I notice there is a parameter of Ext[m] as cos function of time together
with exp part. I am so confused about it.
Can anyone help explain this : *Ext[m] = cos(Et[m].a[0]*(t-t0))*
exp(-sqr(t-t0)/(2.0*sqr(Et[m].a[2])))
What is meaning of the Ext[m] and its part of exp part?
I am want to make sure whether gromacs-4.5.3 support the time-dependent E
field calculation?*
And how to control these parameters in mdp file if it supports?
*
*------------------------
        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);
            }
        }

-----------------------------

Thanks,

Tom,

PS: the source code of of sim_util.c for the part of calc_f_el is attached
here:

 * calc_f_el calculates forces due to an electric field.
 *
 * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
 *
 * Et[] contains the parameters for the time dependent
 * part of the field (not yet used).
 * Ex[] contains the parameters for
 * the spatial dependent part of the field. You can have cool periodic
 * fields in principle, but only a constant field is supported
 * now.
 * The function should return the energy due to the electric field
 * (if any) but for now returns 0.
 *
 * WARNING:
 * There can be problems with the virial.
 * Since the field is not self-consistent this is unavoidable.
 * For neutral molecules the virial is correct within this approximation.
 * For neutral systems with many charged molecules the error is small.
 * But for systems with a net charge or a few charged molecules
 * the error can be significant when the field is high.
 * Solution: implement a self-consitent electric field into PME.
 */
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);
    }
}