[gmx-developers] Verlet Algorithm

[David van der Spoel]

On 04/10/16 15:15, Elena della Valle wrote:
>
> 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
>
>
>
Just a quick question: which version did you use? All development should 
be in the master version and this looks like uses something older.



-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se