I submitted this to gmx-users earlier but thought this list will me more appropriate.<br><br>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.<br><br>My algorithm:<br>
<br>1. Divide the MD box in 12 slabs along the x-axis. <br>2. Pick slab 3 as my source and slab 9 as sink.<br>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)<br>
4. Exchange their velocity vectors.<br><br>I achieve this by making following changes to the update.c routine:<br><br>**In static void do_update_md(<br><br>*I comment out the position update line.<br><br> v[n][d] = vb;<br>
/*xprime[n][d] = x[n][d]+vb*dt;*/<br><br>*then I exchange velocities as mentioned above<br><br>*and now I update the positions based on these new velocities:<br><br>xprime[n][d] = x[n][d]+v[n][d]*dt;<br><br>
-----------<br>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.<br>
Note that the average temperature does not increase in the equilibrium NVE with same '.mdp' parameters (no velocity exchange).<br><br>Any thoughts?<br><br>Thanks,<br><font color="#888888">Gaurav</font>