[gmx-developers] Verlet Algorithm

[Elena della Valle]

Hi to all,

I'm Elena della Valle, a Ph.D. student coming from la Sapienza 
University of Rome.

I'm writing because I have some questions about some modifications that 
I have done to the update of velocities and positions in the verlet 
algorithm. My aim is to implement the magnetic field in gromacs by 
introducing the therm of the larmor frequency to the velocities and 
positions. By literature I read that this has been done by the update of 
the verlet velocities and positions with the frequency larmor therm. 
These are  the equations in order to be more clear on what i wanted to do:

I did this by modifying in the update.c file in gromacs the 
update_do_vv_vel and update_do_vv_pos ad follows:

  if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
             {
                 v[n][0] = mv1*(mv1*v[n][0] + 
0.5*(w_dt*mv2*f[n][0]))+0.5*accel[ga][0]*dt + 
w_dt*charge[n]*campoB*v[n][1]*mv1*mv1 + 
0.25*dt*invmass[n]*charge[n]*campoB*((w_dt*mv2*mv1*f[n][1]) + 
accel[ga][1]*dt -2*w_dt*charge[n]*campoB*mv1*mv1*v[n][0]);
                 v[n][1] = mv1*(mv1*v[n][1] + 
0.5*(w_dt*mv2*f[n][1]))+0.5*accel[ga][1]*dt + 
w_dt*charge[n]*campoB*v[n][0]*mv1*mv1 - 
0.25*dt*invmass[n]*charge[n]*campoB*((w_dt*mv2*mv1*f[n][0]) + 
accel[ga][0]*dt +2*w_dt*charge[n]*campoB*mv1*mv1*v[n][1]);
                 v[n][2] = mv1*(mv1*v[n][2] + 
0.5*(w_dt*mv2*f[n][2]))+0.5*accel[ga][2]*dt;



             //printf("frame %lf \n", mv1);
            // printf("frame %d: %lf\t%lf%lf\n", n, v[n][0], v[n][1], 
v[n][2]);
             }

in the do_update_vv_pos:

if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
             {
                 xprime[n][0]   = mr1*(mr1*x[n][0]+mr2*dt*v[n][0]) + 
mr1*0.5*dt*(w_dt*mr2*f[n][0] + w_dt*charge[n]*campoB*v[n][1]*mr1);
                 xprime[n][1]   = mr1*(mr1*x[n][1]+mr2*dt*v[n][1]) + 
mr1*0.5*dt*(w_dt*mr2*f[n][1] - w_dt*charge[n]*campoB*v[n][0]*mr1);
                 xprime[n][2]   = mr1*(mr1*x[n][2]+mr2*dt*v[n][2]) + 
mr1*0.5*dt*(w_dt*mr2*f[n][2]);
             }


After that in the grompp file I modified the integrator in md-vv and I 
run the simulation and it seems to use the position and velocities 
modifications.

I want know if you think that this way is correct, if I implemented in 
the right way

And Also I would like to know if in the grompp file I have to modify 
other parameters (I use a beredensed coupling)

Thanks in advance
Sorry for bothering you
Best Regards
Elena della Valle

-- 
Elena della Valle
Ph.D. Student in Electronic Engineering

Department of Information Engineering, Electronics and Telecommunications
Sapienza, University of Rome
via Eudossiana, 18 00184 Rome

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20161004/78b9d599/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 58009 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20161004/78b9d599/attachment-0001.png>