[gmx-developers] Nonequilibrium MD- implementing a temperature gradient simulation

[Gaurav Goel]

I submitted this to gmx-users earlier but thought this list will me more
appropriate.

My aim is to establish a temperature gradient in the MD simulation box. I'm
simulating a magnesium silicate system with interatomic potentials (Mg-Mg,
Mg-Si, etc.) and no bond or angle potentials.

My algorithm:

1. Divide the MD box in 12 slabs along the x-axis.
2. Pick slab 3 as my source and slab 9 as sink.
3. Find the atom with highest kinetic energy in slab 3 and the atom with
lowest KE in slab 9 (I amke sure both atoms have same type and hence same
mass)
4. Exchange their velocity vectors.

I achieve this by making following changes to the update.c routine:

**In static void do_update_md(

*I comment out the position update line.

 v[n][d]        = vb;
          /*xprime[n][d]   = x[n][d]+vb*dt;*/

*then I exchange velocities as mentioned above

*and now I update the positions based on these new velocities:

xprime[n][d]   = x[n][d]+v[n][d]*dt;

-----------
I print out the kinetic energy of atoms before and after exchanging and see
no problem there. Total KE of these two atoms is the same before and after.
But, over the length of my simulation I see average temperature of the
system keeps on increasing in the NVE ensemble-- in 1ns temperature of the
system increased from 3000K -> 3500K.
Note that the average temperature does not increase in the equilibrium NVE
with same '.mdp' parameters (no velocity exchange).

Any thoughts?

Thanks,
Gaurav
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20110203/5e427c88/attachment.html>