From d0e374e4c3a8793c89ee9c34d9df454d20d1b470 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 11 Mar 2008 17:15:30 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1591 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/fix_npt_asphere.cpp | 6 +- src/ASPHERE/fix_nvt_asphere.cpp | 4 +- src/compute.cpp | 14 ++++ src/compute.h | 6 ++ src/compute_temp_com.cpp | 31 +++++--- src/compute_temp_com.h | 1 + src/compute_temp_deform.cpp | 19 +++++ src/compute_temp_deform.h | 1 + src/compute_temp_partial.cpp | 20 +++++ src/compute_temp_partial.h | 1 + src/compute_temp_ramp.cpp | 13 ++++ src/compute_temp_ramp.h | 1 + src/compute_temp_region.cpp | 17 ++++ src/compute_temp_region.h | 1 + src/fix.h | 1 + src/fix_langevin.cpp | 88 ++++++++++++++------- src/fix_langevin.h | 6 +- src/fix_npt.cpp | 133 +++++++++++++++++++++++--------- src/fix_npt.h | 2 +- src/fix_nve.cpp | 3 + src/fix_nvt.cpp | 109 ++++++++++++++++++-------- src/fix_nvt.h | 1 + src/fix_nvt_sllod.cpp | 8 +- src/fix_temp_rescale.cpp | 98 +++++++---------------- src/fix_temp_rescale.h | 7 +- src/style.h | 2 + src/temper.cpp | 25 ++---- 27 files changed, 412 insertions(+), 206 deletions(-) diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index f8d66e21ff..5147e0f65f 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -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; diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index 440dcadf28..ec9b22f238 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -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; diff --git a/src/compute.cpp b/src/compute.cpp index 47e1bbb555..d9a9937c71 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -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]; +} + diff --git a/src/compute.h b/src/compute.h index efe24a4d1f..cde54e3174 100644 --- a/src/compute.h +++ b/src/compute.h @@ -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); diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index 5a251595ae..df80cd3873 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -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]; +} diff --git a/src/compute_temp_com.h b/src/compute_temp_com.h index 02ceae2afc..cdb2b261cf 100644 --- a/src/compute_temp_com.h +++ b/src/compute_temp_com.h @@ -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; diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 2544ac7c56..0413c1a382 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -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]; +} diff --git a/src/compute_temp_deform.h b/src/compute_temp_deform.h index 2eba8ed091..a489e82765 100644 --- a/src/compute_temp_deform.h +++ b/src/compute_temp_deform.h @@ -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; diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index df738e7360..29034a1d8f 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -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; + } +} diff --git a/src/compute_temp_partial.h b/src/compute_temp_partial.h index d12a96d078..b1db6c08cf 100644 --- a/src/compute_temp_partial.h +++ b/src/compute_temp_partial.h @@ -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; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 6c7663472b..e88b446a73 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -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]; +} diff --git a/src/compute_temp_ramp.h b/src/compute_temp_ramp.h index ea3d604017..31b624cbc4 100644 --- a/src/compute_temp_ramp.h +++ b/src/compute_temp_ramp.h @@ -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; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index aeaad13877..3d40dcde6e 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -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; + } +} diff --git a/src/compute_temp_region.h b/src/compute_temp_region.h index 32cf5ac1d7..c8944e2aed 100644 --- a/src/compute_temp_region.h +++ b/src/compute_temp_region.h @@ -25,6 +25,7 @@ class ComputeTempRegion : public Compute { void init(); double compute_scalar(); void compute_vector(); + void remove_bias(int, double *); private: int iregion; diff --git a/src/fix.h b/src/fix.h index a83886639a..729bb40c04 100644 --- a/src/fix.h +++ b/src/fix.h @@ -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;} diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index e52069482e..d13768c816 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -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; +} diff --git a/src/fix_langevin.h b/src/fix_langevin.h index 1deeb84926..b8b9e1b73d 100644 --- a/src/fix_langevin.h +++ b/src/fix_langevin.h @@ -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; }; diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 7636739f6c..0321a81a16 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -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,12 +382,27 @@ void FixNPT::initial_integrate(int vflag) if (igroup == atom->firstgroup) nlocal = atom->nfirst; double dtfm; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - 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]; + + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + 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]; + } + } + + } 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]); + } } } @@ -390,7 +410,7 @@ void FixNPT::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) { @@ -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,12 +446,27 @@ void FixNPT::final_integrate() if (igroup == atom->firstgroup) nlocal = atom->nfirst; double dtfm; - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - 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]; + + if (which == NOBIAS) { + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + 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]; + } + } + + } 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]); + } } } @@ -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,19 +590,33 @@ 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 - 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]; + if (which == NOBIAS) { + 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]; + } + } + + } 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,12 +658,26 @@ void FixNPT::final_integrate_respa(int ilevel) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - 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]; + if (which == NOBIAS) { + 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]; + } + } + + } 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]; diff --git a/src/fix_npt.h b/src/fix_npt.h index f40370711f..5494e0c1ba 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -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; diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index 8035cf25f5..ab05cba17f 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -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(); } diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 76b4bc1caa..be975f4d0e 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -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,15 +174,32 @@ void FixNVT::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - 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]; - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + 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]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + 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,12 +220,26 @@ void FixNVT::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - 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]; + if (which == NOBIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + 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]; + } + } + + } 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]); + } } } @@ -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,19 +296,29 @@ 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]]; + 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]; + } + } - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - 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]; + } 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]); + } } } - // innermost level - also update x of only atoms in NVT group - if (ilevel == 0) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -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; diff --git a/src/fix_nvt.h b/src/fix_nvt.h index ff0ade8965..18d6683370 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -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; diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp index d0fcbdbb37..42bed93bd7 100644 --- a/src/fix_nvt_sllod.cpp +++ b/src/fix_nvt_sllod.cpp @@ -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++) { diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp index 1b352e7086..9c2343244b 100644 --- a/src/fix_temp_rescale.cpp +++ b/src/fix_temp_rescale.cpp @@ -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); - } + newarg[2] = (char *) "temp"; + modify->add_compute(3,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() diff --git a/src/fix_temp_rescale.h b/src/fix_temp_rescale.h index 11137bbbc4..4f2f0c2a93 100644 --- a/src/fix_temp_rescale.h +++ b/src/fix_temp_rescale.h @@ -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; diff --git a/src/style.h b/src/style.h index 278b917f76..31be23a950 100644 --- a/src/style.h +++ b/src/style.h @@ -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) diff --git a/src/temper.cpp b/src/temper.cpp index fea8f30e4c..ee07c3480d 100644 --- a/src/temper.cpp +++ b/src/temper.cpp @@ -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,11 +272,8 @@ 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 // root procs update their value if swap took place