diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 9d88d8f0fb..06a134bc17 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -1580,10 +1580,38 @@ void FixNH::reset_dt() p_freq_max = MAX(p_freq_max,p_freq[4]); p_freq_max = MAX(p_freq_max,p_freq[5]); } + + double kt = boltz * t_target; + double nkt = atom->natoms * kt; + + for (int i = 0; i < 3; i++) + if (p_flag[i]) + omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + } + + // masses and initial forces on barostat thermostat variables + + if (mpchain) { + etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - + boltz * t_target) / etap_mass[ich]; + } + pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); } if (tstat_flag) + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); }