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

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1591 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2008-03-11 17:15:30 +00:00
parent 07c6376bcf
commit d0e374e4c3
27 changed files with 412 additions and 206 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/ASPHERE/fix_npt_asphere.cpp
View File

@ -79,7 +79,7 @@ void FixNPTASphere::initial_integrate(int vflag)
}
ang_factor = exp(-dthalf*eta_dot);
// v update only for atoms in NPT group
// v update only for atoms in group
double **x = atom->x;
double **v = atom->v;
@ -107,7 +107,7 @@ void FixNPTASphere::initial_integrate(int vflag)
remap(0);
// x update by full step only for atoms in NPT group
// x update by full step only for atoms in group
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -147,7 +147,7 @@ void FixNPTASphere::final_integrate()
{
int i;
// v update only for atoms in NPT group
// v update only for atoms in group
double **v = atom->v;
double **f = atom->f;

4
src/ASPHERE/fix_nvt_asphere.cpp
View File

@ -61,7 +61,7 @@ void FixNVTASphere::initial_integrate(int vflag)
eta += dtv*eta_dot;
factor = exp(-dthalf*eta_dot);
// update v and x of only atoms in NVT group
// update v and x of only atoms in group
double **x = atom->x;
double **v = atom->v;
@ -106,7 +106,7 @@ void FixNVTASphere::final_integrate()
{
double dtfm;
// update v of only atoms in NVT group
// update v of only atoms in group
double **v = atom->v;
double **f = atom->f;

14
src/compute.cpp
View File

@ -60,6 +60,7 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
scalar_flag = vector_flag = peratom_flag = 0;
tempflag = pressflag = peflag = 0;
pressatomflag = peatomflag = 0;
tempbias = 0;
id_pre = NULL;
timeflag = 0;
invoked = 0;
@ -162,3 +163,16 @@ int Compute::matchstep(int ntimestep)
}
return 0;
}
/* ----------------------------------------------------------------------
add back in velocity bias removed by remove_bias() to leave thermal temp
assume remove_bias() was previously called for this atom
------------------------------------------------------------------------- */
void Compute::restore_bias(double *v)
{
v[0] += vbias[0];
v[1] += vbias[1];
v[2] += vbias[2];
}

6
src/compute.h
View File

@ -45,6 +45,9 @@ class Compute : protected Pointers {
int peflag; // 1 if Compute calculates PE (uses Force energies)
int peatomflag; // 1 if Compute calculates per-atom PE
int tempbias; // 0/1 if has bias routines for isolating thermal temp
double vbias[3]; // stored velocity bias for one atom
char *id_pre; // ID of pre-compute the Compute may store
int timeflag; // 1 if Compute stores list of timesteps it's called on
@ -73,6 +76,9 @@ class Compute : protected Pointers {
virtual int pack_reverse_comm(int, int, double *) {return 0;}
virtual void unpack_reverse_comm(int, int *, double *) {}
virtual void remove_bias(int, double *) {}
void restore_bias(double *);
void addstep(int);
int matchstep(int);

31
src/compute_temp_com.cpp
View File

@ -44,6 +44,8 @@ ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
vbias[0] = vbias[1] = vbias[2] = 0.0;
vector = new double[6];
}
@ -81,12 +83,12 @@ void ComputeTempCOM::recount()
double ComputeTempCOM::compute_scalar()
{
double vcm[3],vthermal[3];
double vthermal[3];
invoked |= INVOKED_SCALAR;
if (dynamic) masstotal = group->mass(igroup);
group->vcm(igroup,masstotal,vcm);
group->vcm(igroup,masstotal,vbias);
double **v = atom->v;
double *mass = atom->mass;
@ -98,9 +100,9 @@ double ComputeTempCOM::compute_scalar()
double t = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
vthermal[0] = v[i][0] - vcm[0];
vthermal[1] = v[i][1] - vcm[1];
vthermal[2] = v[i][2] - vcm[2];
vthermal[0] = v[i][0] - vbias[0];
vthermal[1] = v[i][1] - vbias[1];
vthermal[2] = v[i][2] - vbias[2];
if (mass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
@ -120,12 +122,12 @@ double ComputeTempCOM::compute_scalar()
void ComputeTempCOM::compute_vector()
{
int i;
double vcm[3],vthermal[3];
double vthermal[3];
invoked |= INVOKED_VECTOR;
if (dynamic) masstotal = group->mass(igroup);
group->vcm(igroup,masstotal,vcm);
group->vcm(igroup,masstotal,vbias);
double **x = atom->x;
double **v = atom->v;
@ -140,9 +142,9 @@ void ComputeTempCOM::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
vthermal[0] = v[i][0] - vcm[0];
vthermal[1] = v[i][1] - vcm[1];
vthermal[2] = v[i][2] - vcm[2];
vthermal[0] = v[i][0] - vbias[0];
vthermal[1] = v[i][1] - vbias[1];
vthermal[2] = v[i][2] - vbias[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
@ -157,3 +159,12 @@ void ComputeTempCOM::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ---------------------------------------------------------------------- */
void ComputeTempCOM::remove_bias(int i, double *v)
{
v[0] = v[0] - vbias[0];
v[1] = v[1] - vbias[1];
v[2] = v[2] - vbias[2];
}

1
src/compute_temp_com.h
View File

@ -25,6 +25,7 @@ class ComputeTempCOM : public Compute {
void init();
double compute_scalar();
void compute_vector();
void remove_bias(int, double *);
private:
int fix_dof;

19
src/compute_temp_deform.cpp
View File

@ -46,6 +46,7 @@ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
vector = new double[6];
}
@ -198,3 +199,21 @@ void ComputeTempDeform::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ---------------------------------------------------------------------- */
void ComputeTempDeform::remove_bias(int i, double *v)
{
double lamda[3];
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
domain->x2lamda(atom->x[i],lamda);
vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] +
h_rate[4]*lamda[2] + h_ratelo[0];
vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1];
vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2];
v[0] -= vbias[0];
v[1] -= vbias[1];
v[2] -= vbias[2];
}

1
src/compute_temp_deform.h
View File

@ -25,6 +25,7 @@ class ComputeTempDeform : public Compute {
void init();
double compute_scalar();
void compute_vector();
void remove_bias(int, double *);
private:
int fix_dof;

20
src/compute_temp_partial.cpp
View File

@ -42,6 +42,8 @@ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
vbias[0] = vbias[1] = vbias[2] = 0.0;
vector = new double[6];
}
@ -140,3 +142,21 @@ void ComputeTempPartial::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ---------------------------------------------------------------------- */
void ComputeTempPartial::remove_bias(int i, double *v)
{
if (!xflag) {
vbias[0] = v[0];
v[0] = 0.0;
}
if (!yflag) {
vbias[1] = v[1];
v[1] = 0.0;
}
if (!zflag) {
vbias[2] = v[2];
v[2] = 0.0;
}
}

1
src/compute_temp_partial.h
View File

@ -25,6 +25,7 @@ class ComputeTempPartial : public Compute {
void init();
double compute_scalar();
void compute_vector();
void remove_bias(int, double *);
private:
int xflag,yflag,zflag;

13
src/compute_temp_ramp.cpp
View File

@ -107,6 +107,8 @@ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
vbias[0] = vbias[1] = vbias[2] = 0.0;
vector = new double[6];
}
@ -224,3 +226,14 @@ void ComputeTempRamp::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ---------------------------------------------------------------------- */
void ComputeTempRamp::remove_bias(int i, double *v)
{
double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
fraction = MAX(fraction,0.0);
fraction = MIN(fraction,1.0);
vbias[v_dim] = v_lo + fraction*(v_hi - v_lo);
v[v_dim] -= vbias[v_dim];
}

1
src/compute_temp_ramp.h
View File

@ -25,6 +25,7 @@ class ComputeTempRamp : public Compute {
void init();
double compute_scalar();
void compute_vector();
void remove_bias(int, double *);
private:
int coord_dim;

17
src/compute_temp_region.cpp
View File

@ -40,6 +40,7 @@ ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) :
extscalar = 0;
extvector = 1;
tempflag = 1;
tempbias = 1;
vector = new double[6];
}
@ -137,3 +138,19 @@ void ComputeTempRegion::compute_vector()
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
}
/* ---------------------------------------------------------------------- */
void ComputeTempRegion::remove_bias(int i, double *v)
{
double *x = atom->x[i];
if (atom->mask[i] & groupbit &&
domain->regions[iregion]->match(x[0],x[1],x[2]))
vbias[0] = vbias[1] = vbias[2] = 0.0;
else {
vbias[0] = v[0];
vbias[1] = v[1];
vbias[2] = v[2];
v[0] = v[1] = v[2] = 0.0;
}
}

1
src/compute_temp_region.h
View File

@ -25,6 +25,7 @@ class ComputeTempRegion : public Compute {
void init();
double compute_scalar();
void compute_vector();
void remove_bias(int, double *);
private:
int iregion;

1
src/fix.h
View File

@ -106,6 +106,7 @@ class Fix : protected Pointers {
virtual int dof(int) {return 0;}
virtual void deform(int) {}
virtual void reset_target(double) {}
virtual void reset_dt() {}
virtual int modify_param(int, char **) {return 0;}

88
src/fix_langevin.cpp
View File

@ -19,6 +19,8 @@
#include "atom.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
@ -28,6 +30,8 @@
using namespace LAMMPS_NS;
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
@ -55,19 +59,11 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
// optional args
flagx = flagy = flagz = 1;
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
iregion = -1;
int iarg = 7;
while (iarg < narg) {
if (strcmp(arg[iarg],"axes") == 0) {
if (iarg+4 > narg) error->all("Illegal fix langevin command");
flagx = atoi(arg[iarg+1]);
flagy = atoi(arg[iarg+2]);
flagz = atoi(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"scale") == 0) {
if (strcmp(arg[iarg],"scale") == 0) {
if (iarg+3 > narg) error->all("Illegal fix langevin command");
int itype = atoi(arg[iarg+1]);
double scale = atof(arg[iarg+2]);
@ -75,13 +71,13 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
error->all("Illegal fix langevin command");
ratio[itype] = scale;
iarg += 3;
} else if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all("Illegal fix langevin command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all("Fix langevin region ID does not exist");
iarg += 2;
} else error->all("Illegal fix langevin command");
}
// set temperature = NULL, user can override via fix_modify if wants bias
id_temp = NULL;
temperature = NULL;
}
/* ---------------------------------------------------------------------- */
@ -92,6 +88,7 @@ FixLangevin::~FixLangevin()
delete [] gfactor1;
delete [] gfactor2;
delete [] ratio;
delete [] id_temp;
}
/* ---------------------------------------------------------------------- */
@ -121,6 +118,9 @@ void FixLangevin::init()
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
if (temperature && temperature->tempbias) which = BIAS;
else which = NOBIAS;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
@ -142,6 +142,8 @@ void FixLangevin::setup(int vflag)
void FixLangevin::post_force(int vflag)
{
double gamma1,gamma2;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
@ -154,34 +156,40 @@ void FixLangevin::post_force(int vflag)
double t_target = t_start + delta * (t_stop-t_start);
double tsqrt = sqrt(t_target);
double gamma1,gamma2;
// apply damping and thermostat to appropriate atoms
// apply damping and thermostat to all atoms in fix group
if (iregion == -1) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
if (flagx) f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (flagy) f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (flagz) f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
}
}
// apply damping and thermostat to all atoms in fix group and in region
// invoke temperature since some computes require it to remove bias
// test on v = 0 since some computes mask non-participating atoms via v = 0
} else if (which == BIAS) {
double tmp = temperature->compute_scalar();
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit &&
domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) {
if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]];
gamma2 = gfactor2[type[i]] * tsqrt;
if (flagx) f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (flagy) f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (flagz) f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
temperature->remove_bias(i,v[i]);
if (v[i][0] != 0.0)
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
if (v[i][1] != 0.0)
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
if (v[i][2] != 0.0)
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
temperature->restore_bias(v[i]);
}
}
}
}
@ -209,3 +217,27 @@ void FixLangevin::reset_dt()
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
}
/* ---------------------------------------------------------------------- */
int FixLangevin::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Could not find fix_modify temp ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temp ID does not compute temperature");
if (temperature->igroup != igroup && comm->me == 0)
error->warning("Group for fix_modify temp != fix group");
return 2;
}
return 0;
}

6
src/fix_langevin.h
View File

@ -29,12 +29,16 @@ class FixLangevin : public Fix {
void post_force_respa(int, int, int);
void reset_target(double);
void reset_dt();
int modify_param(int, char **);
private:
int which;
double t_start,t_stop,t_period;
int flagx,flagy,flagz,iregion;
double *gfactor1,*gfactor2,*ratio;
char *id_temp;
class Compute *temperature;
int nlevels_respa;
class RanMars *random;
};

85
src/fix_npt.cpp
View File

@ -36,6 +36,8 @@ using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
@ -257,6 +259,9 @@ void FixNPT::init()
if (icompute < 0) error->all("Temp ID for fix npt does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Press ID for fix npt does not exist");
pressure = modify->compute[icompute];
@ -365,7 +370,7 @@ void FixNPT::initial_integrate(int vflag)
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPT group
// v update only for atoms in group
double **x = atom->x;
double **v = atom->v;
@ -377,6 +382,8 @@ void FixNPT::initial_integrate(int vflag)
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double dtfm;
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -386,11 +393,24 @@ void FixNPT::initial_integrate(int vflag)
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
temperature->restore_bias(v[i]);
}
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
// x update by full step only for atoms in NPT group
// x update by full step only for atoms in group
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -415,7 +435,7 @@ void FixNPT::final_integrate()
{
int i;
// v update only for atoms in NPT group
// v update only for atoms in group
double **v = atom->v;
double **f = atom->f;
@ -426,6 +446,8 @@ void FixNPT::final_integrate()
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double dtfm;
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -435,6 +457,19 @@ void FixNPT::final_integrate()
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
temperature->restore_bias(v[i]);
}
}
}
// compute new T,P
t_current = temperature->compute_scalar();
@ -505,7 +540,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
// outermost level - update eta_dot and omega_dot, apply to v, remap box
// all other levels - NVE update of v
// x,v updates only performed for atoms in NPT group
// x,v updates only performed for atoms in group
if (ilevel == nlevels_respa-1) {
@ -538,7 +573,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPT group
// v update only for atoms in group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -555,8 +590,9 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
} else {
// v update only for atoms in NPT group
// v update only for atoms in group
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -565,9 +601,22 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
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];
temperature->restore_bias(v[i]);
}
}
}
}
// innermost level - also update x only for atoms in NPT group
// innermost level - also update x only for atoms in group
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {
@ -594,12 +643,12 @@ void FixNPT::final_integrate_respa(int ilevel)
// outermost level - update eta_dot and omega_dot,
// apply to v via final_integrate()
// all other levels - NVE update of v
// v update only performed for atoms in NPT group
// v update only performed for atoms in group
if (ilevel == nlevels_respa-1) final_integrate();
else {
// v update only for atoms in NPT group
// v update only for atoms in group
double **v = atom->v;
double **f = atom->f;
@ -609,6 +658,7 @@ void FixNPT::final_integrate_respa(int ilevel)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -617,6 +667,19 @@ void FixNPT::final_integrate_respa(int ilevel)
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
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];
temperature->restore_bias(v[i]);
}
}
}
}
}
@ -763,7 +826,7 @@ int FixNPT::modify_param(int narg, char **arg)
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify temp ID");
temperature = modify->compute[icompute];
@ -793,7 +856,7 @@ int FixNPT::modify_param(int narg, char **arg)
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(id_press);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify press ID");
pressure = modify->compute[icompute];

2
src/fix_npt.h
View File

@ -36,7 +36,7 @@ class FixNPT : public Fix {
void reset_dt();
protected:
int dimension;
int dimension,which;
double dtv,dtf,dtq,dthalf;
double boltz,nktv2p;
double vol0;

3
src/fix_nve.cpp
View File

@ -146,6 +146,9 @@ void FixNVE::initial_integrate_respa(int vflag, int ilevel, int flag)
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
// innermost level - NVE update of v and x
// all other levels - NVE update of v
if (ilevel == 0) initial_integrate(vflag);
else final_integrate();
}

67
src/fix_nvt.cpp
View File

@ -31,6 +31,8 @@
using namespace LAMMPS_NS;
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
@ -117,6 +119,9 @@ void FixNVT::init()
if (icompute < 0) error->all("Temp ID for fix nvt does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
// set timesteps and frequencies
dtv = update->dt;
@ -158,7 +163,7 @@ void FixNVT::initial_integrate(int vflag)
eta += dtv*eta_dot;
factor = exp(-dthalf*eta_dot);
// update v and x of only atoms in NVT group
// update v and x of only atoms in group
double **x = atom->x;
double **v = atom->v;
@ -169,6 +174,7 @@ void FixNVT::initial_integrate(int vflag)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -180,6 +186,22 @@ void FixNVT::initial_integrate(int vflag)
x[i][2] += dtv * v[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(v[i]);
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
@ -188,7 +210,7 @@ void FixNVT::final_integrate()
{
double dtfm;
// update v of only atoms in NVT group
// update v of only atoms in group
double **v = atom->v;
double **f = atom->f;
@ -198,6 +220,7 @@ void FixNVT::final_integrate()
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]] * factor;
@ -207,6 +230,19 @@ void FixNVT::final_integrate()
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]] * factor;
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(v[i]);
}
}
}
// compute current T
t_current = temperature->compute_scalar();
@ -242,8 +278,9 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// outermost level - update eta_dot and apply to v
// all other levels - NVE update of v
// outermost level - update eta_dot and apply to v with factor
// all other levels - NVE update of v (factor = 1)
// innermost level - also update x
if (ilevel == nlevels_respa-1) {
double delta = update->ntimestep - update->beginstep;
@ -259,8 +296,7 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
factor = exp(-dthalf*eta_dot);
} else factor = 1.0;
// update v of only atoms in NVT group
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
@ -270,7 +306,18 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
}
}
// innermost level - also update x of only atoms in NVT group
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(v[i]);
}
}
}
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {
@ -287,10 +334,9 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
void FixNVT::final_integrate_respa(int ilevel)
{
double dtfm;
// set timesteps by level
double dtfm;
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
@ -299,9 +345,6 @@ void FixNVT::final_integrate_respa(int ilevel)
if (ilevel == nlevels_respa-1) final_integrate();
else {
// update v of only atoms in NVT group
double **v = atom->v;
double **f = atom->f;
double *mass = atom->mass;

1
src/fix_nvt.h
View File

@ -37,6 +37,7 @@ class FixNVT : public Fix {
void reset_dt();
protected:
int which;
double t_start,t_stop;
double t_current,t_target;
double t_freq,drag,drag_factor;

8
src/fix_nvt_sllod.cpp
View File

@ -78,7 +78,7 @@ void FixNVTSlodd::initial_integrate(int vflag)
eta += dtv*eta_dot;
factor = exp(-dthalf*eta_dot);
// update vthermal and x of only atoms in NVT group
// update vthermal and x of only atoms in group
// lamda = 0-1 triclinic lamda coords
// vstream = streaming velocity = Hrate*lamda + Hratelo
// vthermal = thermal velocity = v - vstream
@ -134,7 +134,7 @@ void FixNVTSlodd::final_integrate()
{
double dtfm;
// update vthermal of only atoms in NVT group
// update vthermal of only atoms in group
// lamda = 0-1 triclinic lamda coords
// vstream = streaming velocity = Hrate*lamda + Hratelo
// vthermal = thermal velocity = v - vstream
@ -232,7 +232,7 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
factor = exp(-dthalf*eta_dot);
} else factor = 1.0;
// update v of only atoms in NVT group
// update v of only atoms in group
double *h_rate = domain->h_rate;
double *h_ratelo = domain->h_ratelo;
@ -265,7 +265,7 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
}
}
// innermost level - also update x of only atoms in NVT group
// innermost level - also update x of only atoms in group
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {

94
src/fix_temp_rescale.cpp
View File

@ -28,7 +28,7 @@
using namespace LAMMPS_NS;
enum{STANDARD,REGION,PARTIAL};
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
@ -45,41 +45,11 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
t_start = atof(arg[4]);
t_end = atof(arg[5]);
t_stop = atof(arg[5]);
t_window = atof(arg[6]);
fraction = atof(arg[7]);
// optional args
iregion = -1;
partial = 0;
xflag = yflag = zflag = 1;
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all("Illegal fix temp/rescale command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all("Fix temp/rescale region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"partial") == 0) {
if (iarg+4 > narg) error->all("Illegal fix temp/rescale command");
xflag = atoi(arg[iarg+1]);
yflag = atoi(arg[iarg+2]);
zflag = atoi(arg[iarg+3]);
partial = 1;
iarg += 4;
} else error->all("Illegal fix temp/rescale command");
}
if (iregion == -1 && partial == 0) type = STANDARD;
else if (iregion >= 0 && partial == 0) type = REGION;
else if (iregion == -1 && partial) type = PARTIAL;
else
error->all("Cannot use both region, partial options in fix temp/rescale");
// create a new compute temp or temp/region or temp/partial
// create a new compute temp
// id = fix-ID + temp, compute group = fix group
int n = strlen(id) + 6;
@ -90,23 +60,8 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
char **newarg = new char*[6];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
if (type == STANDARD) {
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
} else if (type == REGION) {
newarg[2] = (char *) "temp/region";
newarg[3] = domain->regions[iregion]->id;
modify->add_compute(4,newarg);
} else if (type == PARTIAL) {
newarg[2] = (char *) "temp/partial";
if (xflag) newarg[3] = (char *) "1";
else newarg[3] = (char *) "0";
if (yflag) newarg[4] = (char *) "1";
else newarg[4] = (char *) "0";
if (zflag) newarg[5] = (char *) "1";
else newarg[5] = (char *) "0";
modify->add_compute(6,newarg);
}
delete [] newarg;
tflag = 1;
@ -141,6 +96,9 @@ void FixTempRescale::init()
if (icompute < 0) error->all("Temp ID for fix temp/rescale does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
temperature->init(); // not yet called by Modify::init()
efactor = (0.5 * force->boltz * temperature->dof);
energy = 0.0;
@ -151,20 +109,24 @@ void FixTempRescale::init()
void FixTempRescale::end_of_step()
{
double t_current = temperature->compute_scalar();
if (t_current == 0.0)
error->all("Computed temperature for fix temp/berendsen cannot be 0.0");
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
double t_target = t_start + delta * (t_end-t_start);
double t_target = t_start + delta * (t_stop-t_start);
// rescale velocity of appropriate atoms if outside window
if (fabs(t_current-t_target) > t_window) {
t_target = t_current - fraction*(t_current-t_target);
double factor = sqrt(t_target/t_current);
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (type == STANDARD) {
if (which == NOBIAS) {
energy += (t_current-t_target) * efactor;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -174,25 +136,15 @@ void FixTempRescale::end_of_step()
}
}
} else if (type == REGION) {
efactor = (0.5 * force->boltz * temperature->dof);
energy += (t_current-t_target) * efactor;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit &&
domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) {
v[i][0] *= factor;
v[i][1] *= factor;
v[i][2] *= factor;
}
}
} else {
} else if (which == BIAS) {
energy += (t_current-t_target) * efactor;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (xflag) v[i][0] *= factor;
if (yflag) v[i][1] *= factor;
if (zflag) v[i][2] *= factor;
temperature->remove_bias(i,v[i]);
v[i][0] *= factor;
v[i][1] *= factor;
v[i][2] *= factor;
temperature->restore_bias(v[i]);
}
}
}
@ -228,6 +180,14 @@ int FixTempRescale::modify_param(int narg, char **arg)
return 0;
}
/* ---------------------------------------------------------------------- */
void FixTempRescale::reset_target(double t_new)
{
t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
double FixTempRescale::compute_scalar()

7
src/fix_temp_rescale.h
View File

@ -25,12 +25,13 @@ class FixTempRescale : public Fix {
int setmask();
void init();
void end_of_step();
double compute_scalar();
int modify_param(int, char **);
void reset_target(double);
double compute_scalar();
private:
int type,iregion,partial,xflag,yflag,zflag;
double t_start,t_end,t_window;
int which;
double t_start,t_stop,t_window;
double fraction,energy,efactor;
char *id_temp;

2
src/style.h
View File

@ -180,6 +180,7 @@ DumpStyle(xyz,DumpXYZ)
#include "fix_spring.h"
#include "fix_spring_rg.h"
#include "fix_spring_self.h"
#include "fix_temp_berendsen.h"
#include "fix_temp_rescale.h"
#include "fix_tmd.h"
#include "fix_viscosity.h"
@ -233,6 +234,7 @@ FixStyle(SHEAR_HISTORY,FixShearHistory)
FixStyle(spring,FixSpring)
FixStyle(spring/rg,FixSpringRG)
FixStyle(spring/self,FixSpringSelf)
FixStyle(temp/berendsen,FixTempBerendsen)
FixStyle(temp/rescale,FixTempRescale)
FixStyle(tmd,FixTMD)
FixStyle(viscosity,FixViscosity)

23
src/temper.cpp
View File

@ -30,8 +30,6 @@
#include "output.h"
#include "thermo.h"
#include "fix.h"
#include "fix_nvt.h"
#include "fix_langevin.h"
#include "random_park.h"
#include "finish.h"
#include "timer.h"
@ -40,8 +38,6 @@
using namespace LAMMPS_NS;
#define NVT 1
#define LANGEVIN 2
// #define TEMPER_DEBUG 1
/* ---------------------------------------------------------------------- */
@ -97,10 +93,11 @@ void Temper::command(int narg, char **arg)
// fix style must be appropriate for temperature control
if (strcmp(modify->fix[whichfix]->style,"nvt") == 0) fixstyle = NVT;
else if (strcmp(modify->fix[whichfix]->style,"langevin") == 0)
fixstyle = LANGEVIN;
else error->universe_all("Tempering fix is not valid");
if ((strcmp(modify->fix[whichfix]->style,"nvt") != 0) &&
(strcmp(modify->fix[whichfix]->style,"langevin") != 0) &&
(strcmp(modify->fix[whichfix]->style,"temp/berendsen") != 0) &&
(strcmp(modify->fix[whichfix]->style,"temp/rescale") != 0))
error->universe_all("Tempering temperature fix is not valid");
// setup for long tempering run
@ -173,10 +170,7 @@ void Temper::command(int narg, char **arg)
if (narg == 7) {
double new_temp = set_temp[my_set_temp];
if (fixstyle == NVT)
((FixNVT *) modify->fix[whichfix])->reset_target(new_temp);
else if (fixstyle == LANGEVIN)
((FixLangevin *) modify->fix[whichfix])->reset_target(new_temp);
modify->fix[whichfix]->reset_target(new_temp);
}
// setup tempering runs
@ -278,10 +272,7 @@ void Temper::command(int narg, char **arg)
if (swap) {
new_temp = set_temp[partner_set_temp];
if (fixstyle == NVT)
((FixNVT *) modify->fix[whichfix])->reset_target(new_temp);
else if (fixstyle == LANGEVIN)
((FixLangevin *) modify->fix[whichfix])->reset_target(new_temp);
modify->fix[whichfix]->reset_target(new_temp);
}
// update my_set_temp and temp2world on every proc