lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Added GJF Langevin method 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10045 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps 2013-06-06 18:10:38 +00:00
parent 8f2bc5c721
commit fe71627776
2 changed files with 363 additions and 257 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

598
src/fix_langevin.cpp
View File

@ -12,7 +12,8 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Carolyn Phillips (U Mich), reservoir energy tally
Contributing authors: Carolyn Phillips (U Mich), reservoir energy tally
Aidan Thompson (SNL) GJF formulation
------------------------------------------------------------------------- */
#include "mpi.h"
@ -90,9 +91,10 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
// optional args
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
oflag = 0;
ascale = 0.0;
tally = 0;
gjfflag = 0;
oflag = 0;
tallyflag = 0;
zeroflag = 0;
int iarg = 7;
@ -102,6 +104,12 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg+1],"no") == 0) ascale = 0.0;
else ascale = force->numeric(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"gjf") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) gjfflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) gjfflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"omega") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
@ -118,8 +126,8 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
iarg += 3;
} else if (strcmp(arg[iarg],"tally") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) tally = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) tally = 1;
if (strcmp(arg[iarg+1],"no") == 0) tallyflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) tallyflag = 1;
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"zero") == 0) {
@ -141,8 +149,30 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
energy = 0.0;
flangevin = NULL;
franprev = NULL;
tforce = NULL;
maxatom1 = maxatom2 = 0;
// Setup atom-based array for franprev
// register with Atom class
// No need to set peratom_flag
// as this data is for internal use only
if (gjfflag) {
nvalues = 3;
grow_arrays(atom->nlocal);
atom->add_callback(0);
// initialize franprev to zero
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
}
}
}
/* ---------------------------------------------------------------------- */
@ -157,6 +187,11 @@ FixLangevin::~FixLangevin()
delete [] id_temp;
memory->destroy(flangevin);
memory->destroy(tforce);
if (gjfflag) {
memory->destroy(franprev);
atom->delete_callback(id,0);
}
}
/* ---------------------------------------------------------------------- */
@ -166,6 +201,7 @@ int FixLangevin::setmask()
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= POST_INTEGRATE;
mask |= END_OF_STEP;
mask |= THERMO_ENERGY;
return mask;
@ -232,11 +268,14 @@ void FixLangevin::init()
}
}
if (temperature && temperature->tempbias) which = BIAS;
else which = NOBIAS;
if (temperature && temperature->tempbias) tbiasflag = BIAS;
else tbiasflag = NOBIAS;
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (gjfflag) gjffac = 1.0/(1.0+update->dt/2.0/t_period);
}
/* ---------------------------------------------------------------------- */
@ -256,8 +295,143 @@ void FixLangevin::setup(int vflag)
void FixLangevin::post_force(int vflag)
{
if (tally) post_force_tally();
else post_force_no_tally();
double *rmass = atom->rmass;
// enumerate all 2^6 possibilities for template parameters
// this avoids testing them inside inner loop:
// TSTYLEATOM, GJF, TALLY, BIAS, RMASS, ZERO
#ifdef TEMPLATED_FIX_LANGEVIN
if (tstyle == ATOM)
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,1,1,1,1>();
else post_force_templated<1,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,1,0,1>();
else post_force_templated<1,1,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,1,0,1,1>();
else post_force_templated<1,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,1,0,0,1>();
else post_force_templated<1,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,0,1,1,1>();
else post_force_templated<1,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,1,0,1>();
else post_force_templated<1,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,1,0,0,1,1>();
else post_force_templated<1,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,1,0,0,0,1>();
else post_force_templated<1,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,1,1,1,1>();
else post_force_templated<1,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,1,0,1>();
else post_force_templated<1,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,1,0,1,1>();
else post_force_templated<1,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,1,0,0,1>();
else post_force_templated<1,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,0,0,1,1,1>();
else post_force_templated<1,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,1,0,1>();
else post_force_templated<1,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<1,0,0,0,1,1>();
else post_force_templated<1,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<1,0,0,0,0,1>();
else post_force_templated<1,0,0,0,0,0>();
else
if (gjfflag)
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,1,1,1,1>();
else post_force_templated<0,1,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,1,0,1>();
else post_force_templated<0,1,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,1,0,1,1>();
else post_force_templated<0,1,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,1,0,0,1>();
else post_force_templated<0,1,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,1,0,1,1,1>();
else post_force_templated<0,1,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,1,0,1>();
else post_force_templated<0,1,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,1,0,0,1,1>();
else post_force_templated<0,1,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,1,0,0,0,1>();
else post_force_templated<0,1,0,0,0,0>();
else
if (tallyflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,1,1,1,1>();
else post_force_templated<0,0,1,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,1,0,1>();
else post_force_templated<0,0,1,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,1,0,1,1>();
else post_force_templated<0,0,1,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,1,0,0,1>();
else post_force_templated<0,0,1,0,0,0>();
else
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<0,0,0,1,1,1>();
else post_force_templated<0,0,0,1,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,1,0,1>();
else post_force_templated<0,0,0,1,0,0>();
else
if (rmass)
if (zeroflag) post_force_templated<0,0,0,0,1,1>();
else post_force_templated<0,0,0,0,1,0>();
else
if (zeroflag) post_force_templated<0,0,0,0,0,1>();
else post_force_templated<0,0,0,0,0,0>();
#else
post_force_untemplated(int(tstyle==ATOM), gjfflag, tallyflag,
int(tbiasflag==BIAS), int(rmass!=NULL), zeroflag);
#endif
}
/* ---------------------------------------------------------------------- */
@ -267,9 +441,19 @@ void FixLangevin::post_force_respa(int vflag, int ilevel, int iloop)
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
modify forces using one of the many Langevin styles
------------------------------------------------------------------------- */
void FixLangevin::post_force_no_tally()
#ifdef TEMPLATED_FIX_LANGEVIN
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void FixLangevin::post_force_templated()
#else
void FixLangevin::post_force_untemplated(int Tp_TSTYLEATOM,
int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO)
#endif
{
double gamma1,gamma2;
@ -280,161 +464,133 @@ void FixLangevin::post_force_no_tally()
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
// set current t_target and t_sqrt
// if variable temp, evaluate variable, wrap with clear/add
// reallocate tforce array if necessary
if (tstyle == CONSTANT) {
t_target = t_start + delta * (t_stop-t_start);
tsqrt = sqrt(t_target);
} else {
modify->clearstep_compute();
if (tstyle == EQUAL) {
t_target = input->variable->compute_equal(tvar);
if (t_target < 0.0)
error->one(FLERR,"Fix langevin variable returned negative temperature");
tsqrt = sqrt(t_target);
} else {
if (nlocal > maxatom2) {
maxatom2 = atom->nmax;
memory->destroy(tforce);
memory->create(tforce,maxatom2,"langevin:tforce");
}
input->variable->compute_atom(tvar,igroup,tforce,1,0);
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");
}
modify->addstep_compute(update->ntimestep + 1);
}
// apply damping and thermostat to atoms in group
// for BIAS:
// for Tp_TSTYLEATOM:
// use per-atom per-coord target temperature
// for Tp_GJF:
// use Gronbech-Jensen/Farago algorithm
// else use regular algorithm
// for Tp_TALLY:
// store drag plus random forces in flangevin[nlocal][3]
// for Tp_BIAS:
// calculate temperature since some computes require temp
// computed on current nlocal atoms to remove bias
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
// for ZEROFLAG:
// for Tp_RMASS:
// use per-atom masses
// else use per-type masses
// for Tp_ZERO:
// sum random force over all atoms in group
// subtract sum/count from each atom in group
double fran[3],fsum[3],fsumall[3];
fsum[0] = fsum[1] = fsum[2] = 0.0;
double fdrag[3],fran[3],fsum[3],fsumall[3];
bigint count;
double fswap;
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
if (zeroflag) {
compute_target();
if (Tp_ZERO) {
fsum[0] = fsum[1] = fsum[2] = 0.0;
count = group->count(igroup);
if (count == 0)
error->all(FLERR,"Cannot zero Langevin force of 0 atoms");
}
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
f[i][0] += gamma1*v[i][0] + fran[0];
f[i][1] += gamma1*v[i][1] + fran[1];
f[i][2] += gamma1*v[i][2] + fran[2];
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
}
}
// reallocate flangevin if necessary
} else if (which == BIAS) {
temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
temperature->remove_bias(i,v[i]);
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + fran[0];
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + fran[1];
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + fran[2];
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
temperature->restore_bias(i,v[i]);
}
}
if (Tp_TALLY) {
if (atom->nlocal > maxatom1) {
memory->destroy(flangevin);
maxatom1 = atom->nmax;
memory->create(flangevin,maxatom1,3,"langevin:flangevin");
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
f[i][0] += gamma1*v[i][0] + fran[0];
f[i][1] += gamma1*v[i][1] + fran[1];
f[i][2] += gamma1*v[i][2] + fran[2];
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (Tp_TSTYLEATOM) tsqrt = sqrt(tforce[i]);
if (Tp_RMASS) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
} else {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
}
} else if (which == BIAS) {
temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + fran[0];
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + fran[1];
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + fran[2];
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
temperature->restore_bias(i,v[i]);
}
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
if (Tp_BIAS) {
temperature->remove_bias(i,v[i]);
fdrag[1] = gamma1*v[i][1];
fdrag[2] = gamma1*v[i][2];
fdrag[0] = gamma1*v[i][0];
if (v[i][0] != 0.0) f[i][0] += fdrag[0] + fran[0];
else fran[0] = 0.0;
if (v[i][1] != 0.0) f[i][1] += fdrag[1] + fran[1];
else fran[1] = 0.0;
if (v[i][2] != 0.0) f[i][2] += fdrag[2] + fran[2];
else fran[2] = 0.0;
temperature->restore_bias(i,v[i]);
} else {
fdrag[0] = gamma1*v[i][0];
fdrag[1] = gamma1*v[i][1];
fdrag[2] = gamma1*v[i][2];
}
if (Tp_GJF) {
fswap = 0.5*(fran[0]+franprev[i][0]);
franprev[i][0] = fran[0];
fran[0] = fswap;
fswap = 0.5*(fran[1]+franprev[i][1]);
franprev[i][1] = fran[1];
fran[1] = fswap;
fswap = 0.5*(fran[2]+franprev[i][2]);
franprev[i][2] = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f[i][0] *= gjffac;
f[i][1] *= gjffac;
f[i][2] *= gjffac;
}
f[i][0] += fdrag[0] + fran[0];
f[i][1] += fdrag[1] + fran[1];
f[i][2] += fdrag[2] + fran[2];
if (Tp_TALLY) {
flangevin[i][0] = fdrag[0] + fran[0];
flangevin[i][1] = fdrag[1] + fran[1];
flangevin[i][2] = fdrag[2] + fran[2];
}
if (Tp_ZERO) {
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
}
}
}
// set total force to zero
if (zeroflag) {
if (Tp_ZERO) {
MPI_Allreduce(fsum,fsumall,3,MPI_DOUBLE,MPI_SUM,world);
fsumall[0] /= count;
fsumall[1] /= count;
@ -454,31 +610,18 @@ void FixLangevin::post_force_no_tally()
if (ascale) angmom_thermostat();
}
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
set current t_target and t_sqrt
------------------------------------------------------------------------- */
void FixLangevin::post_force_tally()
void FixLangevin::compute_target()
{
double gamma1,gamma2;
// reallocate flangevin if necessary
if (atom->nlocal > maxatom1) {
memory->destroy(flangevin);
maxatom1 = atom->nmax;
memory->create(flangevin,maxatom1,3,"langevin:flangevin");
}
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double delta = update->ntimestep - update->beginstep;
if (delta != 0.0) delta /= update->endstep - update->beginstep;
// set current t_target and t_sqrt
// if variable temp, evaluate variable, wrap with clear/add
// reallocate tforce array if necessary
@ -507,104 +650,6 @@ void FixLangevin::post_force_tally()
}
modify->addstep_compute(update->ntimestep + 1);
}
// apply damping and thermostat to appropriate atoms
// for BIAS:
// calculate temperature since some computes require temp
// computed on current nlocal atoms to remove bias
// test v = 0 since some computes mask non-participating atoms via v = 0
// and added force has extra term not multiplied by v = 0
double boltz = force->boltz;
double dt = update->dt;
double mvv2e = force->mvv2e;
double ftm2v = force->ftm2v;
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
} else if (which == BIAS) {
temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
temperature->remove_bias(i,v[i]);
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0;
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
} else if (which == BIAS) {
temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (tstyle == ATOM) tsqrt = sqrt(tforce[i]);
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
temperature->remove_bias(i,v[i]);
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
else flangevin[i][0] = 0.0;
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
else flangevin[i][1] = 0.0;
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
else flangevin[i][2] = 0.0;
temperature->restore_bias(i,v[i]);
}
}
}
}
// thermostat omega and angmom
if (oflag) omega_thermostat();
if (ascale) angmom_thermostat();
}
/* ----------------------------------------------------------------------
@ -714,7 +759,7 @@ void FixLangevin::angmom_thermostat()
void FixLangevin::end_of_step()
{
if (!tally) return;
if (!tallyflag) return;
double **v = atom->v;
int *mask = atom->mask;
@ -781,7 +826,7 @@ int FixLangevin::modify_param(int narg, char **arg)
double FixLangevin::compute_scalar()
{
if (!tally || flangevin == NULL) return 0.0;
if (!tallyflag || flangevin == NULL) return 0.0;
// capture the very first energy transfer to thermal reservoir
@ -798,6 +843,8 @@ double FixLangevin::compute_scalar()
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;
@ -825,7 +872,50 @@ void *FixLangevin::extract(const char *str, int &dim)
double FixLangevin::memory_usage()
{
double bytes = 0.0;
if (tally) double bytes = atom->nmax*3 * sizeof(double);
if (tforce) bytes = atom->nmax * sizeof(double);
if (gjfflag) bytes += atom->nmax*3 * sizeof(double);
if (tallyflag) bytes += atom->nmax*3 * sizeof(double);
if (tforce) bytes += atom->nmax * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array for franprev
------------------------------------------------------------------------- */
void FixLangevin::grow_arrays(int nmax)
{
if (!gjfflag) return;
memory->grow(franprev,nmax,3,"fix_langevin:franprev");
array = franprev;
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixLangevin::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nvalues; m++)
array[j][m] = array[i][m];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixLangevin::pack_exchange(int i, double *buf)
{
for (int m = 0; m < nvalues; m++) buf[m] = array[i][m];
return nvalues;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixLangevin::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < nvalues; m++) array[nlocal][m] = buf[m];
return nvalues;
}

22
src/fix_langevin.h
View File

@ -40,15 +40,20 @@ class FixLangevin : public Fix {
virtual double compute_scalar();
double memory_usage();
virtual void *extract(const char *, int &);
void grow_arrays(int);
void copy_arrays(int, int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
protected:
int which,tally,zeroflag,oflag;
int gjfflag,oflag,tallyflag,zeroflag,tbiasflag;
double ascale;
double t_start,t_stop,t_period,t_target;
double *gfactor1,*gfactor2,*ratio;
double energy,energy_onestep;
double tsqrt;
int tstyle,tvar;
double gjffac;
char *tstr;
class AtomVecEllipsoid *avec;
@ -56,6 +61,8 @@ class FixLangevin : public Fix {
int maxatom1,maxatom2;
double **flangevin;
double *tforce;
double **franprev, **array;
int nvalues;
char *id_temp;
class Compute *temperature;
@ -63,10 +70,19 @@ class FixLangevin : public Fix {
int nlevels_respa;
class RanMars *random;
virtual void post_force_no_tally();
virtual void post_force_tally();
// comment next line to turn off templating
#define TEMPLATED_FIX_LANGEVIN
#ifdef TEMPLATED_FIX_LANGEVIN
template < int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO >
void post_force_templated();
#else
void post_force_untemplated(int, int, int,
int, int, int);
#endif
void omega_thermostat();
void angmom_thermostat();
void compute_target();
};
}