ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

changes to fixes that use THERMO_ENERGY

Browse Source
This commit is contained in:
Plimpton
2021-01-21 11:32:11 -07:00
parent f54fd8fa72
commit 182eb35f1a
19 changed files with 151 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files

172
src/fix_temp_csvr.cpp
View File

@ -33,89 +33,12 @@
#include "random_mars.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{NOBIAS,BIAS};
enum{CONSTANT,EQUAL};
double FixTempCSVR::gamdev(const int ia)
{
int j;
double am,e,s,v1,v2,x,y;
if (ia < 1) return 0.0;
if (ia < 6) {
x=1.0;
for (j=1; j<=ia; j++)
x *= random->uniform();
// make certain, that -log() doesn't overflow.
if (x < 2.2250759805e-308)
x = 708.4;
else
x = -log(x);
} else {
restart:
do {
do {
do {
v1 = random->uniform();
v2 = 2.0*random->uniform() - 1.0;
} while (v1*v1 + v2*v2 > 1.0);
y=v2/v1;
am=ia-1;
s=sqrt(2.0*am+1.0);
x=s*y+am;
} while (x <= 0.0);
if (am*log(x/am)-s*y < -700 || v1<0.00001) {
goto restart;
}
e=(1.0+y*y)*exp(am*log(x/am)-s*y);
} while (random->uniform() > e);
}
return x;
}
/* -------------------------------------------------------------------
returns the sum of n independent gaussian noises squared
(i.e. equivalent to summing the square of the return values of nn
calls to gasdev)
---------------------------------------------------------------------- */
double FixTempCSVR::sumnoises(int nn) {
if (nn == 0) {
return 0.0;
} else if (nn == 1) {
const double rr = random->gaussian();
return rr*rr;
} else if (nn % 2 == 0) {
return 2.0 * gamdev(nn / 2);
} else {
const double rr = random->gaussian();
return 2.0 * gamdev((nn-1) / 2) + rr*rr;
}
}
/* -------------------------------------------------------------------
returns the scaling factor for velocities to thermalize
the system so it samples the canonical ensemble
---------------------------------------------------------------------- */
double FixTempCSVR::resamplekin(double ekin_old, double ekin_new) {
const double tdof = temperature->dof;
const double c1 = exp(-update->dt/t_period);
const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof;
const double r1 = random->gaussian();
const double r2 = sumnoises(tdof - 1);
const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2);
return sqrt(scale);
}
/* ---------------------------------------------------------------------- */
FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
@ -129,6 +52,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
restart_global = 1;
nevery = 1;
scalar_flag = 1;
ecouple_flag = 1;
global_freq = nevery;
dynamic_group_allow = 1;
extscalar = 1;
@ -168,7 +92,7 @@ FixTempCSVR::FixTempCSVR(LAMMPS *lmp, int narg, char **arg) :
tflag = 1;
nmax = -1;
energy = 0.0;
ecouple = 0.0;
}
/* ---------------------------------------------------------------------- */
@ -192,7 +116,6 @@ int FixTempCSVR::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
mask |= THERMO_ENERGY;
return mask;
}
@ -200,7 +123,6 @@ int FixTempCSVR::setmask()
void FixTempCSVR::init()
{
// check variable
if (tstr) {
@ -224,7 +146,6 @@ void FixTempCSVR::init()
void FixTempCSVR::end_of_step()
{
// set current t_target
// if variable temp, evaluate variable, wrap with clear/add
@ -285,7 +206,7 @@ void FixTempCSVR::end_of_step()
// tally the kinetic energy transferred between heat bath and system
energy += ekin_old * (1.0 - lamda*lamda);
ecouple += ekin_old * (1.0 - lamda*lamda);
}
/* ---------------------------------------------------------------------- */
@ -320,6 +241,85 @@ int FixTempCSVR::modify_param(int narg, char **arg)
/* ---------------------------------------------------------------------- */
double FixTempCSVR::gamdev(const int ia)
{
int j;
double am,e,s,v1,v2,x,y;
if (ia < 1) return 0.0;
if (ia < 6) {
x=1.0;
for (j=1; j<=ia; j++)
x *= random->uniform();
// make certain, that -log() doesn't overflow.
if (x < 2.2250759805e-308)
x = 708.4;
else
x = -log(x);
} else {
restart:
do {
do {
do {
v1 = random->uniform();
v2 = 2.0*random->uniform() - 1.0;
} while (v1*v1 + v2*v2 > 1.0);
y=v2/v1;
am=ia-1;
s=sqrt(2.0*am+1.0);
x=s*y+am;
} while (x <= 0.0);
if (am*log(x/am)-s*y < -700 || v1<0.00001) {
goto restart;
}
e=(1.0+y*y)*exp(am*log(x/am)-s*y);
} while (random->uniform() > e);
}
return x;
}
/* -------------------------------------------------------------------
returns the sum of n independent gaussian noises squared
(i.e. equivalent to summing the square of the return values of nn
calls to gasdev)
---------------------------------------------------------------------- */
double FixTempCSVR::sumnoises(int nn) {
if (nn == 0) {
return 0.0;
} else if (nn == 1) {
const double rr = random->gaussian();
return rr*rr;
} else if (nn % 2 == 0) {
return 2.0 * gamdev(nn / 2);
} else {
const double rr = random->gaussian();
return 2.0 * gamdev((nn-1) / 2) + rr*rr;
}
}
/* -------------------------------------------------------------------
returns the scaling factor for velocities to thermalize
the system so it samples the canonical ensemble
---------------------------------------------------------------------- */
double FixTempCSVR::resamplekin(double ekin_old, double ekin_new) {
const double tdof = temperature->dof;
const double c1 = exp(-update->dt/t_period);
const double c2 = (1.0-c1)*ekin_new/ekin_old/tdof;
const double r1 = random->gaussian();
const double r2 = sumnoises(tdof - 1);
const double scale = c1 + c2*(r1*r1+r2) + 2.0*r1*sqrt(c1*c2);
return sqrt(scale);
}
/* ---------------------------------------------------------------------- */
void FixTempCSVR::reset_target(double t_new)
{
t_target = t_start = t_stop = t_new;
@ -329,7 +329,7 @@ void FixTempCSVR::reset_target(double t_new)
double FixTempCSVR::compute_scalar()
{
return energy;
return ecouple;
}
/* ----------------------------------------------------------------------
@ -339,11 +339,11 @@ double FixTempCSVR::compute_scalar()
void FixTempCSVR::write_restart(FILE *fp)
{
const int PRNGSIZE = 98+2+3;
int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + energy
int nsize = PRNGSIZE*comm->nprocs+2; // pRNG state per proc + nprocs + ecouple
double *list = nullptr;
if (comm->me == 0) {
list = new double[nsize];
list[0] = energy;
list[0] = ecouple;
list[1] = comm->nprocs;
}
double state[PRNGSIZE];
@ -366,7 +366,7 @@ void FixTempCSVR::restart(char *buf)
{
double *list = (double *) buf;
energy = list[0];
ecouple = list[0];
int nprocs = (int) list[1];
if (nprocs != comm->nprocs) {
if (comm->me == 0)