[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