[gmx-developers] The do_update_sd2 function for Langevin dynamics
李华平
lihuaping at tju.edu.cn
Fri Nov 29 10:43:48 CET 2013
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20131129/ae07d6d1/attachment-0001.html>
-------------- next part --------------
Hi GROMACS user group!
I've been having a question when I read the code of do_update_sd2 function in update.c file.
......
kT = BOLTZ*ref_t[n];
/* The mass is encounted for later, since this differs per atom */
sig[n].V = sqrt(kT*(1-sdc[n].em));
sig[n].X = sqrt(kT*sqr(tau_t[n])*sdc[n].c);
sig[n].Yv = sqrt(kT*sdc[n].b/sdc[n].c);
sig[n].Yx = sqrt(kT*sqr(tau_t[n])*sdc[n].b/(1-sdc[n].em)); // formula 1
......
And mostly sdc[n].b is calculated in the init_stochd function through the below code:
......
sdc[n].gdt = ir->delta_t/ir->opts.tau_t[n];
sdc[n].eph = exp(sdc[n].gdt/2);
sdc[n].emh = exp(-sdc[n].gdt/2);
sdc[n].em = exp(-sdc[n].gdt);
......
......
if (sdc[n].gdt >= 0.05)
{
sdc[n].b = sdc[n].gdt*(sdc[n].eph*sdc[n].eph - 1)
- 4*(sdc[n].eph - 1)*(sdc[n].eph - 1); // formula 2
sdc[n].c = sdc[n].gdt - 3 + 4*sdc[n].emh - sdc[n].em;
sdc[n].d = 2 - sdc[n].eph - sdc[n].emh;
}
......
I found the code is not consistent with the formula (3.32) in the original paper,
W. F. Van Gunsteren & H. J. C. Berendsen (1988): A Leap-frog Algorithm for Stochastic Dynamics, Molecular Simulation, 1:3, 173-185,
where the B function has a negative parameter. But the above formula 1 (in bold) has a positive parameter.
This contradiction is confusimg me. Is it a bug or just have I made a mistake?
I wonder whether this negligence might influence our simulation results.
Hope somebody corrects me if I am wrong.
Best wishes!
yours sincerely,
Huaping Li
More information about the gromacs.org_gmx-developers
mailing list