jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

avoid division by zero

This commit is contained in:
Axel Kohlmeyer 2020-06-07 15:47:46 -04:00
parent fc1be785fc
commit 4da8ff3d21
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
1 changed files with 2 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/USER-QTB/fix_qbmsst.cpp
View File

@ -520,8 +520,8 @@ void FixQBMSST::initial_integrate(int /*vflag*/)
// decide if the qtb temperature need to be updated or not
if (counter_l == 0) {
t_current -= dtv*fric_coef*eta*beta*(old_eavg-e0)/(3*ntotal*boltz);
if (t_current <= 0.0) break;
old_eavg = 0;//clear old energy average
if (t_current < 0.0) t_current=0;
// load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H
for (int k = 0; k < 2*N_f; k++) {
@ -529,7 +529,7 @@ void FixQBMSST::initial_integrate(int /*vflag*/)
if(k == N_f) {
omega_H[k]=sqrt(force->boltz * t_current);
} else {
double energy_k= force->hplanck * fabs(f_k);
double energy_k= force->hplanck * fabs(f_k);
omega_H[k]=sqrt( energy_k * (0.5+1.0/( exp(energy_k/(force->boltz * t_current)) - 1.0 )) );
omega_H[k]*=alpha*sin((k-N_f)*MY_PI/(2*alpha*N_f))/sin((k-N_f)*MY_PI/(2*N_f));
}