[gmx-users] Re: SD temperature coupling groups

[A.Rzepiela]

Hey!

I think Your problem is from function init_stochd in update.c 
(src/mdlib).
When You set tau_t to 0 You got division by 0 on line 288:

for(n=0; n<ngtc; n++) {
       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);

I think the easiest solution, I checked it for sd1, is to change in 
do_update_sd1 (update.c) following lines

     for(d=0; d<DIM; d++) {
       if((ptype[n] != eptVSite) && (ptype[n] != eptShell) && 
!nFreeze[gf][d]) {
         sd_V = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);

         v[n][d] = v[n][d]*sdc[gt].em
           + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - 
sdc[gt].em)
           + sd_V;

         xprime[n][d] = x[n][d] + v[n][d]*dt;
       } else {
         v[n][d]      = 0.0;
         xprime[n][d] = x[n][d];
       }
     }


to

    for(d=0; d<DIM; d++) {
       if((ptype[n] != eptVSite) && (ptype[n] != eptShell) && 
!nFreeze[gf][d] && (tau_t[gt] != 0.0)) {
         sd_V = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);
         v[n][d] = v[n][d]*sdc[gt].em
           + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - 
sdc[gt].em)
           + sd_V;

         xprime[n][d] = x[n][d] + v[n][d]*dt;
       } else  {
               if((ptype[n] != eptVSite) && (ptype[n] != eptShell) && 
!nFreeze[gf][d]){
                  v[n][d]      = v[n][d] + (invmass[n]*f[n][d] + 
accel[ga][d])*dt;
                  xprime[n][d] = x[n][d] + v[n][d]*dt ;
               }
               else {
                  v[n][d] = 0;
                  xprime[n][d] = x[n][d]; }
       }
     }

You can do similar change  for sd2.
Hope somebody corrects me if I am wrong.
Greetings
Andrzej