# [gmx-users] Re: SD temperature coupling groups

A.Rzepiela A.Rzepiela at rug.nl
Tue Jun 30 13:29:11 CEST 2009

```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

```