[gmx-developers] External forces or field

[S. Joseph]

Hi,

I tried modifying  gromacs-4.03 to add  an external electric field
that varies with time and space. I checked the mailing lists about
adding external forces and there was an email from Gerrit Groenhoff
who mentioned that external forces could be added in force.c.  I
changed the function calc_f_el called by do_force in mdlib/sim_util.c.
Instead of the Ext field used by gromacs, I added the force due to my
external field using positions rvec x and time t. I assumed the
positions were the actual positions (not shifted by PBC).  Now I am
testing the system that has  a single Cl- ion with my external field
acting on it with a cut-off (no PME) for electrostatics. I give it an
initial velocity of 0.5 nm/ps in the conf.gro file used by grompp.

If I have a frequency of oscillation of 100s of GHz of the field and
from the field I gave I expect the ion trajectory to stabilize at a
particular point based on theory. (This simulation was done and
published by another person with NAMD). In my case I see this
stabilization only if I give an extermely large timestep of (0.5-2
ps). When I gave a timestep of 0.002 ps, depending on the magnitude of
the field, the ion either did not move at all or moved in some
trajectory which was radically different from expected (In NAMD after
a time of 1 ns, it converged to a point even at a timestep of 2fs, but
with GROMACS the time of convergence depended on my timestep- for
larger timesteps-2ps, it converged within 1 ns, but smaller the time
step, it took much longer - even microseconds until ultimately it
stopped moving at the order of fs).

Also g_energy does not give any kinetic energy for the ion. The
results changed little whether I started from 0 velocity or 0.5 nm/ps
velocity (That is, different temperatures did not give different
results)

Am I doing something wrong in adding forces or is the non neutrality
of the system causing problems. At the beginning of calc_f_el, it is
mentioned that the system has to be neutral for the correct virial,
but I have an NVT system. I changed from single to double precision
with no difference.

I am planning to add water and other molecules to the system (NVT), so
a timestep of 0.5ps would be unthinkable. Any ideas on what is going
wrong?

Regards,
Joseph