diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index 5147e0f65f..a5f2752bcc 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -27,9 +27,12 @@ #include "kspace.h" #include "update.h" #include "domain.h" +#include "error.h" using namespace LAMMPS_NS; +enum{NOBIAS,BIAS}; + /* ---------------------------------------------------------------------- */ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : @@ -41,6 +44,14 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : "quat, angmom, torque, shape"); } +/* ---------------------------------------------------------------------- */ + +void FixNPTASphere::init() +{ + FixNPT::init(); + dtq = 0.5 * update->dt; +} + /* ---------------------------------------------------------------------- 1st half of Verlet update ------------------------------------------------------------------------- */ @@ -48,6 +59,7 @@ FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) : void FixNPTASphere::initial_integrate(int vflag) { int i; + double dtfm; double delta = update->ntimestep - update->beginstep; delta /= update->endstep - update->beginstep; @@ -77,9 +89,9 @@ void FixNPTASphere::initial_integrate(int vflag) factor[i] = exp(-dthalf*(eta_dot+omega_dot[i])); dilation[i] = exp(dthalf*omega_dot[i]); } - ang_factor = exp(-dthalf*eta_dot); + factor_rotate = exp(-dthalf*eta_dot); - // v update only for atoms in group + // update v of only atoms in group double **x = atom->x; double **v = atom->v; @@ -93,13 +105,25 @@ void FixNPTASphere::initial_integrate(int vflag) int nlocal = atom->nlocal; 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]); + } } } @@ -120,11 +144,12 @@ void FixNPTASphere::initial_integrate(int vflag) // update angular momentum by 1/2 step // update quaternion a full step via Richardson iteration // returns new normalized quaternion + for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - angmom[i][0] = angmom[i][0] * ang_factor + dtf * torque[i][0]; - angmom[i][1] = angmom[i][1] * ang_factor + dtf * torque[i][1]; - angmom[i][2] = angmom[i][2] * ang_factor + dtf * torque[i][2]; + angmom[i][0] = angmom[i][0]*factor_rotate + dtf*torque[i][0]; + angmom[i][1] = angmom[i][1]*factor_rotate + dtf*torque[i][1]; + angmom[i][2] = angmom[i][2]*factor_rotate + dtf*torque[i][2]; double inertia[3]; calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); @@ -146,8 +171,9 @@ void FixNPTASphere::initial_integrate(int vflag) void FixNPTASphere::final_integrate() { int i; + double dtfm; - // v update only for atoms in group + // update v of only atoms in group double **v = atom->v; double **f = atom->f; @@ -159,16 +185,31 @@ void FixNPTASphere::final_integrate() int nlocal = atom->nlocal; 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]; - angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * ang_factor; - angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * ang_factor; - angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * ang_factor; + 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]; + angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * factor_rotate; + angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * factor_rotate; + angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate; + } + } + } 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]); + angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * factor_rotate; + angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * factor_rotate; + angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate; + } } } @@ -208,6 +249,14 @@ void FixNPTASphere::final_integrate() } } +/* ---------------------------------------------------------------------- */ + +void FixNPTASphere::reset_dt() +{ + FixNPT::reset_dt(); + dtq = 0.5 * update->dt; +} + /* ---------------------------------------------------------------------- Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ diff --git a/src/ASPHERE/fix_npt_asphere.h b/src/ASPHERE/fix_npt_asphere.h index 8cad21b0bc..e7e19b07f4 100755 --- a/src/ASPHERE/fix_npt_asphere.h +++ b/src/ASPHERE/fix_npt_asphere.h @@ -21,12 +21,14 @@ namespace LAMMPS_NS { class FixNPTASphere : public FixNPT { public: FixNPTASphere(class LAMMPS *, int, char **); - ~FixNPTASphere() {} + void init(); void initial_integrate(int); void final_integrate(); + void reset_dt(); private: - double ang_factor; + double dtq; + double factor_rotate; void richardson(double *, double *, double *); void omega_from_mq(double *, double *, double *, double *); diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index ec9b22f238..2b29a7df0a 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -32,6 +32,8 @@ using namespace LAMMPS_NS; +enum{NOBIAS,BIAS}; + /* ---------------------------------------------------------------------- */ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : @@ -45,6 +47,14 @@ FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ +void FixNVTASphere::init() +{ + FixNVT::init(); + dtq = 0.5 * update->dt; +} + +/* ---------------------------------------------------------------------- */ + void FixNVTASphere::initial_integrate(int vflag) { double dtfm; @@ -75,27 +85,56 @@ void FixNVTASphere::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]; + + // update angular momentum by 1/2 step + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + + angmom[i][0] = angmom[i][0]*factor + dtf*torque[i][0]; + angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1]; + angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2]; + + double inertia[3]; + calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); + richardson(quat[i],angmom[i],inertia); + } + } - // update angular momentum by 1/2 step - // update quaternion a full step via Richardson iteration - // returns new normalized quaternion - - angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0]; - angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1]; - angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2]; - - double inertia[3]; - calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); - richardson(quat[i],angmom[i],inertia); + } 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]; + + // update angular momentum by 1/2 step + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + + angmom[i][0] = angmom[i][0]*factor + dtf*torque[i][0]; + angmom[i][1] = angmom[i][1]*factor + dtf*torque[i][1]; + angmom[i][2] = angmom[i][2]*factor + dtf*torque[i][2]; + + double inertia[3]; + calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia); + richardson(quat[i],angmom[i],inertia); + } } } } @@ -118,16 +157,32 @@ void FixNVTASphere::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]; + angmom[i][0] = (angmom[i][0] + dtf*torque[i][0]) * factor; + angmom[i][1] = (angmom[i][1] + dtf*torque[i][1]) * factor; + angmom[i][2] = (angmom[i][2] + dtf*torque[i][2]) * factor; + } + } - angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0]; - angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1]; - angmom[i][2] = angmom[i][2] * factor + dtf * torque[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]); + angmom[i][0] = (angmom[i][0] + dtf*torque[i][0]) * factor; + angmom[i][1] = (angmom[i][1] + dtf*torque[i][1]) * factor; + angmom[i][2] = (angmom[i][2] + dtf*torque[i][2]) * factor; + } } } @@ -142,6 +197,14 @@ void FixNVTASphere::final_integrate() eta_dot *= drag_factor; } +/* ---------------------------------------------------------------------- */ + +void FixNVTASphere::reset_dt() +{ + FixNVT::reset_dt(); + dtq = 0.5 * update->dt; +} + /* ---------------------------------------------------------------------- Richardson iteration to update quaternion accurately ------------------------------------------------------------------------- */ diff --git a/src/ASPHERE/fix_nvt_asphere.h b/src/ASPHERE/fix_nvt_asphere.h index 8f037d541b..5dd4de37b5 100755 --- a/src/ASPHERE/fix_nvt_asphere.h +++ b/src/ASPHERE/fix_nvt_asphere.h @@ -21,11 +21,14 @@ namespace LAMMPS_NS { class FixNVTASphere : public FixNVT { public: FixNVTASphere(class LAMMPS *, int, char **); - ~FixNVTASphere() {} + void init(); void initial_integrate(int); void final_integrate(); + void reset_dt(); private: + double dtq; + void richardson(double *, double *, double *); void omega_from_mq(double *, double *, double *, double *); void calculate_inertia(double mass, double *shape, double *inertia); diff --git a/src/compute.cpp b/src/compute.cpp index d9a9937c71..83edb2d933 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -175,4 +175,3 @@ void Compute::restore_bias(double *v) v[1] += vbias[1]; v[2] += vbias[2]; } - diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index 0321a81a16..a377fb6dae 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -178,6 +178,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : if (strcmp(style,"npt") == 0) newarg[2] = (char *) "temp"; else if (strcmp(style,"npt/asphere") == 0) newarg[2] = (char *) "temp/asphere"; + else if (strcmp(style,"npt/sphere") == 0) + newarg[2] = (char *) "temp/sphere"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; @@ -270,7 +272,6 @@ void FixNPT::init() dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; dthalf = 0.5 * update->dt; double freq = MAX(p_freq[0],p_freq[1]); @@ -392,7 +393,6 @@ void FixNPT::initial_integrate(int vflag) 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) { @@ -456,7 +456,6 @@ void FixNPT::final_integrate() 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) { @@ -601,7 +600,6 @@ 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) { @@ -667,7 +665,6 @@ 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) { @@ -894,7 +891,6 @@ void FixNPT::reset_dt() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; dthalf = 0.5 * update->dt; double freq = MAX(p_freq[0],p_freq[1]); diff --git a/src/fix_npt.h b/src/fix_npt.h index 5494e0c1ba..c770e9fd16 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -23,7 +23,7 @@ class FixNPT : public Fix { FixNPT(class LAMMPS *, int, char **); virtual ~FixNPT(); int setmask(); - void init(); + virtual void init(); void setup(int); void initial_integrate(int); virtual void final_integrate(); @@ -33,11 +33,11 @@ class FixNPT : public Fix { void write_restart(FILE *); void restart(char *); int modify_param(int, char **); - void reset_dt(); + virtual void reset_dt(); protected: int dimension,which; - double dtv,dtf,dtq,dthalf; + double dtv,dtf,dthalf; double boltz,nktv2p; double vol0; diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index be975f4d0e..1ea854955e 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -80,6 +80,8 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) : newarg[2] = (char *) "temp/asphere"; else if (strcmp(style,"nvt/sllod") == 0) newarg[2] = (char *) "temp/deform"; + else if (strcmp(style,"nvt/sphere") == 0) + newarg[2] = (char *) "temp/sphere"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; @@ -126,7 +128,6 @@ void FixNVT::init() dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; dthalf = 0.5 * update->dt; drag_factor = 1.0 - (update->dt * t_freq * drag); @@ -435,7 +436,6 @@ void FixNVT::reset_dt() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - dtq = 0.5 * update->dt; dthalf = 0.5 * update->dt; drag_factor = 1.0 - (update->dt * t_freq * drag); diff --git a/src/fix_nvt.h b/src/fix_nvt.h index 18d6683370..b172813a35 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -23,7 +23,7 @@ class FixNVT : public Fix { FixNVT(class LAMMPS *, int, char **); virtual ~FixNVT(); int setmask(); - void init(); + virtual void init(); void setup(int); virtual void initial_integrate(int); virtual void final_integrate(); @@ -34,7 +34,7 @@ class FixNVT : public Fix { void restart(char *); int modify_param(int, char **); void reset_target(double); - void reset_dt(); + virtual void reset_dt(); protected: int which; @@ -42,7 +42,7 @@ class FixNVT : public Fix { double t_current,t_target; double t_freq,drag,drag_factor; double f_eta,eta_dot,eta,factor; - double dtv,dtf,dtq,dthalf; + double dtv,dtf,dthalf; int nlevels_respa; double *step_respa; diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp index 42bed93bd7..423e6468e9 100644 --- a/src/fix_nvt_sllod.cpp +++ b/src/fix_nvt_sllod.cpp @@ -35,6 +35,7 @@ using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp +enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ @@ -47,6 +48,9 @@ void FixNVTSlodd::init() { FixNVT::init(); + if (!temperature->tempbias) + error->all("Temperature for fix nvt/sllod does not have a bias"); + // check fix deform remap settings int i; @@ -78,11 +82,10 @@ void FixNVTSlodd::initial_integrate(int vflag) eta += dtv*eta_dot; factor = exp(-dthalf*eta_dot); - // 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 - // vdelu = Hrate*Hinv*vthermal + // update v and x of only atoms in group + // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo + // thermostat thermal velocity only + // vdelu = SLLOD correction = Hrate*Hinv*vthermal double **x = atom->x; double **v = atom->v; @@ -93,34 +96,21 @@ void FixNVTSlodd::initial_integrate(int vflag) int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double *h_rate = domain->h_rate; - double *h_ratelo = domain->h_ratelo; - double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; - MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + double h_two[6],vdelu[3]; + MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2]; + vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2]; + vdelu[2] = h_two[2]*v[i][2]; dtfm = dtf / mass[type[i]]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + temperature->restore_bias(v[i]); - domain->x2lamda(x[i],lamda); - vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + - h_rate[4]*lamda[2] + h_ratelo[0]; - vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; - vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; - vthermal[0] = v[i][0] - vstream[0]; - vthermal[1] = v[i][1] - vstream[1]; - vthermal[2] = v[i][2] - vstream[2]; - vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + - h_two[4]*vthermal[2]; - vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; - vdelu[2] = h_two[2]*vthermal[2]; - - v[i][0] = vstream[0] + - vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; - v[i][1] = vstream[1] + - vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; - v[i][2] = vstream[2] + - vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; x[i][0] += dtv * v[i][0]; x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; @@ -134,11 +124,10 @@ void FixNVTSlodd::final_integrate() { double dtfm; - // update vthermal of only atoms in group - // lamda = 0-1 triclinic lamda coords - // vstream = streaming velocity = Hrate*lamda + Hratelo - // vthermal = thermal velocity = v - vstream - // vdelu = Hrate*Hinv*vthermal + // update v of only atoms in group + // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo + // thermostat thermal velocity only + // vdelu = SLLOD correction = Hrate*Hinv*vthermal double **x = atom->x; double **v = atom->v; @@ -149,34 +138,21 @@ void FixNVTSlodd::final_integrate() int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double *h_rate = domain->h_rate; - double *h_ratelo = domain->h_ratelo; - double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; - MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + double h_two[6],vdelu[3]; + MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2]; + vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2]; + vdelu[2] = h_two[2]*v[i][2]; + //dtfm = dtf / mass[type[i]] * factor; dtfm = dtf / mass[type[i]]; - - domain->x2lamda(x[i],lamda); - vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + - h_rate[4]*lamda[2] + h_ratelo[0]; - vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; - vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; - vthermal[0] = v[i][0] - vstream[0]; - vthermal[1] = v[i][1] - vstream[1]; - vthermal[2] = v[i][2] - vstream[2]; - vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + - h_two[4]*vthermal[2]; - vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; - vdelu[2] = h_two[2]*vthermal[2]; - - v[i][0] = vstream[0] + - vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; - v[i][1] = vstream[1] + - vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; - v[i][2] = vstream[2] + - vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + temperature->restore_bias(v[i]); } } @@ -234,34 +210,20 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag) // update v of only atoms in group - double *h_rate = domain->h_rate; - double *h_ratelo = domain->h_ratelo; - double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3]; - MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two); + double h_two[6],vdelu[3]; + MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2]; + vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2]; + vdelu[2] = h_two[2]*v[i][2]; dtfm = dtf / mass[type[i]]; - - domain->x2lamda(x[i],lamda); - vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + - h_rate[4]*lamda[2] + h_ratelo[0]; - vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; - vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; - vthermal[0] = v[i][0] - vstream[0]; - vthermal[1] = v[i][1] - vstream[1]; - vthermal[2] = v[i][2] - vstream[2]; - vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + - h_two[4]*vthermal[2]; - vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2]; - vdelu[2] = h_two[2]*vthermal[2]; - - v[i][0] = vstream[0] + - vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; - v[i][1] = vstream[1] + - vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; - v[i][2] = vstream[2] + - vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + v[i][0] = v[i][0]*factor + dtfm*f[i][0] - dthalf*vdelu[0]; + v[i][1] = v[i][1]*factor + dtfm*f[i][1] - dthalf*vdelu[1]; + v[i][2] = v[i][2]*factor + dtfm*f[i][2] - dthalf*vdelu[2]; + temperature->restore_bias(v[i]); } } diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h index 7adaac67eb..b0a5955c22 100644 --- a/src/fix_nvt_sllod.h +++ b/src/fix_nvt_sllod.h @@ -21,7 +21,6 @@ namespace LAMMPS_NS { class FixNVTSlodd : public FixNVT { public: FixNVTSlodd(class LAMMPS *, int, char **); - ~FixNVTSlodd() {} void init(); void initial_integrate(int); void final_integrate(); diff --git a/src/style.h b/src/style.h index ec8952d41a..a255fdef8c 100644 --- a/src/style.h +++ b/src/style.h @@ -164,12 +164,14 @@ DumpStyle(xyz,DumpXYZ) #include "fix_momentum.h" #include "fix_nph.h" #include "fix_npt.h" +#include "fix_npt_sphere.h" #include "fix_nve.h" #include "fix_nve_limit.h" #include "fix_nve_noforce.h" #include "fix_nve_sphere.h" #include "fix_nvt.h" #include "fix_nvt_sllod.h" +#include "fix_nvt_sphere.h" #include "fix_plane_force.h" #include "fix_press_berendsen.h" #include "fix_print.h" @@ -220,12 +222,14 @@ FixStyle(momentum,FixMomentum) FixStyle(msd,FixMSD) FixStyle(nph,FixNPH) FixStyle(npt,FixNPT) +FixStyle(npt/sphere,FixNPTSphere) FixStyle(nve,FixNVE) FixStyle(nve/limit,FixNVELimit) FixStyle(nve/noforce,FixNVENoforce) FixStyle(nve/sphere,FixNVESphere) FixStyle(nvt,FixNVT) FixStyle(nvt/sllod,FixNVTSlodd) +FixStyle(nvt/sphere,FixNVTSphere) FixStyle(orient/fcc,FixOrientFCC) FixStyle(press/berendsen,FixPressBerendsen) FixStyle(planeforce,FixPlaneForce)