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
_______________________________________________
gmx-users mailing list [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
