diff --git a/src/EXTRA-FIX/fix_langevin_gjf.cpp b/src/EXTRA-FIX/fix_langevin_gjf.cpp index 7f38a4d201..59d5cff1c5 100644 --- a/src/EXTRA-FIX/fix_langevin_gjf.cpp +++ b/src/EXTRA-FIX/fix_langevin_gjf.cpp @@ -44,8 +44,8 @@ enum { CONSTANT, EQUAL, ATOM }; /* ---------------------------------------------------------------------- */ -FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), gjfflag(0), gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), +FixLangevinGJF::FixLangevinGJF(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), gjfflag(0), tstr(nullptr), flangevin(nullptr), tforce(nullptr), lv(nullptr), id_temp(nullptr), random(nullptr) { if (narg < 8) error->all(FLERR, "Illegal fix langevin/gjf command"); @@ -78,16 +78,10 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : random = new RanMars(lmp, seed + comm->me); - // allocate per-type arrays for force prefactors - - // gfactor1 = new double[atom->ntypes + 1]; - // gfactor2 = new double[atom->ntypes + 1]; - // ratio = new double[atom->ntypes + 1]; int GJmethods = 8 // number of currently implemented GJ methods // optional args - for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0; osflag = 0; GJmethod = 0; @@ -105,39 +99,28 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg], "method") == 0) { if (iarg + 2 > narg) error->all(FLERR, "Illegal fix langevin/gjf command"); GJmethod = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); - if (GJmethod <= 0 || GJmethod > GJmethods) error->all(FLERR, "Invalid GJ method choice in langevin/gjf command"); + if (GJmethod < 0 || GJmethod > GJmethods) error->all(FLERR, "Invalid GJ method choice in langevin/gjf command"); iarg += 2; } else error->all(FLERR, "Illegal fix langevin/gjf command"); } // set temperature = nullptr, user can override via fix_modify if wants bias - id_temp = nullptr; temperature = nullptr; - energy = 0.0; - - // flangevin is unallocated until first call to setup() - // compute_scalar checks for this and returns 0.0 - // if flangevin_allocated is not set - - flangevin = nullptr; - flangevin_allocated = 0; lv = nullptr; tforce = nullptr; - maxatom1 = maxatom2 = 0; // setup atom-based array for lv // register with Atom class // no need to set peratom_flag, b/c data is for internal use only - FixLangevin::grow_arrays(atom->nmax); + FixLangevinGJF::grow_arrays(atom->nmax); atom->add_callback(Atom::GROW); // initialize lv to zero - int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { lv[i][0] = 0.0; @@ -148,17 +131,13 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixLangevin::~FixLangevin() +FixLangevinGJF::~FixLangevinGJF() { if (copymode) return; delete random; delete[] tstr; - delete[] gfactor1; - delete[] gfactor2; - delete[] ratio; delete[] id_temp; - memory->destroy(flangevin); memory->destroy(tforce); memory->destroy(lv); @@ -167,18 +146,18 @@ FixLangevin::~FixLangevin() /* ---------------------------------------------------------------------- */ -int FixLangevin::setmask() +int FixLangevinGJF::setmask() { int mask = 0; mask |= INITIAL_INTEGRATE; mask |= FINAL_INTEGRATE; - mask |= END_OF_STEP; + if (!osflag) mask |= END_OF_STEP; return mask; } /* ---------------------------------------------------------------------- */ -void FixLangevin::init() +void FixLangevinGJF::init() { if (id_temp) { temperature = modify->get_compute_by_id(id_temp); @@ -193,23 +172,13 @@ void FixLangevin::init() if (tstr) { tvar = input->variable->find(tstr); - if (tvar < 0) error->all(FLERR, "Variable name {} for fix langevin does not exist", tstr); + if (tvar < 0) error->all(FLERR, "Variable name {} for fix langevin/gjf does not exist", tstr); if (input->variable->equalstyle(tvar)) tstyle = EQUAL; else if (input->variable->atomstyle(tvar)) tstyle = ATOM; else - error->all(FLERR, "Variable {} for fix langevin is invalid style", tstr); - } - - // set force prefactors - - if (!atom->rmass) { - for (int i = 1; i <= atom->ntypes; i++) { - gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v; - gfactor2[i] = sqrt(atom->mass[i]) / force->ftm2v; - gfactor2[i] *= sqrt(2.0 * update->dt * force->boltz / t_period / force->mvv2e); // gjfflag - } + error->all(FLERR, "Variable {} for fix langevin/gjf is invalid style", tstr); } if (temperature && temperature->tempbias) @@ -219,44 +188,44 @@ void FixLangevin::init() if (utils::strmatch(update->integrate_style, "^respa")) { nlevels_respa = (static_cast(update->integrate))->nlevels; - if (gjfflag) error->all(FLERR, "Fix langevin gjf and run style respa are not compatible"); - } - - if (gjfflag) { - gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period); - gjfc1 = 1.0 / (1.0 + update->dt / 2.0 / t_period); + error->all(FLERR, "Fix langevin gjf and run style respa are not compatible"); } + // Complete set of thermostats is given in Gronbech-Jensen, Molecular Physics, 118 (2020) switch (GJmethod) { case 1: gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period); - gjfc1 = 1.0 / (1.0 + update->dt / 2.0 / t_period); break; case 2: - // Insert logic for method 2 + gjfc2 = exp(-update->dt / t_period); break; case 3: - // Insert logic for method 3 + gjfc2 = 1.0 - update->dt / t_period; break; case 4: - // Insert logic for method 4 + gjfc2 = ( sqrt(1.0 + 4.0 * (update->dt / t_period) ) - 1.0 ) / ( 2.0 * update->dt / t_period ); break; case 5: - // Insert logic for method 5 + gjfc2 = 1.0 / (1.0 + update->dt / t_period); break; case 6: - // Insert logic for method 6 + gjfc2 = (1.0 / (1.0 + update->dt / 2.0 / t_period)) * (1.0 / (1.0 + update->dt / 2.0 / t_period)); break; - case 7: - // Insert logic for method 7 + case 7: // provided in Finkelstein (2021) + gjfc2 = 1; // TODO: correct this break; - case 8: - // Insert logic for method 8 + case 8: // provided in Gronbech-Jensen (2024) + gjfc2 = sqrt( (update->dt / t_period)*(update->dt / t_period) + 1.0 ) - update->dt / t_period; + break; + case 0: + gjfc2 = 0.0; break; default: error->all(FLERR, "Fix langevin/gjf method not found"); break; -} + } + gjfc1 = (1.0 + gjfc2) / 2.0; + gjfc3 = (1.0 - gjfc2) * t_period / update->dt; } /* ---------------------------------------------------------------------- @@ -268,10 +237,9 @@ void FixLangevin::init() 4. Velocity Choice in end_of_step() ------------------------------------------------------------------------- */ -void FixLangevin::initial_integrate(int /* vflag */) +void FixLangevinGJF::initial_integrate(int /* vflag */) { - double gamma1,gamma2; - + // This function provides the integration of the GJ formulation 24a-e double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -290,79 +258,21 @@ void FixLangevin::initial_integrate(int /* vflag */) double dtf = 0.5 * dt * ftm2v; double dtfm; - double c1sqrt = sqrt(gjfc1); - - // NVE integrates position and velocity according to Eq. 8a, 8b - // This function embeds the GJF formulation into the NVE framework, which corresponds to the GJF case c1=c3. + double c1sq = sqrt(gjfc1); + double c3sq = sqrt(gjfc3); + double csq = sqrt(gjfc3 / gjfc1); + double m, beta; - //NVE - if (rmass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dt * v[i][0]; - x[i][1] += dt * v[i][1]; - x[i][2] += dt * v[i][2]; - } - - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; - x[i][0] += dt * v[i][0]; - x[i][1] += dt * v[i][1]; - x[i][2] += dt * v[i][2]; - } - } - - // The initial NVE integration should always use the on-site velocity. Therefore, a velocity correction - // must be done when using the half-step option. - //---------- + // If user elected vhalf, v needs to be reassigned to onsite velocity for integration if (!osflag) { - if (rmass) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - // Undo NVE integration - x[i][0] -= dt * v[i][0]; - x[i][1] -= dt * v[i][1]; - x[i][2] -= dt * v[i][2]; - // Obtain Eq. 24a. lv[][] stores on-site velocity from previous timestep - v[i][0] = lv[i][0] + dtfm * f[i][0]; - v[i][1] = lv[i][1] + dtfm * f[i][1]; - v[i][2] = lv[i][2] + dtfm * f[i][2]; - // Redo NVE integration with correct velocity - x[i][0] += dt * v[i][0]; - x[i][1] += dt * v[i][1]; - x[i][2] += dt * v[i][2]; - } - - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - // Undo NVE integration - x[i][0] -= dt * v[i][0]; - x[i][1] -= dt * v[i][1]; - x[i][2] -= dt * v[i][2]; - // Obtain Eq. 24a. lv[][] stores on-site velocity from previous timestep - v[i][0] = lv[i][0] + dtfm * f[i][0]; - v[i][1] = lv[i][1] + dtfm * f[i][1]; - v[i][2] = lv[i][2] + dtfm * f[i][2]; - // Redo NVE integration with correct velocity - x[i][0] += dt * v[i][0]; - x[i][1] += dt * v[i][1]; - x[i][2] += dt * v[i][2]; - } - } + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + // lv is Eq. 24f from previous time step + v[i][0] = lv[i][0]; + v[i][1] = lv[i][1]; + v[i][2] = lv[i][2]; + } } - //---------- compute_target(); @@ -370,59 +280,58 @@ void FixLangevin::initial_integrate(int /* vflag */) for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { + if (tstyle == ATOM) tsqrt = sqrt(tforce[i]); - if (rmass) { - gamma2 = sqrt(rmass[i]) * sqrt(2.0*dt*boltz/t_period/mvv2e) / ftm2v; - gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt; + if (rmass) { + m = rmass[i]; + beta = tsqrt * sqrt(2.0*dt*rmass[i]*boltz/t_period/mvv2e) / ftm2v; } else { - gamma2 = gfactor2[type[i]] * tsqrt; + m = mass[type[i]]; + beta = tsqrt * sqrt(2.0*dt*atom->mass[i]*boltz/t_period/mvv2e) / ftm2v; } - fran[0] = gamma2*random->gaussian(); - fran[1] = gamma2*random->gaussian(); - fran[2] = gamma2*random->gaussian(); + + fran[0] = beta*random->gaussian(); + fran[1] = beta*random->gaussian(); + fran[2] = beta*random->gaussian(); + + // First integration delivers Eq. 24a and 24b: + dtfm = dtf / m; + v[i][0] += csq * dtfm * f[i][0]; + v[i][1] += csq * dtfm * f[i][1]; + v[i][2] += csq * dtfm * f[i][2]; + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; - // NVE integrator delivers Eq. 24a, but also overshoots position integration. Calculate Eq. 24b: - x[i][0] -= 0.5 * dt * v[i][0]; - x[i][1] -= 0.5 * dt * v[i][1]; - x[i][2] -= 0.5 * dt * v[i][2]; // Calculate Eq. 24c: if (tbiasflag == BIAS) temperature->remove_bias(i,v[i]); - if (rmass) { - lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c1sqrt / (2.0 * rmass[i])) * fran[0]; - lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c1sqrt / (2.0 * rmass[i])) * fran[1]; - lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c1sqrt / (2.0 * rmass[i])) * fran[2]; - } else { - lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c1sqrt / (2.0 * mass[type[i]])) * fran[0]; - lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c1sqrt / (2.0 * mass[type[i]])) * fran[1]; - lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c1sqrt / (2.0 * mass[type[i]])) * fran[2]; - } + lv[i][0] = c1sqrt*v[i][0] + ftm2v * (c3sqrt / (2.0 * m)) * fran[0]; + lv[i][1] = c1sqrt*v[i][1] + ftm2v * (c3sqrt / (2.0 * m)) * fran[1]; + lv[i][2] = c1sqrt*v[i][2] + ftm2v * (c3sqrt / (2.0 * m)) * fran[2]; + if (tbiasflag == BIAS) temperature->restore_bias(i,v[i]); if (tbiasflag == BIAS) temperature->restore_bias(i,lv[i]); // Calculate Eq. 24d - if (tbiasflag == BIAS) temperature->remove_bias(i, lv[i]); - if (atom->rmass) { - v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * (0.5 / rmass[i]) * fran[0]; - v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * (0.5 / rmass[i]) * fran[1]; - v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * (0.5 / rmass[i]) * fran[2]; - } else { - v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * (0.5 / mass[type[i]]) * fran[0]; - v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * (0.5 / mass[type[i]]) * fran[1]; - v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * (0.5 / mass[type[i]]) * fran[2]; - } - if (tbiasflag == BIAS) temperature->restore_bias(i, lv[i]); - // Calculate Eq. 24e. NVE integrator then calculates Eq. 24f. - x[i][0] += 0.5 * dt * v[i][0]; - x[i][1] += 0.5 * dt * v[i][1]; - x[i][2] += 0.5 * dt * v[i][2]; + if (tbiasflag == BIAS) temperature->remove_bias(i, v[i]); + + v[i][0] = (gjfc2 / c1sqrt) * lv[i][0] + ftm2v * csq * (0.5 / m) * fran[0]; + v[i][1] = (gjfc2 / c1sqrt) * lv[i][1] + ftm2v * csq * (0.5 / m) * fran[1]; + v[i][2] = (gjfc2 / c1sqrt) * lv[i][2] + ftm2v * csq * (0.5 / m) * fran[2]; + if (tbiasflag == BIAS) temperature->restore_bias(i, v[i]); + + // Calculate Eq. 24e. Final integrator then calculates Eq. 24f after force update. + x[i][0] += 0.5 * csq * dt * v[i][0]; + x[i][1] += 0.5 * csq * dt * v[i][1]; + x[i][2] += 0.5 * csq * dt * v[i][2]; } } } -void FixLangevin::final_integrate() +void FixLangevinGJF::final_integrate() { double dtfm; double dt = update->dt; @@ -444,18 +353,18 @@ void FixLangevin::final_integrate() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dtfm = dtf / rmass[i]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; + v[i][0] += csq * dtfm * f[i][0]; + v[i][1] += csq * dtfm * f[i][1]; + v[i][2] += csq * dtfm * f[i][2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm * f[i][0]; - v[i][1] += dtfm * f[i][1]; - v[i][2] += dtfm * f[i][2]; + v[i][0] += csq * dtfm * f[i][0]; + v[i][1] += csq * dtfm * f[i][1]; + v[i][2] += csq * dtfm * f[i][2]; } } } @@ -464,7 +373,7 @@ void FixLangevin::final_integrate() set current t_target and t_sqrt ------------------------------------------------------------------------- */ -void FixLangevin::compute_target() +void FixLangevinGJF::compute_target() { int *mask = atom->mask; int nlocal = atom->nlocal; @@ -483,7 +392,7 @@ void FixLangevin::compute_target() if (tstyle == EQUAL) { t_target = input->variable->compute_equal(tvar); if (t_target < 0.0) - error->one(FLERR, "Fix langevin variable returned negative temperature"); + error->one(FLERR, "Fix langevin/gjf variable returned negative temperature"); tsqrt = sqrt(t_target); } else { if (atom->nmax > maxatom2) { @@ -495,86 +404,94 @@ void FixLangevin::compute_target() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (tforce[i] < 0.0) - error->one(FLERR, "Fix langevin variable returned negative temperature"); + error->one(FLERR, "Fix langevin/gjf variable returned negative temperature"); } modify->addstep_compute(update->ntimestep + 1); } } /* ---------------------------------------------------------------------- - tally energy transfer to thermal reservoir, select velocity for GJF + select velocity for GJF ------------------------------------------------------------------------- */ -void FixLangevin::end_of_step() +void FixLangevinGJF::end_of_step() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; - energy_onestep = 0.0; - - if (tallyflag) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + - flangevin[i][2]*v[i][2]; - } - - energy += energy_onestep*update->dt; - - // After the NVE integrator delivers 24f, either the on-site or half-step + // After the final integrator delivers 24f, either the on-site or half-step // velocity is used in remaining simulation tasks, depending on user input - if (gjfflag && !osflag) { - double tmp[3]; - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - // v is Eq. 24f - tmp[0] = v[i][0]; - tmp[1] = v[i][1]; - tmp[2] = v[i][2]; - // Move on with half-step velocity - v[i][0] = lv[i][0]; - v[i][1] = lv[i][1]; - v[i][2] = lv[i][2]; - // store Eq. 24f in lv for next timestep - lv[i][0] = tmp[0]; - lv[i][1] = tmp[1]; - lv[i][2] = tmp[2]; - } - } + double tmp[3]; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + // v is Eq. 24f + tmp[0] = v[i][0]; + tmp[1] = v[i][1]; + tmp[2] = v[i][2]; + // Move on with half-step velocity + v[i][0] = lv[i][0]; + v[i][1] = lv[i][1]; + v[i][2] = lv[i][2]; + // store Eq. 24f in lv for next timestep + lv[i][0] = tmp[0]; + lv[i][1] = tmp[1]; + lv[i][2] = tmp[2]; + } } // clang-format on /* ---------------------------------------------------------------------- */ -void FixLangevin::reset_target(double t_new) +void FixLangevinGJF::reset_target(double t_new) { t_target = t_start = t_stop = t_new; } /* ---------------------------------------------------------------------- */ -void FixLangevin::reset_dt() +void FixLangevinGJF::reset_dt() { - if (atom->mass) { - for (int i = 1; i <= atom->ntypes; i++) { - gfactor2[i] = sqrt(atom->mass[i]) / force->ftm2v; - if (gjfflag) - gfactor2[i] *= sqrt(2.0 * update->dt * force->boltz / t_period / force->mvv2e); // sqrt(2*alpha*kT*dt) - else - gfactor2[i] *= sqrt(24.0 * force->boltz / t_period / update->dt / force->mvv2e); - gfactor2[i] *= 1.0 / sqrt(ratio[i]); - } - } - if (gjfflag) { - gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period); - gjfc1 = 1.0 / (1.0 + update->dt / 2.0 / t_period); + // Complete set of thermostats is given in Gronbech-Jensen, Molecular Physics, 118 (2020) + switch (GJmethod) { + case 1: + gjfc2 = (1.0 - update->dt / 2.0 / t_period) / (1.0 + update->dt / 2.0 / t_period); + break; + case 2: + gjfc2 = exp(-update->dt / t_period); + break; + case 3: + gjfc2 = 1.0 - update->dt / t_period; + break; + case 4: + gjfc2 = ( sqrt(1.0 + 4.0 * (update->dt / t_period) ) - 1.0 ) / ( 2.0 * update->dt / t_period ); + break; + case 5: + gjfc2 = 1.0 / (1.0 + update->dt / t_period); + break; + case 6: + gjfc2 = (1.0 / (1.0 + update->dt / 2.0 / t_period)) * (1.0 / (1.0 + update->dt / 2.0 / t_period)); + break; + case 7: // provided in Finkelstein (2021) + gjfc2 = 1; // TODO: correct this + break; + case 8: // provided in Gronbech-Jensen (2024) + gjfc2 = sqrt( (update->dt / t_period)*(update->dt / t_period) + 1.0 ) - update->dt / t_period; + break; + case 0: + gjfc2 = 0.0; + break; + default: + error->all(FLERR, "Fix langevin/gjf method not found"); + break; } + gjfc1 = (1.0 + gjfc2) / 2.0; + gjfc3 = (1.0 - gjfc2) * t_period / update->dt; } /* ---------------------------------------------------------------------- */ -int FixLangevin::modify_param(int narg, char **arg) +int FixLangevinGJF::modify_param(int narg, char **arg) { if (strcmp(arg[0], "temp") == 0) { if (narg < 2) utils::missing_cmd_args(FLERR, "fix_modify", error); @@ -594,52 +511,11 @@ int FixLangevin::modify_param(int narg, char **arg) return 0; } -/* ---------------------------------------------------------------------- */ - -double FixLangevin::compute_scalar() -{ - if (!tallyflag || !flangevin_allocated) return 0.0; - - // capture the very first energy transfer to thermal reservoir - - double **v = atom->v; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (update->ntimestep == update->beginstep) { - energy_onestep = 0.0; - if (!gjfflag) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - energy_onestep += - flangevin[i][0] * v[i][0] + flangevin[i][1] * v[i][1] + flangevin[i][2] * v[i][2]; - energy = 0.5 * energy_onestep * update->dt; - } else { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (tbiasflag) temperature->remove_bias(i, lv[i]); - energy_onestep += - flangevin[i][0] * lv[i][0] + flangevin[i][1] * lv[i][1] + flangevin[i][2] * lv[i][2]; - if (tbiasflag) temperature->restore_bias(i, lv[i]); - } - energy = -0.5 * energy_onestep * update->dt; - } - } - - // convert midstep energy back to previous fullstep energy - - double energy_me = energy - 0.5 * energy_onestep * update->dt; - - double energy_all; - MPI_Allreduce(&energy_me, &energy_all, 1, MPI_DOUBLE, MPI_SUM, world); - return -energy_all; -} - /* ---------------------------------------------------------------------- extract thermostat properties ------------------------------------------------------------------------- */ -void *FixLangevin::extract(const char *str, int &dim) +void *FixLangevinGJF::extract(const char *str, int &dim) { dim = 0; if (strcmp(str, "t_target") == 0) { return &t_target; } @@ -650,11 +526,10 @@ void *FixLangevin::extract(const char *str, int &dim) memory usage of tally array ------------------------------------------------------------------------- */ -double FixLangevin::memory_usage() +double FixLangevinGJF::memory_usage() { double bytes = 0.0; - if (gjfflag) bytes += (double) atom->nmax * 3 * sizeof(double); - if (tallyflag || osflag) bytes += (double) atom->nmax * 3 * sizeof(double); + bytes += (double) atom->nmax * 3 * sizeof(double); if (tforce) bytes += (double) atom->nmax * sizeof(double); return bytes; } @@ -663,16 +538,16 @@ double FixLangevin::memory_usage() allocate atom-based array for lv ------------------------------------------------------------------------- */ -void FixLangevin::grow_arrays(int nmax) +void FixLangevinGJF::grow_arrays(int nmax) { - memory->grow(lv, nmax, 3, "fix_langevin:lv"); + memory->grow(lv, nmax, 3, "fix_langevin_gjf:lv"); } /* ---------------------------------------------------------------------- copy values within local atom-based array ------------------------------------------------------------------------- */ -void FixLangevin::copy_arrays(int i, int j, int /*delflag*/) +void FixLangevinGJF::copy_arrays(int i, int j, int /*delflag*/) { lv[j][0] = lv[i][0]; lv[j][1] = lv[i][1]; @@ -683,7 +558,7 @@ void FixLangevin::copy_arrays(int i, int j, int /*delflag*/) pack values in local atom-based array for exchange with another proc ------------------------------------------------------------------------- */ -int FixLangevin::pack_exchange(int i, double *buf) +int FixLangevinGJF::pack_exchange(int i, double *buf) { int n = 0; buf[n++] = lv[i][0]; @@ -696,7 +571,7 @@ int FixLangevin::pack_exchange(int i, double *buf) unpack values in local atom-based array from exchange with another proc ------------------------------------------------------------------------- */ -int FixLangevin::unpack_exchange(int nlocal, double *buf) +int FixLangevinGJF::unpack_exchange(int nlocal, double *buf) { int n = 0; lv[nlocal][0] = buf[n++]; diff --git a/src/EXTRA-FIX/fix_langevin_gjf.h b/src/EXTRA-FIX/fix_langevin_gjf.h index c54772c457..df6b80e1b7 100644 --- a/src/EXTRA-FIX/fix_langevin_gjf.h +++ b/src/EXTRA-FIX/fix_langevin_gjf.h @@ -24,21 +24,19 @@ FixStyle(langevin,FixLangevin); namespace LAMMPS_NS { -class FixLangevin : public Fix { +class FixLangevinGJF : public Fix { public: - FixLangevin(class LAMMPS *, int, char **); - ~FixLangevin() override; + FixLangevinGJF(class LAMMPS *, int, char **); + ~FixLangevinGJF() override; int setmask() override; void init() override; void setup(int) override; void initial_integrate(int) override; - void post_force(int) override; - void post_force_respa(int, int, int) override; + void final_integrate() override; void end_of_step() override; void reset_target(double) override; void reset_dt() override; int modify_param(int, char **) override; - double compute_scalar() override; double memory_usage() override; void *extract(const char *, int &) override; void grow_arrays(int) override; @@ -48,19 +46,13 @@ class FixLangevin : public Fix { protected: int osflag, GJmethod; - int flangevin_allocated; double t_start, t_stop, t_period, t_target; - double *gfactor1, *gfactor2, *ratio; - double energy, energy_onestep; + double *gfactor2; double tsqrt; - double gjfc1, gjfc2; + double gjfc1, gjfc2, gjfc3; int tstyle, tvar; char *tstr; - class AtomVecEllipsoid *avec; - - int maxatom1, maxatom2; - double **flangevin; double *tforce; double **lv; //half step velocity @@ -71,11 +63,6 @@ class FixLangevin : public Fix { class RanMars *random; int seed; - template - void post_force_templated(); - - void omega_thermostat(); - void angmom_thermostat(); void compute_target(); };