Merge pull request #1707 from charlessievers/lammps_gjf

Lammps gjf small updates
This commit is contained in:
Axel Kohlmeyer
2019-10-08 23:06:36 +02:00
committed by GitHub
2 changed files with 20 additions and 11 deletions

View File

@ -246,11 +246,11 @@ velocity given in "Gronbech-Jensen/Farago"_#Gronbech-Jensen; this velocity
is shown to be systematically lower than the target temperature by a small is shown to be systematically lower than the target temperature by a small
amount, which grows quadratically with the timestep. amount, which grows quadratically with the timestep.
The {gjf} option {vhalf} outputs the 2GJ half-step velocity given in The {gjf} option {vhalf} outputs the 2GJ half-step velocity given in
"Gronbech Jensen/Gronbech-Jensen"_#2Gronbech-Jensen; this velocity is shown "Gronbech Jensen/Gronbech-Jensen"_#2Gronbech-Jensen; for linear systems,
to not have any linear statistical errors for any stable time step. this velocity is shown to not have any statistical errors for any stable time step.
An overview of statistically correct Boltzmann and Maxwell-Boltzmann An overview of statistically correct Boltzmann and Maxwell-Boltzmann
sampling of true on-site and true half-step velocities is given in sampling of true on-site and true half-step velocities is given in
"Gronbech-Jensen_#1Gronbech-Jensen. "Gronbech-Jensen"_#1Gronbech-Jensen.
Regardless of the choice of output velocity, the sampling of the configurational Regardless of the choice of output velocity, the sampling of the configurational
distribution of atom positions is the same, and linearly consistent with the distribution of atom positions is the same, and linearly consistent with the
target temperature. target temperature.

View File

@ -240,7 +240,7 @@ void FixLangevin::init()
if (strcmp(id,modify->fix[i]->id) == 0) before = 0; if (strcmp(id,modify->fix[i]->id) == 0) before = 0;
else if ((modify->fmask[i] && utils::strmatch(modify->fix[i]->style,"^nve")) && before) flag = 1; else if ((modify->fmask[i] && utils::strmatch(modify->fix[i]->style,"^nve")) && before) flag = 1;
} }
if (flag && comm->me == 0) if (flag)
error->all(FLERR,"Fix langevin gjf should come before fix nve"); error->all(FLERR,"Fix langevin gjf should come before fix nve");
} }
@ -295,12 +295,12 @@ void FixLangevin::init()
gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v; gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
if (gjfflag) if (gjfflag)
gfactor2[i] = sqrt(atom->mass[i]) * gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) / sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v; force->ftm2v;
else else
gfactor2[i] = sqrt(atom->mass[i]) * gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) / sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v; force->ftm2v;
gfactor1[i] *= 1.0/ratio[i]; gfactor1[i] *= 1.0/ratio[i];
gfactor2[i] *= 1.0/sqrt(ratio[i]); gfactor2[i] *= 1.0/sqrt(ratio[i]);
} }
@ -1007,12 +1007,21 @@ void FixLangevin::reset_dt()
{ {
if (atom->mass) { if (atom->mass) {
for (int i = 1; i <= atom->ntypes; i++) { for (int i = 1; i <= atom->ntypes; i++) {
gfactor2[i] = sqrt(atom->mass[i]) * if (gjfflag)
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) / gfactor2[i] = sqrt(atom->mass[i]) *
force->ftm2v; sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
else
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
gfactor2[i] *= 1.0/sqrt(ratio[i]); gfactor2[i] *= 1.0/sqrt(ratio[i]);
} }
} }
if (gjfflag) {
gjfa = (1.0-update->dt/2.0/t_period)/(1.0+update->dt/2.0/t_period);
gjfsib = sqrt(1.0+update->dt/2.0/t_period);
}
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */