[gmx-users] Re: Calculating electric field vector at snapshots

[Aaron Fafarman]

This is related to a topic from the 15th of July 2008. To recap, I
would like to calculate the electric field vector at a point in space
in a protein at many snapshots during an MD trajectory. I attempted to
do this by constructing a vsite with a very small "test" charge. I got
help adding code to md.c to output the force on the virtual site at
every timestep. This works great for coulombtype = reaction-field:  I
get perfectly coulombic behavior as I move a test charge away from the
vsite. BUT if I employ PME, the results deviate from coulombic
behavior at distances greater than 0.2 nm, quickly decaying to zero
much faster than 1/R^2. (Imagine, if you will,  my dismay at
discovering this after many weeks worth of production runs...).

I am at a loss both as to why this is happening and how to fix it. I
see in md.c that the spread_vsite_f function is called differently in
the case of ewald-based electrostatics, but I'm outputting the force
on the vsite before any force spreading.

I'm ready to give-up on using the vsite facility to calculate forces
unless any saint out there wants to peruse the code below to uncover
the problem.

As an alternative, I'm interested in following up on what David van
der Spoel suggested below, to simply extract the electrostatic forces
on the atoms of interest in every MD step:

>
> The proper way to do this would be to recompute the non-bonded forces
> (short range only) with LJ parameters set to zero. The contribution to
> potential and field from PME can be extracted without too much trouble.
>

My first question is why a "recompute" step would be required?
Shouldn't it be possible to just take the forces from the PME
calculation directly? Isn't the calculation of electrostatic force
done independantly and then added to the LJ and bonded forces? Second,
why short range only? I want the total electrostatic force on the
system including that from the long distance (i.e. reciprocal space)
term.

Thanks in advance for your help.

-Aaron



Excerpt of modified md.c  (see section between comments):

-------------------
    if (bTCR && bFirstStep) {
      tcr=init_coupling(log,nfile,fnm,cr,fr,mdatoms,&(top->idef));
      fprintf(log,"Done init_coupling\n");
      fflush(log);
    }

    /* Dan Ensign 08-25/08 */
    /* the forces are calculated; before spreading the vsite forces to
parent atoms,
       print the force on the atom specified in userint3 (ie, a test charge). */
    rvec_sub( state->x[inputrec->userint1-1],
state->x[inputrec->userint2-1], diff );
    forceproj = iprod( f[inputrec->userint3-1], diff );
    fprintf( log, "force on atom %d %f %f %f\n", inputrec->userint3-1,
f[inputrec->userint3-1][XX], f[inputrec->userint3-1][YY],
f[inputrec->userint3-1][ZZ] );
    fprintf( log, "dot product is %f\n" , forceproj );

    /* Now we have the energies and forces corresponding to the
     * coordinates at time t. We must output all of this before
     * the update.
     * for RerunMD t is read from input trajectory
     */

    if (bVsites)
      spread_vsite_f(log,state->x,f,&mynrnb,&top->idef,
                     fr,graph,state->box,vsitecomm,cr);

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

md.mdp :

---------------------
integrator      = md
nsteps          = 1
dt              = 0.002
nstlist         = 10
nstcomm         = 1
rlist           = 0.9
coulombtype     = pme
fourierspacing  = 0.12
pme_order       = 6
ewald_rtol      = 1e-5
rcoulomb        = 0.9
vdw-type        = cut-off
rvdw            = 0.9
tcoupl          = Nose-Hoover
tc-grps         = protein non-protein
tau-t           = 0.5 0.5
ref-t           = 298 298
nstfout         = 1
nstxout         = 100
nstxtcout       = 100
nstenergy       = 100
userint1        = 1
userint2        = 2
userint3        = 4

_____________