diff --git a/src/fix.h b/src/fix.h index 0698ffcb3c..41fe83a0cc 100644 --- a/src/fix.h +++ b/src/fix.h @@ -125,6 +125,7 @@ class Fix : protected Pointers { virtual void min_clearstore() {} virtual void min_pushstore() {} virtual void min_popstore() {} + virtual int min_reset_ref() {return 0;} virtual void min_step(double, double *) {} virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp index c8530fefa3..f8cdb09d82 100644 --- a/src/fix_box_relax.cpp +++ b/src/fix_box_relax.cpp @@ -27,13 +27,17 @@ #include "modify.h" #include "compute.h" #include "error.h" +#include "math_extra.h" using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) -enum{XYZ,XY,YZ,XZ,ANISO}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO,TRICLINIC}; + +#define MAX_LIFO_DEPTH 2 // 3 box0 arrays in *.h dimensioned to this /* ---------------------------------------------------------------------- */ @@ -42,66 +46,118 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : { if (narg < 4) error->all("Illegal fix box/relax command"); + scalar_flag = 1; + extscalar = 1; + global_freq = 1; box_change = 1; + no_change_box = 1; - if (strcmp(arg[3],"xyz") == 0) { - if (narg < 5) error->all("Illegal fix box/relax command"); - press_couple = XYZ; - p_target[0] = p_target[1] = p_target[2] = atof(arg[4]); - p_flag[0] = p_flag[1] = p_flag[2] = 1; - if (domain->dimension == 2) p_flag[2] = 0; - - } else { - if (strcmp(arg[3],"xy") == 0) press_couple = XY; - else if (strcmp(arg[3],"yz") == 0) press_couple = YZ; - else if (strcmp(arg[3],"xz") == 0) press_couple = XZ; - else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO; - else error->all("Illegal fix box/relax command"); - - if (narg < 7) error->all("Illegal fix box/relax command"); - - if (domain->dimension == 2 && - (press_couple == XY || press_couple == YZ || press_couple == XZ)) - error->all("Invalid fix box/relax command for a 2d simulation"); - - if (strcmp(arg[4],"NULL") == 0) { - p_target[0] = 0.0; - p_flag[0] = 0; - } else { - p_target[0] = atof(arg[4]); - p_flag[0] = 1; - } - if (strcmp(arg[5],"NULL") == 0) { - p_target[1] = 0.0; - p_flag[1] = 0; - } else { - p_target[1] = atof(arg[5]); - p_flag[1] = 1; - } - if (strcmp(arg[6],"NULL") == 0) { - p_target[2] = 0.0; - p_flag[2] = 0; - } else { - if (domain->dimension == 2) - error->all("Invalid fix box/relax command for a 2d simulation"); - p_target[2] = atof(arg[6]); - p_flag[2] = 1; - } - } - - pflagsum = p_flag[0] + p_flag[1] + p_flag[2]; - - // process extra keywords + // default values + pcouple = NONE; allremap = 1; vmax = 0.0001; + deviatoric_flag = 0; + nreset_h0 = 0; - int iarg; - if (press_couple == XYZ) iarg = 5; - else iarg = 7; + p_target[0] = p_target[1] = p_target[2] = + p_target[3] = p_target[4] = p_target[5] = 0.0; + + // process keywords + + dimension = domain->dimension; + + int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"dilate") == 0) { + if (strcmp(arg[iarg],"iso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + pcouple = XYZ; + p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_target[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + pcouple = NONE; + p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_target[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"tri") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + pcouple = NONE; + p_target[0] = p_target[1] = p_target[2] = atof(arg[iarg+1]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + p_target[3] = p_target[4] = p_target[5] = 0.0; + p_flag[3] = p_flag[4] = p_flag[5] = 1; + if (dimension == 2) { + p_target[2] = p_target[3] = p_target[4] = 0.0; + p_flag[2] = p_flag[3] = p_flag[4] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[0] = atof(arg[iarg+1]); + p_flag[0] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[1] = atof(arg[iarg+1]); + p_flag[1] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[2] = atof(arg[iarg+1]); + p_flag[2] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix box/relax command for a 2d simulation"); + + } else if (strcmp(arg[iarg],"yz") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[3] = atof(arg[iarg+1]); + p_flag[3] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix box/relax command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xz") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[4] = atof(arg[iarg+1]); + p_flag[4] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix box/relax command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xy") == 0) { + if (iarg+4 > narg) error->all("Illegal fix box/relax command"); + p_target[5] = atof(arg[iarg+1]); + p_flag[5] = 1; + deviatoric_flag = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all("Illegal fix box/relax command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all("Illegal fix box/relax command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"dilate") == 0) { if (iarg+2 > narg) error->all("Illegal fix box/relax command"); if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0; @@ -111,24 +167,29 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all("Illegal fix box/relax command"); vmax = atof(arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"nreset") == 0) { + if (iarg+2 > narg) error->all("Illegal fix box/relax command"); + nreset_h0 = atoi(arg[iarg+1]); + if (nreset_h0 < 0) error->all("Illegal fix box/relax command"); + iarg += 2; } else error->all("Illegal fix box/relax command"); } // error checks - if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Illegal fix box/relax command"); - if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) - error->all("Illegal fix box/relax command"); - if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) - error->all("Illegal fix box/relax command"); + if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4])) + error->all("Invalid fix box/relax command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all("Invalid fix box/relax command for a 2d simulation"); - if (press_couple == XY && p_target[0] != p_target[1]) - error->all("Illegal fix box/relax command"); - if (press_couple == YZ && p_target[1] != p_target[2]) - error->all("Illegal fix box/relax command"); - if (press_couple == XZ && p_target[0] != p_target[2]) - error->all("Illegal fix box/relax command"); + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0 || p_flag[2] == 0)) + error->all("Invalid fix box/relax command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all("Invalid fix box/relax command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all("Invalid fix box/relax command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all("Invalid fix box/relax command pressure settings"); if (p_flag[0] && domain->xperiodic == 0) error->all("Cannot use fix box/relax on a non-periodic dimension"); @@ -136,9 +197,37 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : error->all("Cannot use fix box/relax on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) error->all("Cannot use fix box/relax on a non-periodic dimension"); + if (p_flag[3] && domain->zperiodic == 0) + error->all("Cannot use fix box/relax on a 2nd non-periodic dimension"); + if (p_flag[4] && domain->zperiodic == 0) + error->all("Cannot use fix box/relax on a 2nd non-periodic dimension"); + if (p_flag[5] && domain->yperiodic == 0) + error->all("Cannot use fix box/relax on a 2nd non-periodic dimension"); + + if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5])) + error->all("Can not specify Pxy/Pxz/Pyz in " + "fix box/relax with non-triclinic box"); + + if (pcouple == XYZ && + (p_target[0] != p_target[1] || p_target[0] != p_target[2])) + error->all("Invalid fix box/relax pressure settings"); + if (pcouple == XY && p_target[0] != p_target[1]) + error->all("Invalid fix box/relax pressure settings"); + if (pcouple == YZ && p_target[1] != p_target[2]) + error->all("Invalid fix box/relax pressure settings"); + if (pcouple == XZ && p_target[0] != p_target[2]) + error->all("Invalid fix box/relax pressure settings"); if (vmax <= 0.0) error->all("Illegal fix box/relax command"); + // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof + // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof + // else pstyle = ANISO -> 3 dof + + if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC; + else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; + else pstyle = ANISO; + // create a new compute temp style // id = fix-ID + temp // compute group = all since pressure is always global (group all) @@ -221,13 +310,6 @@ void FixBoxRelax::init() if (icompute < 0) error->all("Pressure ID for fix box/relax does not exist"); pressure = modify->compute[icompute]; - // initial box dimensions - - xprdinit = domain->xprd; - yprdinit = domain->yprd; - zprdinit = domain->zprd; - if (dimension == 3) volinit = domain->xprd*domain->yprd*domain->zprd; - else volinit = domain->xprd*domain->yprd; pv2e = 1.0 / force->nktv2p; // detect if any rigid fixes exist so rigid bodies move when box is remapped @@ -245,6 +327,19 @@ void FixBoxRelax::init() for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; } + + // initial box dimensions + + xprdinit = domain->xprd; + yprdinit = domain->yprd; + zprdinit = domain->zprd; + if (dimension == 2) zprdinit = 1.0; + vol0 = xprdinit * yprdinit * zprdinit; + + // hydrostatic target pressure and deviatoric target stress + + compute_press_target(); + if (deviatoric_flag) compute_sigma(); } /* ---------------------------------------------------------------------- @@ -253,12 +348,11 @@ void FixBoxRelax::init() double FixBoxRelax::min_energy(double *fextra) { - double eng,scale,scalex,scaley,scalez; + double eng,scale,scalex,scaley,scalez,scalevol; double t_current = temperature->compute_scalar(); - if (press_couple == XYZ) { - double tmp = pressure->compute_scalar(); - } else { + if (pstyle == ISO) double tmp = pressure->compute_scalar(); + else { temperature->compute_vector(); pressure->compute_vector(); } @@ -272,14 +366,14 @@ double FixBoxRelax::min_energy(double *fextra) // returned eng = PV must be in units of energy // returned fextra must likewise be in units of energy - if (press_couple == XYZ) { + if (pstyle == ISO) { scale = domain->xprd/xprdinit; if (dimension == 3) { - eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*volinit; - fextra[0] = pv2e * (p_current[0] - p_target[0])*3.0*scale*scale*volinit; + eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*vol0; + fextra[0] = pv2e * (p_current[0] - p_target[0])*3.0*scale*scale*vol0; } else { - eng = pv2e * p_target[0] * (scale*scale-1.0)*volinit; - fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*volinit; + eng = pv2e * p_target[0] * (scale*scale-1.0)*vol0; + fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*vol0; } } else { @@ -288,15 +382,38 @@ double FixBoxRelax::min_energy(double *fextra) if (p_flag[0]) scalex = domain->xprd/xprdinit; if (p_flag[1]) scaley = domain->yprd/yprdinit; if (p_flag[2]) scalez = domain->zprd/zprdinit; - eng = pv2e * (p_flag[0]*p_target[0] + p_flag[1]*p_target[1] + - p_flag[2]*p_target[2])/pflagsum * - (scalex*scaley*scalez-1.0)*volinit; + scalevol = scalex*scaley*scalez; + eng = pv2e * p_hydro * (scalevol-1.0)*vol0; if (p_flag[0]) - fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*scalez*volinit; + fextra[0] = pv2e * (p_current[0] - p_hydro)*scaley*scalez*vol0; if (p_flag[1]) - fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*scalez*volinit; + fextra[1] = pv2e * (p_current[1] - p_hydro)*scalex*scalez*vol0; if (p_flag[2]) - fextra[2] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit; + fextra[2] = pv2e * (p_current[2] - p_hydro)*scalex*scaley*vol0; + + if (pstyle == TRICLINIC) { + fextra[3] = fextra[4] = fextra[5] = 0.0; + if (p_flag[3]) + fextra[3] = pv2e*p_current[3]*scaley*yprdinit*scalex*xprdinit*yprdinit; + if (p_flag[4]) + fextra[4] = pv2e*p_current[4]*scalex*xprdinit*scaley*yprdinit*xprdinit; + if (p_flag[5]) + fextra[5] = pv2e*p_current[5]*scalex*xprdinit*scalez*zprdinit*xprdinit; + } + + if (deviatoric_flag) { + compute_deviatoric(); + if (p_flag[0]) fextra[0] -= fdev[0]*xprdinit; + if (p_flag[1]) fextra[1] -= fdev[1]*yprdinit; + if (p_flag[2]) fextra[2] -= fdev[2]*zprdinit; + if (pstyle == TRICLINIC) { + if (p_flag[3]) fextra[3] -= fdev[3]*yprdinit; + if (p_flag[4]) fextra[4] -= fdev[4]*xprdinit; + if (p_flag[5]) fextra[5] -= fdev[5]*xprdinit; + } + + eng += compute_strain_energy(); + } } return eng; @@ -319,6 +436,11 @@ void FixBoxRelax::min_store() s0[0] = (boxhi0[current_lifo][0]-boxlo0[current_lifo][0])/xprdinit; s0[1] = (boxhi0[current_lifo][1]-boxlo0[current_lifo][1])/yprdinit; s0[2] = (boxhi0[current_lifo][2]-boxlo0[current_lifo][2])/zprdinit; + if (pstyle == TRICLINIC) { + boxtilt0[current_lifo][0] = domain->yz; + boxtilt0[current_lifo][1] = domain->xz; + boxtilt0[current_lifo][2] = domain->xy; + } } /* ---------------------------------------------------------------------- @@ -340,7 +462,6 @@ void FixBoxRelax::min_pushstore() error->all("Attempt to push beyond stack limit in fix box/relax"); return; } - current_lifo++; } @@ -355,24 +476,51 @@ void FixBoxRelax::min_popstore() error->all("Attempt to pop empty stack in fix box/relax"); return; } - current_lifo--; } +/* ---------------------------------------------------------------------- + check if time to reset reference state. If so, do so. +------------------------------------------------------------------------- */ + +int FixBoxRelax::min_reset_ref() +{ + int itmp = 0; + + // if nreset_h0 > 0, reset reference box + // every nreset_h0 timesteps + // only needed for deviatoric external stress + + if (deviatoric_flag && nreset_h0 > 0) { + int delta = update->ntimestep - update->beginstep; + if (delta % nreset_h0 == 0) { + compute_sigma(); + itmp = 1; + } + } + return itmp; +} + /* ---------------------------------------------------------------------- change the box dimensions by fraction ds = alpha*hextra ------------------------------------------------------------------------- */ void FixBoxRelax::min_step(double alpha, double *hextra) { - if (press_couple == XYZ) { + if (pstyle == ISO) { ds[0] = ds[1] = ds[2] = alpha*hextra[0]; } else { + ds[0] = ds[1] = ds[2] = 0.0; if (p_flag[0]) ds[0] = alpha*hextra[0]; if (p_flag[1]) ds[1] = alpha*hextra[1]; if (p_flag[2]) ds[2] = alpha*hextra[2]; + if (pstyle == TRICLINIC) { + ds[3] = ds[4] = ds[5] = 0.0; + if (p_flag[3]) ds[3] = alpha*hextra[3]; + if (p_flag[4]) ds[4] = alpha*hextra[4]; + if (p_flag[5]) ds[5] = alpha*hextra[5]; + } } - remap(); } @@ -382,12 +530,17 @@ void FixBoxRelax::min_step(double alpha, double *hextra) double FixBoxRelax::max_alpha(double *hextra) { - double alpha = 0.0; - if (press_couple == XYZ) alpha = vmax/fabs(hextra[0]); + double alpha = 1.0; + if (pstyle == ISO) alpha = vmax/fabs(hextra[0]); else { - alpha = vmax/fabs(hextra[0]); - alpha = MIN(alpha,vmax/fabs(hextra[1])); - alpha = MIN(alpha,vmax/fabs(hextra[2])); + if (p_flag[0]) alpha = MIN(alpha,vmax/fabs(hextra[0])); + if (p_flag[1]) alpha = MIN(alpha,vmax/fabs(hextra[1])); + if (p_flag[2]) alpha = MIN(alpha,vmax/fabs(hextra[2])); + if (pstyle == TRICLINIC) { + if (p_flag[3]) alpha = MIN(alpha,vmax/fabs(hextra[3])); + if (p_flag[4]) alpha = MIN(alpha,vmax/fabs(hextra[4])); + if (p_flag[5]) alpha = MIN(alpha,vmax/fabs(hextra[5])); + } } return alpha; } @@ -398,7 +551,8 @@ double FixBoxRelax::max_alpha(double *hextra) int FixBoxRelax::min_dof() { - if (press_couple == XYZ) return 1; + if (pstyle == ISO) return 1; + if (pstyle == TRICLINIC) return 6; return 3; } @@ -441,8 +595,16 @@ void FixBoxRelax::remap() ctr = 0.5 * (currentBoxLo0 + currentBoxHi0); domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0-ctr)*ds[i]/s0[i]; domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0-ctr)*ds[i]/s0[i]; + if (domain->boxlo[i] >= domain->boxhi[i]) + error->all("fix box/relax generated negative dimension"); } + if (pstyle == TRICLINIC) { + if (p_flag[3]) domain->yz = boxtilt0[current_lifo][0]+ds[3]*yprdinit; + if (p_flag[4]) domain->xz = boxtilt0[current_lifo][1]+ds[4]*xprdinit; + if (p_flag[5]) domain->xy = boxtilt0[current_lifo][2]+ds[5]*xprdinit; + } + domain->set_global_box(); domain->set_local_box(); @@ -466,25 +628,36 @@ void FixBoxRelax::couple() { double *tensor = pressure->vector; - if (press_couple == XYZ) + if (pstyle == ISO) p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == XY) { + else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; - } else if (press_couple == YZ) { + } else if (pcouple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; - } else if (press_couple == XZ) { + } else if (pcouple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; - } else if (press_couple == ANISO) { + } else { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; } + + // switch order from xy-xz-yz to Voigt + + if (pstyle == TRICLINIC) { + p_current[3] = tensor[5]; + p_current[4] = tensor[4]; + p_current[5] = tensor[3]; + } } /* ---------------------------------------------------------------------- */ @@ -540,3 +713,187 @@ int FixBoxRelax::modify_param(int narg, char **arg) } return 0; } + +/* ---------------------------------------------------------------------- + compute sigma tensor (needed whenever reference box is reset) +-----------------------------------------------------------------------*/ + +void FixBoxRelax::compute_sigma() +{ + double pdevmod[3][3],pdeviatoric[3][3],htmp[3][3]; + double tmp1[3][3],sigma_tensor[3][3],h_invtmp[3][3]; + + // reset reference box dimensions + + xprdinit = domain->xprd; + yprdinit = domain->yprd; + zprdinit = domain->zprd; + if (dimension == 2) zprdinit = 1.0; + vol0 = xprdinit * yprdinit * zprdinit; + + h0[0] = domain->h[0]; + h0[1] = domain->h[1]; + h0[2] = domain->h[2]; + h0[3] = domain->h[3]; + h0[4] = domain->h[4]; + h0[5] = domain->h[5]; + + h0_inv[0] = domain->h_inv[0]; + h0_inv[1] = domain->h_inv[1]; + h0_inv[2] = domain->h_inv[2]; + h0_inv[3] = domain->h_inv[3]; + h0_inv[4] = domain->h_inv[4]; + h0_inv[5] = domain->h_inv[5]; + + htmp[0][0] = h0[0]; + htmp[1][1] = h0[1]; + htmp[2][2] = h0[2]; + htmp[1][2] = h0[3]; + htmp[0][2] = h0[4]; + htmp[0][1] = h0[5]; + htmp[2][1] = 0.0; + htmp[2][0] = 0.0; + htmp[1][0] = 0.0; + + h_invtmp[0][0] = h0_inv[0]; + h_invtmp[1][1] = h0_inv[1]; + h_invtmp[2][2] = h0_inv[2]; + h_invtmp[1][2] = h0_inv[3]; + h_invtmp[0][2] = h0_inv[4]; + h_invtmp[0][1] = h0_inv[5]; + h_invtmp[2][1] = 0.0; + h_invtmp[2][0] = 0.0; + h_invtmp[1][0] = 0.0; + + // compute target deviatoric stress tensor pdevmod + + pdeviatoric[0][0] = pdeviatoric[1][1] = pdeviatoric[2][2] = 0.0; + if (p_flag[0]) pdeviatoric[0][0] = p_target[0] - p_hydro; + if (p_flag[1]) pdeviatoric[1][1] = p_target[1] - p_hydro; + if (p_flag[2]) pdeviatoric[2][2] = p_target[2] - p_hydro; + pdeviatoric[1][2] = pdeviatoric[2][1] = p_target[3]; + pdeviatoric[0][2] = pdeviatoric[2][0] = p_target[4]; + pdeviatoric[0][1] = pdeviatoric[1][0] = p_target[5]; + + // Modify to account for off-diagonal terms + // These equations come from the stationarity relation: + // Pdev,sys = Pdev,targ*hinv^t*hdiag + // where: + // Pdev,sys is the system deviatoric stress tensor, + // Pdev,targ = pdeviatoric, effective target deviatoric stress + // hinv^t is the transpose of the inverse h tensor + // hdiag is the diagonal part of the h tensor + + pdeviatoric[1][1] -= pdeviatoric[1][2]*h0_inv[3]*h0[1]; + pdeviatoric[0][1] -= pdeviatoric[0][2]*h0_inv[3]*h0[1]; + pdeviatoric[1][0] = pdeviatoric[0][1]; + pdeviatoric[0][0] -= pdeviatoric[0][1]*h0_inv[5]*h0[0] + + pdeviatoric[0][2]*h0_inv[4]*h0[0]; + + // compute symmetric sigma tensor + + MathExtra::times3(h_invtmp,pdeviatoric,tmp1); + MathExtra::times3_transpose(tmp1,h_invtmp,sigma_tensor); + MathExtra::scalar_times3(vol0,sigma_tensor); + + sigma[0] = sigma_tensor[0][0]; + sigma[1] = sigma_tensor[1][1]; + sigma[2] = sigma_tensor[2][2]; + sigma[3] = sigma_tensor[1][2]; + sigma[4] = sigma_tensor[0][2]; + sigma[5] = sigma_tensor[0][1]; +} + +/* ---------------------------------------------------------------------- + compute strain energy +-----------------------------------------------------------------------*/ + +double FixBoxRelax::compute_strain_energy() +{ + // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units + + double* h = domain->h; + double d0,d1,d2; + + if (dimension == 3) { + d0 = + sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) + + sigma[5]*( h[1]*h[5]+h[3]*h[4]) + + sigma[4]*( h[2]*h[4]); + d1 = + sigma[5]*( h[5]*h[1]+h[4]*h[3]) + + sigma[1]*( h[1]*h[1]+h[3]*h[3]) + + sigma[3]*( h[2]*h[3]); + d2 = + sigma[4]*( h[4]*h[2]) + + sigma[3]*( h[3]*h[2]) + + sigma[2]*( h[2]*h[2]); + } else { + d0 = sigma[0]*(h[0]*h[0]+h[5]*h[5]) + sigma[5]*h[1]*h[5]; + d1 = sigma[5]*h[5]*h[1] + sigma[1]*h[1]*h[1]; + d2 = 0.0; + } + + double energy = 0.5*(d0+d1+d2)*pv2e; + return energy; +} + +/* ---------------------------------------------------------------------- + compute deviatoric barostat force = h*sigma*h^t +-----------------------------------------------------------------------*/ + +void FixBoxRelax::compute_deviatoric() +{ + double* h = domain->h; + + // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] + // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] + // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] + + if (dimension == 3) { + fdev[0] = pv2e*(h[0]*sigma[0]+h[5]*sigma[5]+h[4]*sigma[4]); + fdev[1] = pv2e*(h[1]*sigma[1]+h[3]*sigma[3]); + fdev[2] = pv2e*(h[2]*sigma[2]); + fdev[3] = pv2e*(h[1]*sigma[3]+h[3]*sigma[2]); + fdev[4] = pv2e*(h[0]*sigma[4]+h[5]*sigma[3]+h[4]*sigma[2]); + fdev[5] = pv2e*(h[0]*sigma[5]+h[5]*sigma[1]+h[4]*sigma[3]); + } else { + fdev[0] = pv2e*(h[0]*sigma[0]+h[5]*sigma[5]); + fdev[1] = pv2e*(h[1]*sigma[1]); + fdev[5] = pv2e*(h[0]*sigma[5]+h[5]*sigma[1]); + } +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixBoxRelax::compute_press_target() +{ + pflagsum = p_flag[0] + p_flag[1] + p_flag[2]; + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) p_hydro += p_target[i]; + if (pflagsum) p_hydro /= pflagsum; + + for (int i = 0; i < 3; i++) { + if (p_flag[i] && fabs(p_hydro - p_target[i] > 1.0e-6)) deviatoric_flag = 1; + } + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i] && fabs(p_target[i]) > 1.0e-6) deviatoric_flag = 1; + } +} + +/* ---------------------------------------------------------------------- + compute PV and strain energy for access to the user +/* ---------------------------------------------------------------------- */ + +double FixBoxRelax::compute_scalar() +{ + double ftmp[6] = {0.0,0.0,0.0,0.0,0.0,0.0}; + if (update->ntimestep == 0) return 0.0; + return min_energy(ftmp); +} diff --git a/src/fix_box_relax.h b/src/fix_box_relax.h index 5ca43c7f1b..33230131ff 100644 --- a/src/fix_box_relax.h +++ b/src/fix_box_relax.h @@ -22,8 +22,6 @@ FixStyle(box/relax,FixBoxRelax) #include "fix.h" -const int MAX_LIFO_DEPTH = 2; // to meet the needs of min_hftn - namespace LAMMPS_NS { class FixBoxRelax : public Fix { @@ -38,6 +36,7 @@ class FixBoxRelax : public Fix { void min_clearstore(); void min_pushstore(); void min_popstore(); + int min_reset_ref(); void min_step(double, double *); double max_alpha(double *); int min_dof(); @@ -45,19 +44,19 @@ class FixBoxRelax : public Fix { int modify_param(int, char **); private: - int p_flag[3]; - int press_couple,allremap; + int p_flag[6]; + int pstyle,pcouple,allremap; int dimension; - double p_target[3],p_current[3]; - double dilation[3]; - double volinit,xprdinit,yprdinit,zprdinit; + double p_target[6],p_current[6]; + double vol0,xprdinit,yprdinit,zprdinit; double vmax,pv2e,pflagsum; - int current_lifo; // LIFO stack pointer - double boxlo0[MAX_LIFO_DEPTH][3]; // box bounds at start of line search - double boxhi0[MAX_LIFO_DEPTH][3]; - double s0[6]; // scale matrix in Voigt notation - double ds[6]; // increment in scale matrix + int current_lifo; // LIFO stack pointer + double boxlo0[2][3]; // box bounds at start of line search + double boxhi0[2][3]; + double boxtilt0[2][3]; // xy,xz,yz tilts at start of line search + double s0[3]; // scale matrix at start of line search + double ds[6]; // increment in scale matrix char *id_temp,*id_press; class Compute *temperature,*pressure; @@ -66,8 +65,25 @@ class FixBoxRelax : public Fix { int nrigid; int *rfix; + double sigma[6]; // scaled target stress + double utsigma[3]; // weighting for upper-tri elements + // of modified sigma + int sigmamod_flag; // 1 if modified sigma to be used + double fdev[6]; // Deviatoric force on cell + int deviatoric_flag; // 0 if target stress tensor is hydrostatic + double h0[6]; // h_inv of reference (zero strain) box + double h0_inv[6]; // h_inv of reference (zero strain) box + int nreset_h0; // interval for resetting h0 + double p_hydro; // hydrostatic component of target stress + void remap(); void couple(); + + void compute_sigma(); + void compute_deviatoric(); + double compute_strain_energy(); + void compute_press_target(); + double compute_scalar(); }; } diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 10aeb26130..09a6a001bc 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -181,14 +181,14 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if ((set[0].style && domain->xperiodic == 0) || (set[1].style && domain->yperiodic == 0) || (set[2].style && domain->zperiodic == 0)) - error->all("Cannot fix deform on a non-periodic boundary"); + error->all("Cannot use fix deform on a non-periodic boundary"); - if (set[3].style && (domain->yperiodic == 0 || domain->zperiodic == 0)) - error->all("Cannot fix deform on a non-periodic boundary"); - if (set[4].style && (domain->xperiodic == 0 || domain->zperiodic == 0)) - error->all("Cannot fix deform on a non-periodic boundary"); - if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0)) - error->all("Cannot fix deform on a non-periodic boundary"); + if (set[3].style && domain->zperiodic == 0) + error->all("Cannot use fix deform on a 2nd non-periodic boundary"); + if (set[4].style && domain->zperiodic == 0) + error->all("Cannot use fix deform on a 2nd non-periodic boundary"); + if (set[5].style && domain->yperiodic == 0) + error->all("Cannot use fix deform on a 2nd non-periodic boundary"); // apply scaling to FINAL,DELTA,VEL,WIGGLE since they have distance/vel units diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp new file mode 100644 index 0000000000..cda4187b08 --- /dev/null +++ b/src/fix_nh.cpp @@ -0,0 +1,1872 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Mark Stevens (SNL), Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "fix_nh.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "modify.h" +#include "fix_deform.h" +#include "compute.h" +#include "kspace.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "memory.h" +#include "error.h" +#include "math_extra.h" + +using namespace LAMMPS_NS; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +enum{NOBIAS,BIAS}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO,TRICLINIC}; + +/* ---------------------------------------------------------------------- + NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion + ---------------------------------------------------------------------- */ + +FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg < 4) error->all("Illegal fix nvt/npt/nph command"); + + restart_global = 1; + time_integrate = 1; + scalar_flag = 1; + vector_flag = 1; + global_freq = 1; + extscalar = 1; + extvector = 0; + + // default values + + pcouple = NONE; + drag = 0.0; + allremap = 1; + mtchain = mpchain = 3; + nc_tchain = nc_pchain = 1; + mtk_flag = 1; + deviatoric_flag = 0; + nreset_h0 = 0; + + tstat_flag = 0; + double t_period = 0.0; + + double p_period[6]; + for (int i = 0; i < 6; i++) { + p_start[i] = p_stop[i] = p_period[i] = 0.0; + p_flag[i] = 0; + p_period[i] = 0.0; + } + + // process keywords + + dimension = domain->dimension; + + int iarg = 3; + + while (iarg < narg) { + if (strcmp(arg[iarg],"temp") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + tstat_flag = 1; + t_start = atof(arg[iarg+1]); + t_stop = atof(arg[iarg+2]); + t_period = atof(arg[iarg+3]); + if (t_start < 0.0 || t_stop <= 0.0) + error->all("Target T for fix nvt/npt/nph cannot be 0.0"); + iarg += 4; + + } else if (strcmp(arg[iarg],"iso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + pcouple = XYZ; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + pcouple = NONE; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"tri") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + pcouple = NONE; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + p_start[3] = p_start[4] = p_start[5] = 0.0; + p_stop[3] = p_stop[4] = p_stop[5] = 0.0; + p_period[3] = p_period[4] = p_period[5] = atof(arg[iarg+3]); + p_flag[3] = p_flag[4] = p_flag[5] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + p_start[3] = p_stop[3] = p_period[3] = 0.0; + p_flag[3] = 0; + p_start[4] = p_stop[4] = p_period[4] = 0.0; + p_flag[4] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[0] = atof(arg[iarg+1]); + p_stop[0] = atof(arg[iarg+2]); + p_period[0] = atof(arg[iarg+3]); + p_flag[0] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[1] = atof(arg[iarg+1]); + p_stop[1] = atof(arg[iarg+2]); + p_period[1] = atof(arg[iarg+3]); + p_flag[1] = 1; + deviatoric_flag = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[2] = atof(arg[iarg+1]); + p_stop[2] = atof(arg[iarg+2]); + p_period[2] = atof(arg[iarg+3]); + p_flag[2] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix nvt/npt/nph command for a 2d simulation"); + + } else if (strcmp(arg[iarg],"yz") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[3] = atof(arg[iarg+1]); + p_stop[3] = atof(arg[iarg+2]); + p_period[3] = atof(arg[iarg+3]); + p_flag[3] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix nvt/npt/nph command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xz") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[4] = atof(arg[iarg+1]); + p_stop[4] = atof(arg[iarg+2]); + p_period[4] = atof(arg[iarg+3]); + p_flag[4] = 1; + deviatoric_flag = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix nvt/npt/nph command for a 2d simulation"); + } else if (strcmp(arg[iarg],"xy") == 0) { + if (iarg+4 > narg) error->all("Illegal fix nvt/npt/nph command"); + p_start[5] = atof(arg[iarg+1]); + p_stop[5] = atof(arg[iarg+2]); + p_period[5] = atof(arg[iarg+3]); + p_flag[5] = 1; + deviatoric_flag = 1; + iarg += 4; + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"drag") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + drag = atof(arg[iarg+1]); + if (drag < 0.0) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"dilate") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; + else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0; + else error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"tchain") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + mtchain = atoi(arg[iarg+1]); + if (mtchain < 1) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"pchain") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + mpchain = atoi(arg[iarg+1]); + if (mpchain < 0) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"mtk") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0; + else error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"tloop") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + nc_tchain = atoi(arg[iarg+1]); + if (nc_tchain < 0) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"ploop") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + nc_pchain = atoi(arg[iarg+1]); + if (nc_pchain < 0) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else if (strcmp(arg[iarg],"nreset") == 0) { + if (iarg+2 > narg) error->all("Illegal fix nvt/npt/nph command"); + nreset_h0 = atoi(arg[iarg+1]); + if (nreset_h0 < 0) error->all("Illegal fix nvt/npt/nph command"); + iarg += 2; + } else error->all("Illegal fix nvt/npt/nph command"); + } + + // error checks + + if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4])) + error->all("Invalid fix nvt/npt/nph command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all("Invalid fix nvt/npt/nph command for a 2d simulation"); + + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0 || p_flag[2] == 0)) + error->all("Invalid fix nvt/npt/nph command pressure settings"); + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) + error->all("Invalid fix nvt/npt/nph command pressure settings"); + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all("Invalid fix nvt/npt/nph command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) + error->all("Invalid fix nvt/npt/nph command pressure settings"); + + if (p_flag[0] && domain->xperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a non-periodic dimension"); + if (p_flag[1] && domain->yperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a non-periodic dimension"); + if (p_flag[2] && domain->zperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a non-periodic dimension"); + if (p_flag[3] && domain->zperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); + if (p_flag[4] && domain->zperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); + if (p_flag[5] && domain->yperiodic == 0) + error->all("Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); + + if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5])) + error->all("Can not specify Pxy/Pxz/Pyz in " + "fix nvt/npt/nph with non-triclinic box"); + + if (pcouple == XYZ && + (p_start[0] != p_start[1] || p_start[0] != p_start[2] || + p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[1] || p_period[0] != p_period[2])) + error->all("Invalid fix nvt/npt/nph pressure settings"); + if (pcouple == XY && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all("Invalid fix nvt/npt/nph pressure settings"); + if (pcouple == YZ && + (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || + p_period[1] != p_period[2])) + error->all("Invalid fix nvt/npt/nph pressure settings"); + if (pcouple == XZ && + (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[2])) + error->all("Invalid fix nvt/npt/nph pressure settings"); + + if ((tstat_flag && t_period <= 0.0) || + (p_flag[0] && p_period[0] <= 0.0) || + (p_flag[1] && p_period[1] <= 0.0) || + (p_flag[2] && p_period[2] <= 0.0) || + (p_flag[3] && p_period[3] <= 0.0) || + (p_flag[4] && p_period[4] <= 0.0) || + (p_flag[5] && p_period[5] <= 0.0)) + error->all("Fix nvt/npt/nph damping parameters must be > 0.0"); + + // set pstat_flag and box change variables + + pstat_flag = 0; + for (int i = 0; i < 6; i++) + if (p_flag[i]) pstat_flag = 1; + + if (pstat_flag) { + box_change = 1; + no_change_box = 1; + } + + // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof + // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof + // else pstyle = ANISO -> 3 dof + + if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC; + else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; + else pstyle = ANISO; + + // convert input periods to frequencies + + t_freq = 0.0; + p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0; + + if (tstat_flag) t_freq = 1.0 / t_period; + if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; + if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; + if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + if (p_flag[3]) p_freq[3] = 1.0 / p_period[3]; + if (p_flag[4]) p_freq[4] = 1.0 / p_period[4]; + if (p_flag[5]) p_freq[5] = 1.0 / p_period[5]; + + // Nose/Hoover temp and pressure init + + size_vector = 0; + + if (tstat_flag) { + int ich; + eta = new double[mtchain]; + + // add one extra dummy thermostat, set to zero + + eta_dot = new double[mtchain+1]; + eta_dot[mtchain] = 0.0; + eta_dotdot = new double[mtchain]; + for (ich = 0; ich < mtchain; ich++) { + eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0; + } + eta_mass = new double[mtchain]; + size_vector += 2*2*mtchain; + } + + if (pstat_flag) { + omega[0] = omega[1] = omega[2] = 0.0; + omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0; + omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0; + omega[3] = omega[4] = omega[5] = 0.0; + omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0; + omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0; + if (pstyle == ISO) size_vector += 2*2*1; + else if (pstyle == ANISO) size_vector += 2*2*3; + else if (pstyle == TRICLINIC) size_vector += 2*2*6; + + if (mpchain) { + int ich; + etap = new double[mpchain]; + + // add one extra dummy thermostat, set to zero + + etap_dot = new double[mpchain+1]; + etap_dot[mpchain] = 0.0; + etap_dotdot = new double[mpchain]; + for (ich = 0; ich < mpchain; ich++) { + etap[ich] = etap_dot[ich] = + etap_dotdot[ich] = 0.0; + } + etap_mass = new double[mpchain]; + size_vector += 2*2*mpchain; + } + + if (deviatoric_flag) size_vector += 1; + } + + nrigid = 0; + rfix = NULL; + + // initialize vol0,t0 to zero to signal uninitialized + // values then assigned in init(), if necessary + + vol0 = t0 = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +FixNH::~FixNH() +{ + delete [] rfix; + + // delete temperature and pressure if fix created them + + if (tflag) modify->delete_compute(id_temp); + delete [] id_temp; + + if (tstat_flag) { + delete [] eta; + delete [] eta_dot; + delete [] eta_dotdot; + delete [] eta_mass; + } + + if (pstat_flag) { + if (pflag) modify->delete_compute(id_press); + delete [] id_press; + if (mpchain) { + delete [] etap; + delete [] etap_dot; + delete [] etap_dotdot; + delete [] etap_mass; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNH::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + mask |= THERMO_ENERGY; + mask |= INITIAL_INTEGRATE_RESPA; + mask |= FINAL_INTEGRATE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::init() +{ + // insure no conflict with fix deform + + if (pstat_flag) + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"deform") == 0) { + int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; + if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || + (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) || + (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5])) + error->all("Cannot use fix npt and fix deform on " + "same component of stress tensor"); + } + + // set temperature and pressure ptrs + + int icompute = modify->find_compute(id_temp); + if (icompute < 0) + error->all("Temperature ID for fix nvt/nph/npt does not exist"); + temperature = modify->compute[icompute]; + + if (temperature->tempbias) which = BIAS; + else which = NOBIAS; + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) error->all("Pressure ID for fix npt/nph does not exist"); + pressure = modify->compute[icompute]; + } + + // set timesteps and frequencies + + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dthalf = 0.5 * update->dt; + dt4 = 0.25 * update->dt; + dt8 = 0.125 * update->dt; + + p_freq_max = 0.0; + if (pstat_flag) { + p_freq_max = MAX(p_freq[0],p_freq[1]); + p_freq_max = MAX(p_freq_max,p_freq[2]); + if (pstyle == TRICLINIC) { + p_freq_max = MAX(p_freq_max,p_freq[3]); + p_freq_max = MAX(p_freq_max,p_freq[4]); + p_freq_max = MAX(p_freq_max,p_freq[5]); + } + drag_factor = 1.0 - (update->dt * p_freq_max * drag); + } + + if (tstat_flag) + drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag); + + // tally the number of dimensions that are barostatted + // also compute the initial volume and reference cell + // set initial volume and reference cell, if not already done + + if (pstat_flag) { + pdim = p_flag[0] + p_flag[1] + p_flag[2]; + if (vol0 == 0.0) { + if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; + else vol0 = domain->xprd * domain->yprd; + h0_inv[0] = domain->h_inv[0]; + h0_inv[1] = domain->h_inv[1]; + h0_inv[2] = domain->h_inv[2]; + h0_inv[3] = domain->h_inv[3]; + h0_inv[4] = domain->h_inv[4]; + h0_inv[5] = domain->h_inv[5]; + } + } + + boltz = force->boltz; + nktv2p = force->nktv2p; + + if (force->kspace) kspace_flag = 1; + else kspace_flag = 0; + + if (strcmp(update->integrate_style,"respa") == 0) { + nlevels_respa = ((Respa *) update->integrate)->nlevels; + step_respa = ((Respa *) update->integrate)->step; + } + + // detect if any rigid fixes exist so rigid bodies move when box is remapped + // rfix[] = indices to each fix rigid + + delete [] rfix; + nrigid = 0; + rfix = NULL; + + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigid++; + if (nrigid) { + rfix = new int[nrigid]; + nrigid = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; + } +} + +/* ---------------------------------------------------------------------- + compute T,P before integrator starts +------------------------------------------------------------------------- */ + +void FixNH::setup(int vflag) +{ + // initialize some quantities that were not available earlier + + if (mtk_flag) mtk_factor = 1.0 + 1.0/atom->natoms; + else mtk_factor = 1.0; + tdof = temperature->dof; + + // t_target is used by compute_scalar(), even for NPH + + if (tstat_flag) t_target = t_start; + else if (pstat_flag) { + + // t0 = initial value for piston mass and energy conservation + // cannot be done in init() b/c temperature cannot be called there + // is b/c Modify::init() inits computes after fixes due to dof dependence + // guesstimate a unit-dependent t0 if actual T = 0.0 + // if it was read in from a restart file, leave it be + + if (t0 == 0.0) { + t0 = temperature->compute_scalar(); + if (t0 == 0.0) { + if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0; + else t0 = 300.0; + } + } + t_target = t0; + } + + if (pstat_flag) compute_press_target(); + + t_current = temperature->compute_scalar(); + if (pstat_flag) { + if (pstyle == ISO) double tmp = pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + } + + // initial forces on thermostat variables + + if (tstat_flag) { + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) { + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - + boltz*t_target) / eta_mass[ich]; + } + } + + if (pstat_flag) { + double kt = boltz * t_target; + double nkt = atom->natoms * kt; + + for (int i = 0; i < 3; i++) + if (p_flag[i]) + omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) + if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]); + } + + // initial forces on barostat thermostat variables + + if (mpchain) { + etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max); + for (int ich = 1; ich < mpchain; ich++) + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - + boltz*t_target) / etap_mass[ich]; + } + + // compute appropriately coupled elements of mvv_current + + if (mtk_flag) couple_ke(); + } +} + +/* ---------------------------------------------------------------------- + 1st half of Verlet update +------------------------------------------------------------------------- */ + +void FixNH::initial_integrate(int vflag) +{ + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); + + // update eta_dot + + if (tstat_flag) { + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + t_target = t_start + delta * (t_stop-t_start); + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + nhc_temp_integrate(); + } + + // need to recompute pressure to account for change in KE + // t_current is up-to-date, but compute_temperature is not + // compute appropriately coupled elements of mvv_current + + if (pstat_flag) { + if (pstyle == ISO) { + temperature->compute_scalar(); + double tmp = pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + if (mtk_flag) couple_ke(); + } + + if (pstat_flag) { + compute_press_target(); + nh_omega_dot(); + nh_v_press(); + } + + nve_v(); + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + nve_x(); + + // remap simulation box by 1/2 step + // redo KSpace coeffs since volume has changed + + if (pstat_flag) { + remap(); + if (kspace_flag) force->kspace->setup(); + } +} + +/* ---------------------------------------------------------------------- + 2nd half of Verlet update +------------------------------------------------------------------------- */ + +void FixNH::final_integrate() +{ + nve_v(); + + if (pstat_flag) nh_v_press(); + + // compute new T,P + // compute appropriately coupled elements of mvv_current + + t_current = temperature->compute_scalar(); + if (pstat_flag) { + if (pstyle == ISO) double tmp = pressure->compute_scalar(); + else pressure->compute_vector(); + couple(); + pressure->addstep(update->ntimestep+1); + if (mtk_flag) couple_ke(); + } + + if (pstat_flag) nh_omega_dot(); + + // update eta_dot + + if (tstat_flag) nhc_temp_integrate(); + + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop) +{ + int i; + + // set timesteps by level + + dtv = step_respa[ilevel]; + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dthalf = 0.5 * step_respa[ilevel]; + + // 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 group + + if (ilevel == nlevels_respa-1) { + + // update eta_press_dot + + if (pstat_flag && mpchain) nhc_press_integrate(); + + // update eta_dot + + if (tstat_flag) { + double delta = update->ntimestep - update->beginstep; + delta /= update->endstep - update->beginstep; + t_target = t_start + delta * (t_stop-t_start); + eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq); + for (int ich = 1; ich < mtchain; ich++) + eta_mass[ich] = boltz * t_target / (t_freq*t_freq); + nhc_temp_integrate(); + } + + // recompute pressure to account for change in KE + // t_current is up-to-date, but compute_temperature is not + // compute appropriately coupled elements of mvv_current + + if (pstat_flag) { + if (pstyle == ISO) { + temperature->compute_scalar(); + double tmp = pressure->compute_scalar(); + } else { + temperature->compute_vector(); + pressure->compute_vector(); + } + couple(); + pressure->addstep(update->ntimestep+1); + if (mtk_flag) couple_ke(); + } + + if (pstat_flag) { + compute_press_target(); + nh_omega_dot(); + nh_v_press(); + } + + nve_v(); + + // remap simulation box by 1/2 step + + if (pstat_flag) remap(); + + } else nve_v(); + + // innermost level - also update x only for atoms in group + // if iloop = 0, then is 1st call at innermost level from rRESPA + // perform 2nd half of box remap + // redo KSpace coeffs since volume has changed + + if (ilevel == 0) { + nve_x(); + if (iloop == 0) { + remap(); + if (kspace_flag) force->kspace->setup(); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::final_integrate_respa(int ilevel, int iloop) +{ + // set timesteps by level + + dtf = 0.5 * step_respa[ilevel] * force->ftm2v; + dthalf = 0.5 * step_respa[ilevel]; + + // outermost level - update eta_dot and omega_dot, apply via final_integrate + // all other levels - NVE update of v + + if (ilevel == nlevels_respa-1) final_integrate(); + else nve_v(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::couple() +{ + double *tensor = pressure->vector; + + if (pstyle == ISO) + p_current[0] = p_current[1] = p_current[2] = pressure->scalar; + else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { + double ave = 0.5 * (tensor[0] + tensor[1]); + p_current[0] = p_current[1] = ave; + p_current[2] = tensor[2]; + } else if (pcouple == YZ) { + double ave = 0.5 * (tensor[1] + tensor[2]); + p_current[1] = p_current[2] = ave; + p_current[0] = tensor[0]; + } else if (pcouple == XZ) { + double ave = 0.5 * (tensor[0] + tensor[2]); + p_current[0] = p_current[2] = ave; + p_current[1] = tensor[1]; + } else { + p_current[0] = tensor[0]; + p_current[1] = tensor[1]; + p_current[2] = tensor[2]; + } + + // switch order from xy-xz-yz to Voigt + + if (pstyle == TRICLINIC) { + p_current[3] = tensor[5]; + p_current[4] = tensor[4]; + p_current[5] = tensor[3]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::couple_ke() +{ + double *tensor = temperature->vector; + + if (pstyle == ISO) + mvv_current[0] = mvv_current[1] = mvv_current[2] = + tdof * boltz * t_current/dimension; + else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + mvv_current[0] = mvv_current[1] = mvv_current[2] = ave; + } else if (pcouple == XY) { + double ave = 0.5 * (tensor[0] + tensor[1]); + mvv_current[0] = mvv_current[1] = ave; + mvv_current[2] = tensor[2]; + } else if (pcouple == YZ) { + double ave = 0.5 * (tensor[1] + tensor[2]); + mvv_current[1] = mvv_current[2] = ave; + mvv_current[0] = tensor[0]; + } else if (pcouple == XZ) { + double ave = 0.5 * (tensor[0] + tensor[2]); + mvv_current[0] = mvv_current[2] = ave; + mvv_current[1] = tensor[1]; + } else { + mvv_current[0] = tensor[0]; + mvv_current[1] = tensor[1]; + mvv_current[2] = tensor[2]; + } +} + +/* ---------------------------------------------------------------------- + change box size + remap all atoms or fix group atoms depending on allremap flag + if rigid bodies exist, scale rigid body centers-of-mass +------------------------------------------------------------------------- */ + +void FixNH::remap() +{ + int i; + double oldlo,oldhi,ctr; + + double **x = atom->x; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + // convert pertinent atoms and rigid bodies to lamda coords + + if (allremap) domain->x2lamda(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->x2lamda(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(0); + + // reset global and local box to new size/shape + + if (pstyle == TRICLINIC) { + domain->yz += domain->zprd*dilation[3]; + domain->xz += domain->zprd*dilation[4]; + domain->xy += domain->yprd*dilation[5]; + + if (domain->yz < -0.5*domain->yprd || domain->yz > 0.5*domain->yprd || + domain->xz < -0.5*domain->xprd || domain->xz > 0.5*domain->xprd || + domain->xy < -0.5*domain->xprd || domain->xy > 0.5*domain->xprd) + error->all("fix npt/nph has tilted box beyond 45 degrees"); + } + + for (i = 0; i < 3; i++) { + if (p_flag[i]) { + oldlo = domain->boxlo[i]; + oldhi = domain->boxhi[i]; + ctr = 0.5 * (oldlo + oldhi); + domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr; + domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr; + } + } + + domain->set_global_box(); + domain->set_local_box(); + + // convert pertinent atoms and rigid bodies back to box coords + + if (allremap) domain->lamda2x(nlocal); + else { + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + domain->lamda2x(x[i],x[i]); + } + + if (nrigid) + for (i = 0; i < nrigid; i++) + modify->fix[rfix[i]]->deform(1); +} + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixNH::write_restart(FILE *fp) +{ + int nsize = 2; + if (tstat_flag) nsize += 1 + 2*mtchain; + if (pstat_flag) { + nsize += 16 + 2*mpchain; + if (deviatoric_flag) nsize += 6; + } + + double* list = (double *) memory->smalloc(nsize*sizeof(double),"nh:list"); + + int n = 0; + + list[n++] = tstat_flag; + if (tstat_flag) { + list[n++] = mtchain; + for (int ich = 0; ich < mtchain; ich++) + list[n++] = eta[ich]; + for (int ich = 0; ich < mtchain; ich++) + list[n++] = eta_dot[ich]; + } + + list[n++] = pstat_flag; + if (pstat_flag) { + list[n++] = omega[0]; + list[n++] = omega[1]; + list[n++] = omega[2]; + list[n++] = omega[3]; + list[n++] = omega[4]; + list[n++] = omega[5]; + list[n++] = omega_dot[0]; + list[n++] = omega_dot[1]; + list[n++] = omega_dot[2]; + list[n++] = omega_dot[3]; + list[n++] = omega_dot[4]; + list[n++] = omega_dot[5]; + list[n++] = vol0; + list[n++] = t0; + list[n++] = mpchain; + if (mpchain) { + for (int ich = 0; ich < mpchain; ich++) + list[n++] = etap[ich]; + for (int ich = 0; ich < mpchain; ich++) + list[n++] = etap_dot[ich]; + } + + list[n++] = deviatoric_flag; + if (deviatoric_flag) { + list[n++] = h0_inv[0]; + list[n++] = h0_inv[1]; + list[n++] = h0_inv[2]; + list[n++] = h0_inv[3]; + list[n++] = h0_inv[4]; + list[n++] = h0_inv[5]; + } + } + + if (comm->me == 0) { + int size = nsize * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),nsize,fp); + } + + memory->sfree(list); +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixNH::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + int flag = static_cast (list[n++]); + if (flag) { + int m = static_cast (list[n++]); + if (tstat_flag && m == mtchain) { + for (int ich = 0; ich < mtchain; ich++) + eta[ich] = list[n++]; + for (int ich = 0; ich < mtchain; ich++) + eta_dot[ich] = list[n++]; + } else n += 2*m; + } + flag = static_cast (list[n++]); + if (flag) { + omega[0] = list[n++]; + omega[1] = list[n++]; + omega[2] = list[n++]; + omega[3] = list[n++]; + omega[4] = list[n++]; + omega[5] = list[n++]; + omega_dot[0] = list[n++]; + omega_dot[1] = list[n++]; + omega_dot[2] = list[n++]; + omega_dot[3] = list[n++]; + omega_dot[4] = list[n++]; + omega_dot[5] = list[n++]; + vol0 = list[n++]; + t0 = list[n++]; + int m = static_cast (list[n++]); + if (pstat_flag && m == mpchain) { + for (int ich = 0; ich < mpchain; ich++) + etap[ich] = list[n++]; + for (int ich = 0; ich < mpchain; ich++) + etap_dot[ich] = list[n++]; + } else n+=2*m; + flag = static_cast (list[n++]); + if (flag) { + h0_inv[0] = list[n++]; + h0_inv[1] = list[n++]; + h0_inv[2] = list[n++]; + h0_inv[3] = list[n++]; + h0_inv[4] = list[n++]; + h0_inv[5] = list[n++]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNH::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"temp") == 0) { + if (narg < 2) error->all("Illegal fix_modify command"); + if (tflag) { + modify->delete_compute(id_temp); + tflag = 0; + } + delete [] id_temp; + int n = strlen(arg[1]) + 1; + id_temp = new char[n]; + strcpy(id_temp,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) error->all("Could not find fix_modify temperature ID"); + temperature = modify->compute[icompute]; + + if (temperature->tempflag == 0) + error->all("Fix_modify temperature ID does not compute temperature"); + if (temperature->igroup != 0 && comm->me == 0) + error->warning("Temperature for fix modify is not for group all"); + + // reset id_temp of pressure to new temperature ID + + if (pstat_flag) { + icompute = modify->find_compute(id_press); + if (icompute < 0) + error->all("Pressure ID for fix modify does not exist"); + modify->compute[icompute]->reset_extra_compute_fix(id_temp); + } + + return 2; + + } else if (strcmp(arg[0],"press") == 0) { + if (narg < 2) error->all("Illegal fix_modify command"); + if (!pstat_flag) error->all("Illegal fix_modify command"); + if (pflag) { + modify->delete_compute(id_press); + pflag = 0; + } + delete [] id_press; + int n = strlen(arg[1]) + 1; + id_press = new char[n]; + strcpy(id_press,arg[1]); + + int icompute = modify->find_compute(arg[1]); + if (icompute < 0) error->all("Could not find fix_modify pressure ID"); + pressure = modify->compute[icompute]; + + if (pressure->pressflag == 0) + error->all("Fix_modify pressure ID does not compute pressure"); + return 2; + } + + return 0; +} + +/* ---------------------------------------------------------------------- */ + +double FixNH::compute_scalar() +{ + int i; + double volume; + double energy; + double kt = boltz * t_target; + double lkt = tdof * kt; + double lkt_press = kt; + int ich; + if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; + else volume = domain->xprd * domain->yprd; + + energy = 0.0; + + // thermostat chain energy is equivalent to Eq. (2) in + // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 + // Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M), + // where L = tdof + // M = mtchain + // p_eta_k = Q_k*eta_dot[k-1] + // Q_1 = L*k*T/t_freq^2 + // Q_k = k*T/t_freq^2, k > 1 + + if (tstat_flag) { + energy += lkt * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; + for (ich = 1; ich < mtchain; ich++) + energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; + } + + if (pstat_flag) { + // barostat energy is equivalent to Eq. (8) in + // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117 + // Sum(0.5*p_omega^2/W + P*V), + // where N = natoms + // p_omega = W*omega_dot + // W = N*k*T/p_freq^2 + // sum is over barostatted dimensions + + for (i = 0; i < 3; i++) + if (p_flag[i]) + energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] + + p_hydro*(volume-vol0) / (pdim*nktv2p); + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) + energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i]; + } + + // extra contributions from thermostat chain for barostat + + if (mpchain) { + energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; + for (ich = 1; ich < mpchain; ich++) + energy += kt * etap[ich] + + 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; + } + + // extra contribution from strain energy + + if (deviatoric_flag) energy += compute_strain_energy(); + } + + return energy; +} + +/* ---------------------------------------------------------------------- + return a single element of the following vectors, in this order: + eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof] + etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain] + PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain] + PE_strain[1] + if no thermostat exists, related quantities are omitted from the list + if no barostat exists, related quantities are omitted from the list + ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI +------------------------------------------------------------------------- */ + +double FixNH::compute_vector(int n) +{ + int ilen; + + if (tstat_flag) { + ilen = mtchain; + if (n < ilen) return eta[n]; + n -= ilen; + ilen = mtchain; + if (n < ilen) return eta_dot[n]; + n -= ilen; + } + + if (pstat_flag) { + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) return omega[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) return omega[n]; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) return omega[n]; + n -= ilen; + } + + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) return omega_dot[n]; + n -= ilen; + } + + if (mpchain) { + ilen = mpchain; + if (n < ilen) return etap[n]; + n -= ilen; + ilen = mpchain; + if (n < ilen) return etap_dot[n]; + n -= ilen; + } + } + + int i; + double volume; + double kt = boltz * t_target; + double lkt = tdof * kt; + double lkt_press = kt; + int ich; + if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; + else volume = domain->xprd * domain->yprd; + + if (tstat_flag) { + ilen = mtchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return lkt * eta[0]; + else + return kt * eta[ich]; + } + n -= ilen; + ilen = mtchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0]; + else + return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich]; + } + n -= ilen; + } + + if (pstat_flag) { + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) + return p_hydro*(volume-vol0) / nktv2p; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) + if (p_flag[n]) + return p_hydro*(volume-vol0) / (pdim*nktv2p); + else + return 0.0; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) + if (n > 2) return 0.0; + else if (p_flag[n]) + return p_hydro*(volume-vol0) / (pdim*nktv2p); + else + return 0.0; + n -= ilen; + } + + if (pstyle == ISO) { + ilen = 1; + if (n < ilen) + return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + n -= ilen; + } else if (pstyle == ANISO) { + ilen = 3; + if (n < ilen) + if (p_flag[n]) + return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + else return 0.0; + n -= ilen; + } else { + ilen = 6; + if (n < ilen) + if (p_flag[n]) + return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n]; + else return 0.0; + n -= ilen; + } + + if (mpchain) { + ilen = mpchain; + if (n < ilen) { + ich = n; + if (ich == 0) return lkt_press * etap[0]; + else return kt * etap[ich]; + } + n -= ilen; + ilen = mpchain; + if (n < ilen) { + ich = n; + if (ich == 0) + return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0]; + else + return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich]; + } + n -= ilen; + } + + if (deviatoric_flag) { + ilen = 1; + if (n < ilen) + return compute_strain_energy(); + n -= ilen; + } + } + + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixNH::reset_dt() +{ + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + dthalf = 0.5 * update->dt; + dt4 = 0.25 * update->dt; + dt8 = 0.125 * update->dt; + + p_freq_max = 0.0; + if (pstat_flag) { + p_freq_max = MAX(p_freq[0],p_freq[1]); + p_freq_max = MAX(p_freq_max,p_freq[2]); + if (pstyle == TRICLINIC) { + p_freq_max = MAX(p_freq_max,p_freq[3]); + p_freq_max = MAX(p_freq_max,p_freq[4]); + p_freq_max = MAX(p_freq_max,p_freq[5]); + } + drag_factor = 1.0 - (update->dt * p_freq_max * drag); + } + + if (tstat_flag) + drag_factor = 1.0 - (update->dt * MAX(p_freq_max,t_freq) * drag); +} + +/* ---------------------------------------------------------------------- + perform half-step update of chain thermostat variables +------------------------------------------------------------------------- */ + +void FixNH::nhc_temp_integrate() +{ + int ich; + double expfac; + + double lkt = tdof * boltz * t_target; + double kecurrent = tdof * boltz * t_current; + eta_dotdot[0] = (kecurrent - lkt)/eta_mass[0]; + + double ncfac = 1.0/nc_tchain; + for (int iloop = 0; iloop < nc_tchain; iloop++) { + + for (ich = mtchain-1; ich > 0; ich--) { + expfac = exp(-ncfac*dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; + eta_dot[ich] *= drag_factor; + eta_dot[ich] *= expfac; + } + + expfac = exp(-ncfac*dt8*eta_dot[1]); + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0] * ncfac*dt4; + eta_dot[0] *= drag_factor; + eta_dot[0] *= expfac; + + factor_eta = exp(-ncfac*dthalf*eta_dot[0]); + nh_v_temp(); + + // rescale temperature due to velocity scaling + // should not be necessary to explicitly recompute the temperature + + t_current *= factor_eta*factor_eta; + kecurrent = tdof * boltz * t_current; + eta_dotdot[0] = (kecurrent - lkt)/eta_mass[0]; + + for (ich = 0; ich < mtchain; ich++) + eta[ich] += ncfac*dthalf*eta_dot[ich]; + + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0] * ncfac*dt4; + eta_dot[0] *= expfac; + + for (ich = 1; ich < mtchain; ich++) { + expfac = exp(-ncfac*dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] + - boltz * t_target)/eta_mass[ich]; + eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4; + eta_dot[ich] *= expfac; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of chain thermostat variables for barostat + scale barostat velocities +------------------------------------------------------------------------- */ + +void FixNH::nhc_press_integrate() +{ + int ich,i; + double expfac,factor_etap,wmass,kecurrent; + double kt = boltz * t_target; + double lkt_press = kt; + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + } + + etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; + + double ncfac = 1.0/nc_pchain; + for (int iloop = 0; iloop < nc_pchain; iloop++) { + + for (ich = mpchain-1; ich > 0; ich--) { + expfac = exp(-ncfac*dt8*etap_dot[ich+1]); + etap_dot[ich] *= expfac; + etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; + etap_dot[ich] *= drag_factor; + etap_dot[ich] *= expfac; + } + + expfac = exp(-ncfac*dt8*etap_dot[1]); + etap_dot[0] *= expfac; + etap_dot[0] += etap_dotdot[0] * ncfac*dt4; + etap_dot[0] *= drag_factor; + etap_dot[0] *= expfac; + + for (ich = 0; ich < mpchain; ich++) + etap[ich] += ncfac*dthalf*etap_dot[ich]; + + factor_etap = exp(-ncfac*dthalf*etap_dot[0]); + for (i = 0; i < 3; i++) + if (p_flag[i]) omega_dot[i] *= factor_etap; + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) omega_dot[i] *= factor_etap; + } + + kecurrent = 0.0; + for (i = 0; i < 3; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + + if (pstyle == TRICLINIC) { + for (i = 3; i < 6; i++) + if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i]; + } + + etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; + + etap_dot[0] *= expfac; + etap_dot[0] += etap_dotdot[0] * ncfac*dt4; + etap_dot[0] *= expfac; + + for (ich = 1; ich < mpchain; ich++) { + expfac = exp(-ncfac*dt8*etap_dot[ich+1]); + etap_dot[ich] *= expfac; + etap_dotdot[ich] = + (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) / + etap_mass[ich]; + etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4; + etap_dot[ich] *= expfac; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step barostat scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNH::nh_v_press() +{ + double **v = atom->v; + int *mask = atom->mask; + 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) { + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + if (pstyle == TRICLINIC) { + v[i][0] += v[i][1]*factor[5] + v[i][2]*factor[4]; + v[i][1] += v[i][2]*factor[3]; + } + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + v[i][0] *= factor[0]; + v[i][1] *= factor[1]; + v[i][2] *= factor[2]; + if (pstyle == TRICLINIC) { + v[i][0] += v[i][1]*factor[5] + v[i][2]*factor[4]; + v[i][1] += v[i][2]*factor[3]; + } + temperature->restore_bias(i,v[i]); + } + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of velocities +-----------------------------------------------------------------------*/ + +void FixNH::nve_v() +{ + double dtfm; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; + } + } + } else { + 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]; + } + } + } +} + +/* ---------------------------------------------------------------------- + perform full-step update of positions +-----------------------------------------------------------------------*/ + +void FixNH::nve_x() +{ + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // x update by full step only for atoms in group + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step thermostat scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNH::nh_v_temp() +{ + double **v = atom->v; + int *mask = atom->mask; + 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) { + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + } + } + } else if (which == BIAS) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + temperature->remove_bias(i,v[i]); + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + temperature->restore_bias(i,v[i]); + } + } + } +} + +/* ---------------------------------------------------------------------- + compute sigma tensor + needed whenever p_target or h0_inv changes +-----------------------------------------------------------------------*/ + +void FixNH::compute_sigma() +{ + // if nreset_h0 > 0, reset vol0 and h0_inv + // every nreset_h0 timesteps + + if (nreset_h0 > 0) { + int delta = update->ntimestep - update->beginstep; + if (delta % nreset_h0 == 0) { + if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; + else vol0 = domain->xprd * domain->yprd; + h0_inv[0] = domain->h_inv[0]; + h0_inv[1] = domain->h_inv[1]; + h0_inv[2] = domain->h_inv[2]; + h0_inv[3] = domain->h_inv[3]; + h0_inv[4] = domain->h_inv[4]; + h0_inv[5] = domain->h_inv[5]; + } + } + + // generate upper-triangular half of + // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t + // units of sigma are are PV/L^2 e.g. atm.A + // + // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] + // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] + // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] + + sigma[0] = + vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] + + p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) + + h0_inv[5]*(p_target[5]*h0_inv[0] + + (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) + + h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] + + (p_target[2]-p_hydro)*h0_inv[4])); + sigma[1] = + vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] + + p_target[3]*h0_inv[3]) + + h0_inv[3]*(p_target[3]*h0_inv[1] + + (p_target[2]-p_hydro)*h0_inv[3])); + sigma[2] = + vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[3] = + vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) + + h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[4] = + vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) + + h0_inv[5]*(p_target[3]*h0_inv[2]) + + h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2])); + sigma[5] = + vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) + + h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) + + h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3])); +} + +/* ---------------------------------------------------------------------- + compute strain energy +-----------------------------------------------------------------------*/ + +double FixNH::compute_strain_energy() +{ + // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units + + double* h = domain->h; + double d0,d1,d2; + + d0 = + sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) + + sigma[5]*( h[1]*h[5]+h[3]*h[4]) + + sigma[4]*( h[2]*h[4]); + d1 = + sigma[5]*( h[5]*h[1]+h[4]*h[3]) + + sigma[1]*( h[1]*h[1]+h[3]*h[3]) + + sigma[3]*( h[2]*h[3]); + d2 = + sigma[4]*( h[4]*h[2]) + + sigma[3]*( h[3]*h[2]) + + sigma[2]*( h[2]*h[2]); + + double energy = 0.5*(d0+d1+d2)/nktv2p; + return energy; +} + +/* ---------------------------------------------------------------------- + compute deviatoric barostat force = h*sigma*h^t +-----------------------------------------------------------------------*/ + +void FixNH::compute_deviatoric() +{ + // generate upper-triangular part of h*sigma*h^t + // units of fdev are are PV, e.g. atm*A^3 + // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ] + // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ] + // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ] + + double* h = domain->h; + + fdev[0] = + h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) + + h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) + + h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]); + fdev[1] = + h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) + + h[3]*( sigma[3]*h[1]+sigma[2]*h[3]); + fdev[2] = + h[2]*( sigma[2]*h[2]); + fdev[3] = + h[1]*( sigma[3]*h[2]) + + h[3]*( sigma[2]*h[2]); + fdev[4] = + h[0]*( sigma[4]*h[2]) + + h[5]*( sigma[3]*h[2]) + + h[4]*( sigma[2]*h[2]); + fdev[5] = + h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) + + h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) + + h[4]*( sigma[3]*h[1]+sigma[2]*h[3]); +} + +/* ---------------------------------------------------------------------- + compute hydrostatic target pressure +-----------------------------------------------------------------------*/ + +void FixNH::compute_press_target() +{ + double delta = update->ntimestep - update->beginstep; + if (update->endstep > update->beginstep) + delta /= update->endstep - update->beginstep; + else delta = 0.0; + + p_hydro = 0.0; + for (int i = 0; i < 3; i++) + if (p_flag[i]) { + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + p_hydro += p_target[i]; + } + p_hydro /= pdim; + + if (pstyle == TRICLINIC) + for (int i = 3; i < 6; i++) + p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); + + // if deviatoric, recompute sigma each time p_target changes + + if (deviatoric_flag) compute_sigma(); +} + +/* ---------------------------------------------------------------------- + update omega_dot, omega, dilation +-----------------------------------------------------------------------*/ + +void FixNH::nh_omega_dot() +{ + double mtk_term; + double f_omega,volume; + if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; + else volume = domain->xprd*domain->yprd; + + if (deviatoric_flag) compute_deviatoric(); + + for (int i = 0; i < 3; i++) { + if (p_flag[i]) { + if (mtk_flag) mtk_term = mvv_current[i]/(atom->natoms*volume) * nktv2p; + else mtk_term = 0.0; + f_omega = (p_current[i]-p_hydro+mtk_term)*volume / + (omega_mass[i] * nktv2p); + if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p); + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= drag_factor; + } + omega[i] += dthalf*omega_dot[i]; + factor[i] = exp(-dthalf*omega_dot[i]*mtk_factor); + dilation[i] = exp(dthalf*omega_dot[i]); + } + + if (pstyle == TRICLINIC) { + for (int i = 3; i < 6; i++) { + if (p_flag[i]) { + f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p); + if (deviatoric_flag) + f_omega -= fdev[i]/(omega_mass[i] * nktv2p); + omega_dot[i] += f_omega*dthalf; + omega_dot[i] *= drag_factor; + } + omega[i] += dthalf*omega_dot[i]; + factor[i] = -dthalf*omega_dot[i]; + dilation[i] = dthalf*omega_dot[i]; + } + } +} diff --git a/src/fix_nh.h b/src/fix_nh.h new file mode 100644 index 0000000000..7b4966f905 --- /dev/null +++ b/src/fix_nh.h @@ -0,0 +1,119 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_NH_H +#define LMP_FIX_NH_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNH : public Fix { + public: + FixNH(class LAMMPS *, int, char **); + virtual ~FixNH(); + int setmask(); + virtual void init(); + void setup(int); + virtual void initial_integrate(int); + virtual void final_integrate(); + void initial_integrate_respa(int, int, int); + void final_integrate_respa(int, int); + double compute_scalar(); + double compute_vector(int); + void write_restart(FILE *); + void restart(char *); + int modify_param(int, char **); + void reset_dt(); + + protected: + int dimension,which; + double dtv,dtf,dthalf,dt4,dt8; + double boltz,nktv2p,tdof; + double vol0,t0; + + double t_start,t_stop; + double t_current,t_target; + double t_freq; + + int tstat_flag; // 1 if control T + int pstat_flag; // 1 if control P + + int pstyle,pcouple,allremap; + int p_flag[6]; // 1 if control P on this dim, 0 if not + double p_start[6],p_stop[6]; + double p_freq[6],p_target[6]; + double omega[6],omega_dot[6]; + double omega_mass[6]; + double p_current[6],dilation[6]; + double drag,drag_factor; // drag factor on particle thermostat + double factor[6]; // velocity scaling due to barostat + int kspace_flag; // 1 if KSpace invoked, 0 if not + int nrigid; // number of rigid fixes + int *rfix; // indices of rigid fixes + + int nlevels_respa; + double *step_respa; + + char *id_temp,*id_press; + class Compute *temperature,*pressure; + int tflag,pflag; + + double *eta,*eta_dot; // chain thermostat for particles + double *eta_dotdot; + double *eta_mass; + int mtchain; // length of chain + + double *etap; // chain thermostat for barostat + double *etap_dot; + double *etap_dotdot; + double *etap_mass; + int mpchain; // length of chain + + int mtk_flag; // 0 if using Hoover barostat + int pdim; // number of barostatted dims + double mvv_current[3]; // diagonal of KE tensor + double mtk_factor; // MTK factor + double p_freq_max; // maximum barostat frequency + + double p_hydro; // hydrostatic target pressure + + int nc_tchain,nc_pchain; + double factor_eta; + double sigma[6]; // scaled target stress + double fdev[6]; // deviatoric force on barostat + int deviatoric_flag; // 0 if target stress tensor is hydrostatic + double h0_inv[6]; // h_inv of reference (zero strain) box + int nreset_h0; // interval for resetting h0 + + void couple(); + void couple_ke(); + void remap(); + void nhc_temp_integrate(); + void nhc_press_integrate(); + + virtual void nve_x(); // may be overwritten by child classes + virtual void nve_v(); + virtual void nh_v_press(); + virtual void nh_v_temp(); + + void compute_sigma(); + void compute_deviatoric(); + double compute_strain_energy(); + void compute_press_target(); + void nh_omega_dot(); +}; + +} + +#endif diff --git a/src/fix_nh_asphere.cpp b/src/fix_nh_asphere.cpp new file mode 100644 index 0000000000..f832918c7e --- /dev/null +++ b/src/fix_nh_asphere.cpp @@ -0,0 +1,238 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "math.h" +#include "math_extra.h" +#include "fix_nh_asphere.h" +#include "atom.h" +#include "atom_vec.h" +#include "group.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +enum{NOBIAS,BIAS}; + +/* ---------------------------------------------------------------------- */ + +FixNHAsphere::FixNHAsphere(LAMMPS *lmp, int narg, char **arg) : + FixNH(lmp, narg, arg) +{ + inertia = + memory->create_2d_double_array(atom->ntypes+1,3,"fix_nvt_asphere:inertia"); + + if (!atom->quat_flag || !atom->angmom_flag || !atom->torque_flag || + !atom->avec->shape_type) + error->all("Fix nvt/nph/npt asphere requires atom attributes " + "quat, angmom, torque, shape"); + if (atom->radius_flag || atom->rmass_flag) + error->all("Fix nvt/nph/npt asphere cannot be used with atom attributes " + "diameter or rmass"); +} + +/* ---------------------------------------------------------------------- */ + +FixNHAsphere::~FixNHAsphere() +{ + memory->destroy_2d_double_array(inertia); +} + +/* ---------------------------------------------------------------------- */ + +void FixNHAsphere::init() +{ + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (shape[type[i]][0] == 0.0) + error->one("Fix nvt/nph/npt asphere requires extended particles"); + + FixNH::init(); + calculate_inertia(); +} + +/* ---------------------------------------------------------------------- + Richardson iteration to update quaternion accurately +------------------------------------------------------------------------- */ + +void FixNHAsphere::richardson(double *q, double *m, double *moments) +{ + // compute omega at 1/2 step from m at 1/2 step and q at 0 + + double w[3]; + omega_from_mq(q,m,moments,w); + + // full update from dq/dt = 1/2 w q + + double wq[4]; + MathExtra::multiply_vec_quat(w,q,wq); + + double qfull[4]; + qfull[0] = q[0] + dtq * wq[0]; + qfull[1] = q[1] + dtq * wq[1]; + qfull[2] = q[2] + dtq * wq[2]; + qfull[3] = q[3] + dtq * wq[3]; + MathExtra::normalize4(qfull); + + // 1st half of update from dq/dt = 1/2 w q + + double qhalf[4]; + qhalf[0] = q[0] + 0.5*dtq * wq[0]; + qhalf[1] = q[1] + 0.5*dtq * wq[1]; + qhalf[2] = q[2] + 0.5*dtq * wq[2]; + qhalf[3] = q[3] + 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step + // recompute wq + + omega_from_mq(qhalf,m,moments,w); + MathExtra::multiply_vec_quat(w,qhalf,wq); + + // 2nd half of update from dq/dt = 1/2 w q + + qhalf[0] += 0.5*dtq * wq[0]; + qhalf[1] += 0.5*dtq * wq[1]; + qhalf[2] += 0.5*dtq * wq[2]; + qhalf[3] += 0.5*dtq * wq[3]; + MathExtra::normalize4(qhalf); + + // corrected Richardson update + + q[0] = 2.0*qhalf[0] - qfull[0]; + q[1] = 2.0*qhalf[1] - qfull[1]; + q[2] = 2.0*qhalf[2] - qfull[2]; + q[3] = 2.0*qhalf[3] - qfull[3]; + MathExtra::normalize4(q); +} + +/* ---------------------------------------------------------------------- + compute omega from angular momentum + w = omega = angular velocity in space frame + wbody = angular velocity in body frame + project space-frame angular momentum onto body axes + and divide by principal moments +------------------------------------------------------------------------- */ + +void FixNHAsphere::omega_from_mq(double *q, double *m, double *inertia, + double *w) +{ + double rot[3][3]; + MathExtra::quat_to_mat(q,rot); + + double wbody[3]; + MathExtra::transpose_times_column3(rot,m,wbody); + wbody[0] /= inertia[0]; + wbody[1] /= inertia[1]; + wbody[2] /= inertia[2]; + MathExtra::times_column3(rot,wbody,w); +} + +/* ---------------------------------------------------------------------- + principal moments of inertia for ellipsoids +------------------------------------------------------------------------- */ + +void FixNHAsphere::calculate_inertia() +{ + double *mass = atom->mass; + double **shape = atom->shape; + + for (int i = 1; i <= atom->ntypes; i++) { + inertia[i][0] = 0.2*mass[i] * + (shape[i][1]*shape[i][1]+shape[i][2]*shape[i][2]); + inertia[i][1] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][2]*shape[i][2]); + inertia[i][2] = 0.2*mass[i] * + (shape[i][0]*shape[i][0]+shape[i][1]*shape[i][1]); + } +} + +/* ---------------------------------------------------------------------- + perform half-step update of velocities +-----------------------------------------------------------------------*/ + +void FixNHAsphere::nve_v() +{ + // standard nhc_nve_v velocity update + + FixNH::nve_v(); + + int *type = atom->type; + double **quat = atom->quat; + double **angmom = atom->angmom; + double **torque = atom->torque; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + + dtq = 0.5 * dtv; + + // update angular momentum by 1/2 step for all particles + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + angmom[i][0] += dtf*torque[i][0]; + angmom[i][1] += dtf*torque[i][1]; + angmom[i][2] += dtf*torque[i][2]; + + richardson(quat[i],angmom[i],inertia[type[i]]); + } + } +} + +/* ---------------------------------------------------------------------- + perform half-step temperature scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNHAsphere::nh_v_temp() +{ + // standard nhc_nh_v velocity update + + FixNH::nh_v_temp(); + + double **angmom = atom->angmom; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + + double factor_rotate = exp(-dthalf*eta_dot[0]); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + angmom[i][0] *= factor_rotate; + angmom[i][1] *= factor_rotate; + angmom[i][2] *= factor_rotate; + } + } +} diff --git a/src/fix_nh_asphere.h b/src/fix_nh_asphere.h new file mode 100644 index 0000000000..6496066541 --- /dev/null +++ b/src/fix_nh_asphere.h @@ -0,0 +1,40 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_NH_ASPHERE_H +#define LMP_FIX_NH_ASPHERE_H + +#include "fix_nh.h" + +namespace LAMMPS_NS { + +class FixNHAsphere : public FixNH { + public: + FixNHAsphere(class LAMMPS *, int, char **); + virtual ~FixNHAsphere(); + void init(); + + protected: + double dtq; + double **inertia; + + void richardson(double *, double *, double *); + void omega_from_mq(double *, double *, double *, double *); + void calculate_inertia(); + void nve_v(); + void nh_v_temp(); +}; + +} + +#endif diff --git a/src/fix_nh_sphere.cpp b/src/fix_nh_sphere.cpp new file mode 100644 index 0000000000..29622a13f5 --- /dev/null +++ b/src/fix_nh_sphere.cpp @@ -0,0 +1,191 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Mike Brown (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "fix_nh_sphere.h" +#include "atom.h" +#include "atom_vec.h" +#include "group.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define INERTIA 0.4 // moment of inertia for sphere + +enum{NOBIAS,BIAS}; + +/* ---------------------------------------------------------------------- */ + +FixNHSphere::FixNHSphere(LAMMPS *lmp, int narg, char **arg) : + FixNH(lmp, narg, arg) +{ + if (!atom->omega_flag || !atom->torque_flag) + error->all("Fix nvt/nph/npt sphere requires " + "atom attributes omega, torque"); + if (!atom->radius_flag && !atom->avec->shape_type) + error->all("Fix nvt/nph/npt sphere requires " + "atom attribute radius or shape"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNHSphere::init() +{ + int i,itype; + + // check that all particles are finite-size and spherical + // no point particles allowed + + if (atom->radius_flag) { + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (radius[i] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + } + + } else { + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + itype = type[i]; + if (shape[itype][0] == 0.0) + error->one("Fix nvt/sphere requires extended particles"); + if (shape[itype][0] != shape[itype][1] || + shape[itype][0] != shape[itype][2]) + error->one("Fix nvt/sphere requires spherical particle shapes"); + } + } + + FixNH::init(); +} + +/* ---------------------------------------------------------------------- + perform half-step update of velocities +-----------------------------------------------------------------------*/ + +void FixNHSphere::nve_v() +{ + // standard nve_v velocity update + + FixNH::nve_v(); + + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + double *mass = atom->mass; + double **shape = atom->shape; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + + double dtfrotate = dtf / INERTIA; + double dtirotate; + int itype; + + // update omega for all particles + // d_omega/dt = torque / inertia + // 4 cases depending on radius vs shape and rmass vs mass + + if (radius) { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); + omega[i][0] += dtirotate*torque[i][0]; + omega[i][1] += dtirotate*torque[i][1]; + omega[i][2] += dtirotate*torque[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]); + omega[i][0] += dtirotate*torque[i][0]; + omega[i][1] += dtirotate*torque[i][1]; + omega[i][2] += dtirotate*torque[i][2]; + } + } + } + + } else { + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); + omega[i][0] += dtirotate*torque[i][0]; + omega[i][1] += dtirotate*torque[i][1]; + omega[i][2] += dtirotate*torque[i][2]; + } + } + } else { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + itype = type[i]; + dtirotate = dtfrotate / + (shape[itype][0]*shape[itype][0]*mass[itype]); + omega[i][0] += dtirotate*torque[i][0]; + omega[i][1] += dtirotate*torque[i][1]; + omega[i][2] += dtirotate*torque[i][2]; + } + } + } + } + +} + +/* ---------------------------------------------------------------------- + perform half-step scaling of velocities +-----------------------------------------------------------------------*/ + +void FixNHSphere::nh_v_temp() +{ + // standard nh_v_temp velocity update + + FixNH::nh_v_temp(); + + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + + double factor_rotate = exp(-dthalf*eta_dot[0]); + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + omega[i][0] *= factor_rotate; + omega[i][1] *= factor_rotate; + omega[i][2] *= factor_rotate; + } + } +} diff --git a/src/fix_nh_sphere.h b/src/fix_nh_sphere.h new file mode 100644 index 0000000000..9a3a57b845 --- /dev/null +++ b/src/fix_nh_sphere.h @@ -0,0 +1,34 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifndef LMP_FIX_NH_SPHERE_H +#define LMP_FIX_NH_SPHERE_H + +#include "fix_nh.h" + +namespace LAMMPS_NS { + +class FixNHSphere : public FixNH { + public: + FixNHSphere(class LAMMPS *, int, char **); + virtual ~FixNHSphere() {} + void init(); + + protected: + void nve_v(); + void nh_v_temp(); +}; + +} + +#endif diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp index 5169ce2943..838a2091e4 100644 --- a/src/fix_nph.cpp +++ b/src/fix_nph.cpp @@ -11,182 +11,39 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing author: Mark Stevens (SNL) -------------------------------------------------------------------------- */ - #include "string.h" -#include "stdlib.h" -#include "math.h" #include "fix_nph.h" -#include "atom.h" -#include "force.h" -#include "comm.h" -#include "domain.h" +#include "group.h" #include "modify.h" -#include "fix_deform.h" -#include "compute.h" -#include "kspace.h" -#include "update.h" -#include "respa.h" -#include "domain.h" #include "error.h" using namespace LAMMPS_NS; -#define MIN(A,B) ((A) < (B)) ? (A) : (B) -#define MAX(A,B) ((A) > (B)) ? (A) : (B) - -enum{XYZ,XY,YZ,XZ,ANISO}; - /* ---------------------------------------------------------------------- */ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + FixNH(lmp, narg, arg) { - if (narg < 4) error->all("Illegal fix nph command"); - - restart_global = 1; - box_change = 1; - time_integrate = 1; - scalar_flag = 1; - global_freq = 1; - extscalar = 1; - - double p_period[3]; - if (strcmp(arg[3],"xyz") == 0) { - if (narg < 7) error->all("Illegal fix nph command"); - - press_couple = XYZ; - p_start[0] = p_start[1] = p_start[2] = atof(arg[4]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[6]); - p_flag[0] = p_flag[1] = p_flag[2] = 1; - if (domain->dimension == 2) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } - - } else { - if (strcmp(arg[3],"xy") == 0) press_couple = XY; - else if (strcmp(arg[3],"yz") == 0) press_couple = YZ; - else if (strcmp(arg[3],"xz") == 0) press_couple = XZ; - else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO; - else error->all("Illegal fix nph command"); - - if (narg < 11) error->all("Illegal fix nph command"); - - if (domain->dimension == 2 && - (press_couple == XY || press_couple == YZ || press_couple == XZ)) - error->all("Invalid fix nph command for a 2d simulation"); - - if (strcmp(arg[4],"NULL") == 0) { - p_start[0] = p_stop[0] = p_period[0] = 0.0; - p_flag[0] = 0; - } else { - p_start[0] = atof(arg[4]); - p_stop[0] = atof(arg[5]); - p_flag[0] = 1; - } - if (strcmp(arg[6],"NULL") == 0) { - p_start[1] = p_stop[1] = p_period[1] = 0.0; - p_flag[1] = 0; - } else { - p_start[1] = atof(arg[6]); - p_stop[1] = atof(arg[7]); - p_flag[1] = 1; - } - if (strcmp(arg[8],"NULL") == 0) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } else { - if (domain->dimension == 2) - error->all("Invalid fix nph command for a 2d simulation"); - p_start[2] = atof(arg[8]); - p_stop[2] = atof(arg[9]); - p_flag[2] = 1; - } - - double period = atof(arg[10]); - if (p_flag[0]) p_period[0] = period; - if (p_flag[1]) p_period[1] = period; - if (p_flag[2]) p_period[2] = period; - } - - // process extra keywords - - drag = 0.0; - allremap = 1; - - int iarg; - if (press_couple == XYZ) iarg = 7; - else iarg = 11; - - while (iarg < narg) { - if (strcmp(arg[iarg],"drag") == 0) { - if (iarg+2 > narg) error->all("Illegal fix nph command"); - drag = atof(arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"dilate") == 0) { - if (iarg+2 > narg) error->all("Illegal fix nph command"); - if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; - else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0; - else error->all("Illegal fix nph command"); - iarg += 2; - } else error->all("Illegal fix nph command"); - } - - // error checks - - if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Invalid fix nph command pressure settings"); - if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) - error->all("Invalid fix nph command pressure settings"); - if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) - error->all("Invalid fix nph command pressure settings"); - - if (press_couple == XY && - (p_start[0] != p_start[1] || p_stop[0] != p_stop[1])) - error->all("Invalid fix nph command pressure settings"); - if (press_couple == YZ && - (p_start[1] != p_start[2] || p_stop[1] != p_stop[2])) - error->all("Invalid fix nph command pressure settings"); - if (press_couple == XZ && - (p_start[0] != p_start[2] || p_stop[0] != p_stop[2])) - error->all("Invalid fix nph command pressure settings"); - - if (p_flag[0] && domain->xperiodic == 0) - error->all("Cannot use fix nph on a non-periodic dimension"); - if (p_flag[1] && domain->yperiodic == 0) - error->all("Cannot use fix nph on a non-periodic dimension"); - if (p_flag[2] && domain->zperiodic == 0) - error->all("Cannot use fix nph on a non-periodic dimension"); - - // convert input periods to frequencies - - if ((p_flag[0] && p_period[0] <= 0.0) || - (p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0)) - error->all("Fix nph periods must be > 0.0"); - - p_freq[0] = p_freq[1] = p_freq[2] = 0.0; - if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; - if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; - if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + if (tstat_flag) + error->all("Temperature control can not be used with fix nph"); + if (!pstat_flag) + error->all("Pressure control must be used with fix nph"); // create a new compute temp style // id = fix-ID + temp // compute group = all since pressure is always global (group all) - // and thus its KE/temperature contribution should use group all + // and thus its KE/temperature contribution should use group all int n = strlen(id) + 6; id_temp = new char[n]; strcpy(id_temp,id); strcat(id_temp,"_temp"); - + char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = (char *) "all"; newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); delete [] newarg; tflag = 1; @@ -199,7 +56,7 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : id_press = new char[n]; strcpy(id_press,id); strcat(id_press,"_press"); - + newarg = new char*[4]; newarg[0] = id_press; newarg[1] = (char *) "all"; @@ -208,659 +65,4 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(4,newarg); delete [] newarg; pflag = 1; - - // Nose/Hoover pressure init - - omega[0] = omega[1] = omega[2] = 0.0; - omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0; - - nrigid = 0; - rfix = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixNPH::~FixNPH() -{ - delete [] rfix; - - // delete temperature and pressure if fix created them - - if (tflag) modify->delete_compute(id_temp); - if (pflag) modify->delete_compute(id_press); - delete [] id_temp; - delete [] id_press; -} - -/* ---------------------------------------------------------------------- */ - -int FixNPH::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= FINAL_INTEGRATE; - mask |= THERMO_ENERGY; - mask |= INITIAL_INTEGRATE_RESPA; - mask |= FINAL_INTEGRATE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixNPH::init() -{ - if (domain->triclinic) error->all("Cannot use fix nph with triclinic box"); - - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"deform") == 0) { - int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; - if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || - (p_flag[2] && dimflag[2])) - error->all("Cannot use fix npt and fix deform on same dimension"); - } - - // set temperature and pressure ptrs - - int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all("Temperature ID for fix nph does not exist"); - temperature = modify->compute[icompute]; - - icompute = modify->find_compute(id_press); - if (icompute < 0) error->all("Pressure ID for fix nph does not exist"); - pressure = modify->compute[icompute]; - - // set timesteps and frequencies - - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; - dthalf = 0.5 * update->dt; - - double freq = MAX(p_freq[0],p_freq[1]); - freq = MAX(freq,p_freq[2]); - drag_factor = 1.0 - (update->dt * freq * drag); - - boltz = force->boltz; - nktv2p = force->nktv2p; - dimension = domain->dimension; - if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; - else vol0 = domain->xprd * domain->yprd; - - if (force->kspace) kspace_flag = 1; - else kspace_flag = 0; - - if (strcmp(update->integrate_style,"respa") == 0) { - nlevels_respa = ((Respa *) update->integrate)->nlevels; - step_respa = ((Respa *) update->integrate)->step; - } - - // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - - delete [] rfix; - nrigid = 0; - rfix = NULL; - - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; - } -} - -/* ---------------------------------------------------------------------- - compute T,P before integrator starts -------------------------------------------------------------------------- */ - -void FixNPH::setup(int vflag) -{ - // Nkt = initial value for piston mass and energy conservation - // cannot be done in init() b/c temperature cannot be called there - // is b/c Modify::init() inits computes after fixes due to dof dependence - // guesstimate a unit-dependent t_initial if actual T = 0.0 - - double t_initial = temperature->compute_scalar(); - if (t_initial == 0.0) { - if (strcmp(update->unit_style,"lj") == 0) t_initial = 1.0; - else t_initial = 300.0; - } - nkt = atom->natoms * boltz * t_initial; - - p_target[0] = p_start[0]; // used by compute_scalar() - p_target[1] = p_start[1]; - p_target[2] = p_start[2]; - - if (press_couple == XYZ) { - double tmp = temperature->compute_scalar(); - tmp = pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - - // trigger virial computation on next timestep - - pressure->addstep(update->ntimestep+1); -} - -/* ---------------------------------------------------------------------- - 1st half of Verlet update -------------------------------------------------------------------------- */ - -void FixNPH::initial_integrate(int vflag) -{ - int i; - double dtfm; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = nkt / volume * nktv2p; - - for (i = 0; i < 3; i++) { - p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - omega[i] += dtv*omega_dot[i]; - factor[i] = exp(-dthalf*omega_dot[i]); - dilation[i] = exp(dthalf*omega_dot[i]); - } - - // v update only for atoms in NPH group - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (rmass) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - 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]; - } - } - } - - // remap simulation box by 1/2 step - - remap(); - - // x update by full step only for atoms in NPH group - - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - - // remap simulation box by 1/2 step - // redo KSpace coeffs since volume has changed - - remap(); - if (kspace_flag) force->kspace->setup(); -} - -/* ---------------------------------------------------------------------- - 2nd half of Verlet update -------------------------------------------------------------------------- */ - -void FixNPH::final_integrate() -{ - int i; - double dtfm; - - // v update only for atoms in NPH group - - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (rmass) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - 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]; - } - } - } - - // compute new pressure - - if (press_couple == XYZ) { - double tmp = temperature->compute_scalar(); - tmp = pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - - // trigger virial computation on next timestep - - pressure->addstep(update->ntimestep+1); - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = nkt / volume * nktv2p; - - for (i = 0; i < 3; i++) { - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPH::initial_integrate_respa(int vflag, int ilevel, int iloop) -{ - // set timesteps by level - - double dtfm; - dtv = step_respa[ilevel]; - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // atom quantities - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - // outermost level - update omega_dot, apply to v, remap box - // all other levels - NVE update of v - // x,v updates only performed for atoms in NPH group - - if (ilevel == nlevels_respa-1) { - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = nkt / volume * nktv2p; - - for (int i = 0; i < 3; i++) { - p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - omega[i] += dtv*omega_dot[i]; - factor[i] = exp(-dthalf*omega_dot[i]); - dilation[i] = exp(dthalf*omega_dot[i]); - } - - // v update only for atoms in NPH group - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - for (int 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]; - } - } - } - - // remap simulation box by 1/2 step - - remap(); - remap2flag = 1; - - } else { - - // v update only for atoms in NPH group - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else { - 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]; - } - } - } - } - - // innermost level - also update x for atoms in group - // if remap2flag: - // this is 1st call at innermost level from rRESPA after 1st half remap - // perform 2nd half of box remap - // redo KSpace coeffs since volume has changed - - if (ilevel == 0) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - if (remap2flag) { - remap(); - if (kspace_flag) force->kspace->setup(); - remap2flag = 0; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPH::final_integrate_respa(int ilevel, int iloop) -{ - double dtfm; - - // set timesteps by level - - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // outermost level - update omega_dot, apply to v via final_integrate() - // all other levels - NVE update of v - // v update only performed for atoms in NPH group - - if (ilevel == nlevels_respa-1) final_integrate(); - else { - - // v update only for atoms in NPH group - - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else { - 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]; - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPH::couple() -{ - double *tensor = pressure->vector; - - if (press_couple == XYZ) - p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == XY) { - double ave = 0.5 * (tensor[0] + tensor[1]); - p_current[0] = p_current[1] = ave; - p_current[2] = tensor[2]; - } else if (press_couple == YZ) { - double ave = 0.5 * (tensor[1] + tensor[2]); - p_current[1] = p_current[2] = ave; - p_current[0] = tensor[0]; - } else if (press_couple == XZ) { - double ave = 0.5 * (tensor[0] + tensor[2]); - p_current[0] = p_current[2] = ave; - p_current[1] = tensor[1]; - } else if (press_couple == ANISO) { - p_current[0] = tensor[0]; - p_current[1] = tensor[1]; - p_current[2] = tensor[2]; - } -} - -/* ---------------------------------------------------------------------- - change box size - remap all atoms or fix group atoms depending on allremap flag - if rigid bodies exist, scale rigid body centers-of-mass -------------------------------------------------------------------------- */ - -void FixNPH::remap() -{ - int i; - double oldlo,oldhi,ctr; - - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - // convert pertinent atoms and rigid bodies to lamda coords - - if (allremap) domain->x2lamda(nlocal); - else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->x2lamda(x[i],x[i]); - } - - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); - - // reset global and local box to new size/shape - - for (i = 0; i < 3; i++) { - if (p_flag[i]) { - oldlo = domain->boxlo[i]; - oldhi = domain->boxhi[i]; - ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr; - domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr; - } - } - - domain->set_global_box(); - domain->set_local_box(); - - // convert pertinent atoms and rigid bodies back to box coords - - if (allremap) domain->lamda2x(nlocal); - else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->lamda2x(x[i],x[i]); - } - - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixNPH::write_restart(FILE *fp) -{ - int n = 0; - double list[6]; - list[n++] = omega[0]; - list[n++] = omega[1]; - list[n++] = omega[2]; - list[n++] = omega_dot[0]; - list[n++] = omega_dot[1]; - list[n++] = omega_dot[2]; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixNPH::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - omega[0] = list[n++]; - omega[1] = list[n++]; - omega[2] = list[n++]; - omega_dot[0] = list[n++]; - omega_dot[1] = list[n++]; - omega_dot[2] = list[n++]; -} - -/* ---------------------------------------------------------------------- */ - -int FixNPH::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all("Illegal fix_modify command"); - if (tflag) { - modify->delete_compute(id_temp); - tflag = 0; - } - 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 temperature ID"); - temperature = modify->compute[icompute]; - - if (temperature->tempflag == 0) - error->all("Fix_modify temperature ID does not compute temperature"); - if (temperature->igroup != 0 && comm->me == 0) - error->warning("Temperature for NPH is not for group all"); - - // reset id_temp of pressure to new temperature ID - - icompute = modify->find_compute(id_press); - if (icompute < 0) error->all("Pressure ID for fix npt does not exist"); - modify->compute[icompute]->reset_extra_compute_fix(id_temp); - - return 2; - - } else if (strcmp(arg[0],"press") == 0) { - if (narg < 2) error->all("Illegal fix_modify command"); - if (pflag) { - modify->delete_compute(id_press); - pflag = 0; - } - delete [] id_press; - int n = strlen(arg[1]) + 1; - id_press = new char[n]; - strcpy(id_press,arg[1]); - - int icompute = modify->find_compute(id_press); - if (icompute < 0) error->all("Could not find fix_modify pressure ID"); - pressure = modify->compute[icompute]; - - if (pressure->pressflag == 0) - error->all("Fix_modify pressure ID does not compute pressure"); - return 2; - } - return 0; -} - -/* ---------------------------------------------------------------------- */ - -double FixNPH::compute_scalar() -{ - double volume; - if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; - else volume = domain->xprd * domain->yprd; - - int pdim = p_flag[0] + p_flag[1] + p_flag[2]; - - double energy = 0.0; - for (int i = 0; i < 3; i++) - if (p_freq[i] > 0.0) - energy += 0.5*nkt*omega_dot[i]*omega_dot[i] / - (p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p); - - return energy; } diff --git a/src/fix_nph.h b/src/fix_nph.h index 1911e1d739..8cdbd8834f 100644 --- a/src/fix_nph.h +++ b/src/fix_nph.h @@ -20,54 +20,14 @@ FixStyle(nph,FixNPH) #ifndef LMP_FIX_NPH_H #define LMP_FIX_NPH_H -#include "fix.h" +#include "fix_nh.h" namespace LAMMPS_NS { -class FixNPH : public Fix { +class FixNPH : public FixNH { public: FixNPH(class LAMMPS *, int, char **); - ~FixNPH(); - int setmask(); - void init(); - void setup(int); - void initial_integrate(int); - void final_integrate(); - void initial_integrate_respa(int, int, int); - void final_integrate_respa(int, int); - double compute_scalar(); - void write_restart(FILE *); - void restart(char *); - int modify_param(int, char **); - - private: - int dimension; - double dtv,dtf,dthalf; - double boltz,nktv2p; - double vol0,nkt; - - int press_couple,allremap; - int p_flag[3]; // 1 if control P on this dim, 0 if not - double p_start[3],p_stop[3]; - double p_freq[3],p_target[3]; - double omega[3],omega_dot[3]; - double p_current[3],dilation[3]; - double drag,drag_factor; - double factor[3]; - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes - - int nlevels_respa; - double *step_respa; - int remap2flag; // flag for performing 2nd half remap() - - char *id_temp,*id_press; - class Compute *temperature,*pressure; - int tflag,pflag; - - void couple(); - void remap(); + ~FixNPH() {} }; } diff --git a/src/fix_nph_sphere.cpp b/src/fix_nph_sphere.cpp new file mode 100644 index 0000000000..134cf06c82 --- /dev/null +++ b/src/fix_nph_sphere.cpp @@ -0,0 +1,68 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "fix_nph_sphere.h" +#include "group.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixNPHSphere::FixNPHSphere(LAMMPS *lmp, int narg, char **arg) : + FixNHSphere(lmp, narg, arg) +{ + if (tstat_flag) + error->all("Temperature control can not be used with fix nph/sphere"); + if (!pstat_flag) + error->all("Pressure control must be used with fix nph/sphere"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/sphere"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pflag = 1; +} diff --git a/src/fix_nph_sphere.h b/src/fix_nph_sphere.h new file mode 100644 index 0000000000..76e4736731 --- /dev/null +++ b/src/fix_nph_sphere.h @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(nph/sphere,FixNPHSphere) + +#else + +#ifndef LMP_FIX_NPH_SPHERE_H +#define LMP_FIX_NPH_SPHERE_H + +#include "fix_nh_sphere.h" + +namespace LAMMPS_NS { + +class FixNPHSphere : public FixNHSphere { + public: + FixNPHSphere(class LAMMPS *, int, char **); + ~FixNPHSphere() {} +}; + +} + +#endif +#endif diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index bbd534ee8e..3963e69ebb 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -11,193 +11,38 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing author: Mark Stevens (SNL) -------------------------------------------------------------------------- */ - #include "string.h" -#include "stdlib.h" -#include "math.h" #include "fix_npt.h" -#include "atom.h" -#include "force.h" -#include "comm.h" #include "modify.h" -#include "fix_deform.h" -#include "compute.h" -#include "kspace.h" -#include "update.h" -#include "respa.h" -#include "domain.h" #include "error.h" using namespace LAMMPS_NS; -#define MIN(A,B) ((A) < (B)) ? (A) : (B) -#define MAX(A,B) ((A) > (B)) ? (A) : (B) - -enum{NOBIAS,BIAS}; -enum{XYZ,XY,YZ,XZ,ANISO}; - /* ---------------------------------------------------------------------- */ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + FixNH(lmp, narg, arg) { - if (narg < 7) error->all("Illegal fix npt command"); - - restart_global = 1; - box_change = 1; - time_integrate = 1; - scalar_flag = 1; - global_freq = 1; - extscalar = 1; - - t_start = atof(arg[3]); - t_stop = atof(arg[4]); - double t_period = atof(arg[5]); - - if (t_start < 0.0 || t_stop <= 0.0) - error->all("Target T for fix npt cannot be 0.0"); - - double p_period[3]; - if (strcmp(arg[6],"xyz") == 0) { - if (narg < 10) error->all("Illegal fix npt command"); - press_couple = XYZ; - p_start[0] = p_start[1] = p_start[2] = atof(arg[7]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[9]); - p_flag[0] = p_flag[1] = p_flag[2] = 1; - if (domain->dimension == 2) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } - - } else { - if (strcmp(arg[6],"xy") == 0) press_couple = XY; - else if (strcmp(arg[6],"yz") == 0) press_couple = YZ; - else if (strcmp(arg[6],"xz") == 0) press_couple = XZ; - else if (strcmp(arg[6],"aniso") == 0) press_couple = ANISO; - else error->all("Illegal fix npt command"); - - if (narg < 14) error->all("Illegal fix npt command"); - - if (domain->dimension == 2 && - (press_couple == XY || press_couple == YZ || press_couple == XZ)) - error->all("Invalid fix npt command for a 2d simulation"); - - if (strcmp(arg[7],"NULL") == 0) { - p_start[0] = p_stop[0] = p_period[0] = 0.0; - p_flag[0] = 0; - } else { - p_start[0] = atof(arg[7]); - p_stop[0] = atof(arg[8]); - p_flag[0] = 1; - } - if (strcmp(arg[9],"NULL") == 0) { - p_start[1] = p_stop[1] = p_period[1] = 0.0; - p_flag[1] = 0; - } else { - p_start[1] = atof(arg[9]); - p_stop[1] = atof(arg[10]); - p_flag[1] = 1; - } - if (strcmp(arg[11],"NULL") == 0) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } else { - if (domain->dimension == 2) - error->all("Invalid fix npt command for a 2d simulation"); - p_start[2] = atof(arg[11]); - p_stop[2] = atof(arg[12]); - p_flag[2] = 1; - } - - double period = atof(arg[13]); - if (p_flag[0]) p_period[0] = period; - if (p_flag[1]) p_period[1] = period; - if (p_flag[2]) p_period[2] = period; - } - - // process extra keywords - - drag = 0.0; - allremap = 1; - - int iarg; - if (press_couple == XYZ) iarg = 10; - else iarg = 14; - - while (iarg < narg) { - if (strcmp(arg[iarg],"drag") == 0) { - if (iarg+2 > narg) error->all("Illegal fix npt command"); - drag = atof(arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"dilate") == 0) { - if (iarg+2 > narg) error->all("Illegal fix npt command"); - if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; - else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0; - else error->all("Illegal fix npt command"); - iarg += 2; - } else error->all("Illegal fix npt command"); - } - - // error checks - - if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Invalid fix npt command pressure settings"); - if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) - error->all("Invalid fix npt command pressure settings"); - if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) - error->all("Invalid fix npt command pressure settings"); - - if (press_couple == XY && - (p_start[0] != p_start[1] || p_stop[0] != p_stop[1])) - error->all("Invalid fix npt command pressure settings"); - if (press_couple == YZ && - (p_start[1] != p_start[2] || p_stop[1] != p_stop[2])) - error->all("Invalid fix npt command pressure settings"); - if (press_couple == XZ && - (p_start[0] != p_start[2] || p_stop[0] != p_stop[2])) - error->all("Invalid fix npt command pressure settings"); - - if (p_flag[0] && domain->xperiodic == 0) - error->all("Cannot use fix npt on a non-periodic dimension"); - if (p_flag[1] && domain->yperiodic == 0) - error->all("Cannot use fix npt on a non-periodic dimension"); - if (p_flag[2] && domain->zperiodic == 0) - error->all("Cannot use fix npt on a non-periodic dimension"); - - // convert input periods to frequencies - - if (t_period <= 0.0 || (p_flag[0] && p_period[0] <= 0.0) || - (p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0)) - error->all("Fix npt periods must be > 0.0"); - - t_freq = 1.0 / t_period; - p_freq[0] = p_freq[1] = p_freq[2] = 0.0; - if (p_flag[0]) p_freq[0] = 1.0 / p_period[0]; - if (p_flag[1]) p_freq[1] = 1.0 / p_period[1]; - if (p_flag[2]) p_freq[2] = 1.0 / p_period[2]; + if (!tstat_flag) + error->all("Temperature control must be used with fix npt"); + if (!pstat_flag) + error->all("Pressure control must be used with fix npt"); // create a new compute temp style // id = fix-ID + temp // compute group = all since pressure is always global (group all) - // and thus its KE/temperature contribution should use group all + // and thus its KE/temperature contribution should use group all int n = strlen(id) + 6; id_temp = new char[n]; strcpy(id_temp,id); strcat(id_temp,"_temp"); - + char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = (char *) "all"; - 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"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); delete [] newarg; tflag = 1; @@ -210,7 +55,7 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : id_press = new char[n]; strcpy(id_press,id); strcat(id_press,"_press"); - + newarg = new char*[4]; newarg[0] = id_press; newarg[1] = (char *) "all"; @@ -219,766 +64,4 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) : modify->add_compute(4,newarg); delete [] newarg; pflag = 1; - - // Nose/Hoover temp and pressure init - - eta = eta_dot = 0.0; - omega[0] = omega[1] = omega[2] = 0.0; - omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0; - - nrigid = 0; - rfix = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixNPT::~FixNPT() -{ - delete [] rfix; - - // delete temperature and pressure if fix created them - - if (tflag) modify->delete_compute(id_temp); - if (pflag) modify->delete_compute(id_press); - delete [] id_temp; - delete [] id_press; -} - -/* ---------------------------------------------------------------------- */ - -int FixNPT::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= FINAL_INTEGRATE; - mask |= THERMO_ENERGY; - mask |= INITIAL_INTEGRATE_RESPA; - mask |= FINAL_INTEGRATE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixNPT::init() -{ - if (domain->triclinic) error->all("Cannot use fix npt with triclinic box"); - - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"deform") == 0) { - int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; - if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || - (p_flag[2] && dimflag[2])) - error->all("Cannot use fix npt and fix deform on same dimension"); - } - - // set temperature and pressure ptrs - - int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all("Temperature 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("Pressure ID for fix npt does not exist"); - pressure = modify->compute[icompute]; - - // set timesteps and frequencies - - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; - dthalf = 0.5 * update->dt; - - double freq = MAX(p_freq[0],p_freq[1]); - freq = MAX(freq,p_freq[2]); - drag_factor = 1.0 - (update->dt * freq * drag); - - boltz = force->boltz; - nktv2p = force->nktv2p; - dimension = domain->dimension; - if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd; - else vol0 = domain->xprd * domain->yprd; - - if (force->kspace) kspace_flag = 1; - else kspace_flag = 0; - - if (strcmp(update->integrate_style,"respa") == 0) { - nlevels_respa = ((Respa *) update->integrate)->nlevels; - step_respa = ((Respa *) update->integrate)->step; - } - - // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - - delete [] rfix; - nrigid = 0; - rfix = NULL; - - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; - } -} - -/* ---------------------------------------------------------------------- - compute T,P before integrator starts -------------------------------------------------------------------------- */ - -void FixNPT::setup(int vflag) -{ - t_target = t_start; // used by compute_scalar() - p_target[0] = p_start[0]; - p_target[1] = p_start[1]; - p_target[2] = p_start[2]; - - t_current = temperature->compute_scalar(); - if (press_couple == XYZ) { - double tmp = pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - - // trigger virial computation on next timestep - - pressure->addstep(update->ntimestep+1); -} - -/* ---------------------------------------------------------------------- - 1st half of Verlet update -------------------------------------------------------------------------- */ - -void FixNPT::initial_integrate(int vflag) -{ - int i; - double dtfm; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - // update eta_dot - - t_target = t_start + delta * (t_stop-t_start); - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = atom->natoms*boltz*t_target / volume * nktv2p; - - for (i = 0; i < 3; i++) { - p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - omega[i] += dtv*omega_dot[i]; - factor[i] = exp(-dthalf*(eta_dot+omega_dot[i])); - dilation[i] = exp(dthalf*omega_dot[i]); - } - - // update v and x of atoms in group - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (rmass) { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - } - } - } - - } else { - 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 { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - } - } - } - } - - // remap simulation box by 1/2 step - - remap(); - - // update x by full step for atoms in group - - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - - // remap simulation box by 1/2 step - // redo KSpace coeffs since volume has changed - - remap(); - if (kspace_flag) force->kspace->setup(); -} - -/* ---------------------------------------------------------------------- - 2nd half of Verlet update -------------------------------------------------------------------------- */ - -void FixNPT::final_integrate() -{ - int i; - double dtfm; - - // update v of atoms in group - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (rmass) { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - } - } - } - - } else { - 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 { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - } - } - } - } - - // compute new T,P - - t_current = temperature->compute_scalar(); - if (press_couple == XYZ) { - double tmp = pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - - // trigger virial computation on next timestep - - pressure->addstep(update->ntimestep+1); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = atom->natoms*boltz*t_target / volume * nktv2p; - - for (i = 0; i < 3; i++) { - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPT::initial_integrate_respa(int vflag, int ilevel, int iloop) -{ - // set timesteps by level - - double dtfm; - dtv = step_respa[ilevel]; - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // atom quantities - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // outermost level - update eta_dot and omega_dot, apply to v, remap box - // all other levels - NVE update of v - // x,v updates performed for atoms in group - - if (ilevel == nlevels_respa-1) { - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - // update eta_dot - - t_target = t_start + delta * (t_stop-t_start); - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = atom->natoms*boltz*t_target / volume * nktv2p; - - for (int i = 0; i < 3; i++) { - p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - omega[i] += dtv*omega_dot[i]; - factor[i] = exp(-dthalf*(eta_dot+omega_dot[i])); - dilation[i] = exp(dthalf*omega_dot[i]); - } - - // update v for atoms in group - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - for (int 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]; - } - } - } - - // remap simulation box by 1/2 step - - remap(); - remap2flag = 1; - - } else { - - // update v for atoms in group - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - - } else { - 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]; - } - } - } - } - - // innermost level - also update x for atoms in group - // if remap2flag: - // this is 1st call at innermost level from rRESPA after 1st half remap - // perform 2nd half of box remap - // redo KSpace coeffs since volume has changed - - if (ilevel == 0) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - if (remap2flag) { - remap(); - if (kspace_flag) force->kspace->setup(); - remap2flag = 0; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPT::final_integrate_respa(int ilevel, int iloop) -{ - double dtfm; - - // set timesteps by level - - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // outermost level - update eta_dot and omega_dot, - // apply to v via final_integrate() - // all other levels - NVE update of v - // update v for atoms in group - - if (ilevel == nlevels_respa-1) final_integrate(); - else { - - // update v for atoms in group - - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - - } else { - 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]; - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNPT::couple() -{ - double *tensor = pressure->vector; - - if (press_couple == XYZ) - p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == XY) { - double ave = 0.5 * (tensor[0] + tensor[1]); - p_current[0] = p_current[1] = ave; - p_current[2] = tensor[2]; - } else if (press_couple == YZ) { - double ave = 0.5 * (tensor[1] + tensor[2]); - p_current[1] = p_current[2] = ave; - p_current[0] = tensor[0]; - } else if (press_couple == XZ) { - double ave = 0.5 * (tensor[0] + tensor[2]); - p_current[0] = p_current[2] = ave; - p_current[1] = tensor[1]; - } else if (press_couple == ANISO) { - p_current[0] = tensor[0]; - p_current[1] = tensor[1]; - p_current[2] = tensor[2]; - } -} - -/* ---------------------------------------------------------------------- - change box size - remap all atoms or fix group atoms depending on allremap flag - if rigid bodies exist, scale rigid body centers-of-mass -------------------------------------------------------------------------- */ - -void FixNPT::remap() -{ - int i; - double oldlo,oldhi,ctr; - - double **x = atom->x; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - // convert pertinent atoms and rigid bodies to lamda coords - - if (allremap) domain->x2lamda(nlocal); - else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->x2lamda(x[i],x[i]); - } - - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); - - // reset global and local box to new size/shape - - for (i = 0; i < 3; i++) { - if (p_flag[i]) { - oldlo = domain->boxlo[i]; - oldhi = domain->boxhi[i]; - ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr; - domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr; - } - } - - domain->set_global_box(); - domain->set_local_box(); - - // convert pertinent atoms and rigid bodies back to box coords - - if (allremap) domain->lamda2x(nlocal); - else { - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - domain->lamda2x(x[i],x[i]); - } - - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixNPT::write_restart(FILE *fp) -{ - int n = 0; - double list[8]; - list[n++] = eta; - list[n++] = eta_dot; - list[n++] = omega[0]; - list[n++] = omega[1]; - list[n++] = omega[2]; - list[n++] = omega_dot[0]; - list[n++] = omega_dot[1]; - list[n++] = omega_dot[2]; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixNPT::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - eta = list[n++]; - eta_dot = list[n++]; - omega[0] = list[n++]; - omega[1] = list[n++]; - omega[2] = list[n++]; - omega_dot[0] = list[n++]; - omega_dot[1] = list[n++]; - omega_dot[2] = list[n++]; -} - -/* ---------------------------------------------------------------------- */ - -int FixNPT::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all("Illegal fix_modify command"); - if (tflag) { - modify->delete_compute(id_temp); - tflag = 0; - } - delete [] id_temp; - int n = strlen(arg[1]) + 1; - id_temp = new char[n]; - strcpy(id_temp,arg[1]); - - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) error->all("Could not find fix_modify temperature ID"); - temperature = modify->compute[icompute]; - - if (temperature->tempflag == 0) - error->all("Fix_modify temperature ID does not compute temperature"); - if (temperature->igroup != 0 && comm->me == 0) - error->warning("Temperature for fix modify is not for group all"); - - // reset id_temp of pressure to new temperature ID - - icompute = modify->find_compute(id_press); - if (icompute < 0) error->all("Pressure ID for fix modify does not exist"); - modify->compute[icompute]->reset_extra_compute_fix(id_temp); - - return 2; - - } else if (strcmp(arg[0],"press") == 0) { - if (narg < 2) error->all("Illegal fix_modify command"); - if (pflag) { - modify->delete_compute(id_press); - pflag = 0; - } - delete [] id_press; - int n = strlen(arg[1]) + 1; - id_press = new char[n]; - strcpy(id_press,arg[1]); - - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) error->all("Could not find fix_modify pressure ID"); - pressure = modify->compute[icompute]; - - if (pressure->pressflag == 0) - error->all("Fix_modify pressure ID does not compute pressure"); - return 2; - } - return 0; -} - -/* ---------------------------------------------------------------------- */ - -double FixNPT::compute_scalar() -{ - double ke = temperature->dof * boltz * t_target; - double keplus = atom->natoms * boltz * t_target; - double volume; - if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd; - else volume = domain->xprd * domain->yprd; - - int pdim = p_flag[0] + p_flag[1] + p_flag[2]; - - double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq)); - for (int i = 0; i < 3; i++) - if (p_freq[i] > 0.0) - energy += 0.5*keplus*omega_dot[i]*omega_dot[i] / - (p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p); - - return energy; -} - -/* ---------------------------------------------------------------------- */ - -void FixNPT::reset_dt() -{ - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; - dthalf = 0.5 * update->dt; - - double freq = MAX(p_freq[0],p_freq[1]); - freq = MAX(freq,p_freq[2]); - drag_factor = 1.0 - (update->dt * freq * drag); } diff --git a/src/fix_npt.h b/src/fix_npt.h index d46619d77f..50ded650b7 100644 --- a/src/fix_npt.h +++ b/src/fix_npt.h @@ -20,60 +20,14 @@ FixStyle(npt,FixNPT) #ifndef LMP_FIX_NPT_H #define LMP_FIX_NPT_H -#include "fix.h" +#include "fix_nh.h" namespace LAMMPS_NS { -class FixNPT : public Fix { +class FixNPT : public FixNH { public: FixNPT(class LAMMPS *, int, char **); - virtual ~FixNPT(); - int setmask(); - virtual void init(); - void setup(int); - virtual void initial_integrate(int); - virtual void final_integrate(); - void initial_integrate_respa(int, int, int); - void final_integrate_respa(int, int); - double compute_scalar(); - void write_restart(FILE *); - void restart(char *); - int modify_param(int, char **); - void reset_dt(); - - protected: - int dimension,which; - double dtv,dtf,dthalf; - double boltz,nktv2p; - double vol0; - - double t_start,t_stop; - double t_current,t_target; - double t_freq; - double f_eta,eta_dot,eta; - - int press_couple,allremap; - int p_flag[3]; // 1 if control P on this dim, 0 if not - double p_start[3],p_stop[3]; - double p_freq[3],p_target[3]; - double omega[3],omega_dot[3]; - double p_current[3],dilation[3]; - double drag,drag_factor; - double factor[3]; - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes - - int nlevels_respa; - double *step_respa; - int remap2flag; // flag for performing 2nd half remap() - - char *id_temp,*id_press; - class Compute *temperature,*pressure; - int tflag,pflag; - - void couple(); - void remap(); + ~FixNPT() {} }; } diff --git a/src/fix_npt_sphere.cpp b/src/fix_npt_sphere.cpp index 9bd08f2c40..2e250fad8d 100644 --- a/src/fix_npt_sphere.cpp +++ b/src/fix_npt_sphere.cpp @@ -12,497 +12,56 @@ ------------------------------------------------------------------------- */ #include "string.h" -#include "stdlib.h" -#include "math.h" #include "fix_npt_sphere.h" -#include "atom.h" -#include "atom_vec.h" -#include "force.h" -#include "compute.h" -#include "kspace.h" -#include "update.h" -#include "domain.h" +#include "modify.h" #include "error.h" using namespace LAMMPS_NS; -#define INERTIA 0.4 // moment of inertia for sphere - -enum{NOBIAS,BIAS}; - /* ---------------------------------------------------------------------- */ FixNPTSphere::FixNPTSphere(LAMMPS *lmp, int narg, char **arg) : - FixNPT(lmp, narg, arg) + FixNHSphere(lmp, narg, arg) { - // error checks + if (!tstat_flag) + error->all("Temperature control must be used with fix npt/sphere"); + if (!pstat_flag) + error->all("Pressure control must be used with fix npt/sphere"); - if (!atom->omega_flag || !atom->torque_flag) - error->all("Fix npt/sphere requires atom attributes omega, torque"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Fix npt/sphere requires atom attribute radius or shape"); -} - -/* ---------------------------------------------------------------------- */ - -void FixNPTSphere::init() -{ - int i,itype; - - // check that all particles are finite-size and spherical - // no point particles allowed - - if (atom->radius_flag) { - double *radius = atom->radius; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (radius[i] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - } - - } else { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Fix nvt/sphere requires spherical particle shapes"); - } - } - - FixNPT::init(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNPTSphere::initial_integrate(int vflag) -{ - int i,itype; - double dtfm,dtirotate; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - - // update eta_dot - - t_target = t_start + delta * (t_stop-t_start); - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = atom->natoms*boltz*t_target / volume * nktv2p; - - for (i = 0; i < 3; i++) { - p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]); - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - omega[i] += dtv*omega_dot[i]; - factor[i] = exp(-dthalf*(eta_dot+omega_dot[i])); - dilation[i] = exp(dthalf*omega_dot[i]); - } - factor_rotate = exp(-dthalf*eta_dot); - - // update v of atoms in group - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (rmass) { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - } - } - } - - } else { - 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 { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - } - } - } - } - - // remap simulation box by 1/2 step - - remap(); - - // update x by full step for atoms in group - - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - - // set timestep here since dt may have changed or come via rRESPA - - double dtfrotate = dtf / INERTIA; - - // update omega for all particles - // d_omega/dt = torque / inertia - // 4 cases depending on radius vs shape and rmass vs mass - - if (radius) { - if (rmass) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; - } - } - } else { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[type[i]]); - omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; - } - } - } - - } else { - if (rmass) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; - } - } - } else { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = omega[i][0]*factor_rotate + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor_rotate + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor_rotate + dtirotate*torque[i][2]; - } - } - } - } - - // remap simulation box by 1/2 step - // redo KSpace coeffs since volume has changed - - remap(); - if (kspace_flag) force->kspace->setup(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNPTSphere::final_integrate() -{ - int i,itype; - double dtfm,dtirotate; - - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // set timestep here since dt may have changed or come via rRESPA - - double dtfrotate = dtf / INERTIA; - - // update v,omega of atoms in group - // d_omega/dt = torque / inertia - // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - if (radius) { - if (rmass) { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } - - } else { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - 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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } else { - for (i = 0; i < nlocal; i++) { - double tmp = temperature->compute_scalar(); - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - 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(i,v[i]); - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } - } - - } else { - if (rmass) { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } - - } else { - if (which == NOBIAS) { - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - 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]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - 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(i,v[i]); - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * - factor_rotate; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * - factor_rotate; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * - factor_rotate; - } - } - } - } - } - - // compute new T,P - - t_current = temperature->compute_scalar(); - if (press_couple == 0) { - double tmp = pressure->compute_scalar(); - } else { - temperature->compute_vector(); - pressure->compute_vector(); - } - couple(); - - // trigger virial computation on next timestep - - pressure->addstep(update->ntimestep+1); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - - // update omega_dot - // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change - - double f_omega,volume; - if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd; - else volume = domain->xprd*domain->yprd; - double denskt = atom->natoms*boltz*t_target / volume * nktv2p; - - for (i = 0; i < 3; i++) { - f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt; - omega_dot[i] += f_omega*dthalf; - omega_dot[i] *= drag_factor; - } + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/sphere"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pflag = 1; } diff --git a/src/fix_npt_sphere.h b/src/fix_npt_sphere.h index 383bd970d0..1bbd59e73b 100644 --- a/src/fix_npt_sphere.h +++ b/src/fix_npt_sphere.h @@ -20,20 +20,14 @@ FixStyle(npt/sphere,FixNPTSphere) #ifndef LMP_FIX_NPT_SPHERE_H #define LMP_FIX_NPT_SPHERE_H -#include "fix_npt.h" +#include "fix_nh_sphere.h" namespace LAMMPS_NS { -class FixNPTSphere : public FixNPT { +class FixNPTSphere : public FixNHSphere { public: FixNPTSphere(class LAMMPS *, int, char **); ~FixNPTSphere() {} - void init(); - void initial_integrate(int); - void final_integrate(); - - private: - double factor_rotate; }; } diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 18a3fd12c0..7817113b05 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -11,693 +11,38 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -/* ---------------------------------------------------------------------- - Contributing authors: Mark Stevens (SNL) - Andy Ballard (U Maryland) for Nose/Hoover chains -------------------------------------------------------------------------- */ - #include "string.h" -#include "stdlib.h" -#include "math.h" #include "fix_nvt.h" -#include "atom.h" -#include "force.h" -#include "comm.h" #include "group.h" -#include "update.h" -#include "respa.h" #include "modify.h" -#include "compute.h" #include "error.h" using namespace LAMMPS_NS; -enum{NOBIAS,BIAS}; - /* ---------------------------------------------------------------------- */ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + FixNH(lmp, narg, arg) { - if (narg < 6) error->all("Illegal fix nvt command"); - - restart_global = 1; - time_integrate = 1; - scalar_flag = 1; - global_freq = 1; - extscalar = 1; - - t_start = atof(arg[3]); - t_stop = atof(arg[4]); - double t_period = atof(arg[5]); - - drag = 0.0; - chain = 1; - - int iarg = 6; - while (iarg < narg) { - if (strcmp(arg[iarg],"drag") == 0) { - if (iarg+2 > narg) error->all("Illegal fix nvt command"); - drag = atof(arg[iarg+1]); - if (drag < 0.0) error->all("Illegal fix nvt command"); - iarg += 2; - } else if (strcmp(arg[iarg],"chain") == 0) { - if (iarg+2 > narg) error->all("Illegal fix nvt command"); - if (strcmp(arg[iarg+1],"yes") == 0) chain = 1; - else if (strcmp(arg[iarg+1],"no") == 0) chain = 0; - else error->all("Illegal fix nvt command"); - iarg += 2; - } else error->all("Illegal fix nvt command"); - } - - // error checks - // convert input period to frequency - - if (t_start < 0.0 || t_stop <= 0.0) - error->all("Target T for fix nvt cannot be 0.0"); - if (t_period <= 0.0) error->all("Fix nvt period must be > 0.0"); - t_freq = 1.0 / t_period; + if (!tstat_flag) + error->all("Temperature control must be used with fix nvt"); + if (pstat_flag) + error->all("Pressure control can not be used with fix nvt"); // create a new compute temp style - // id = fix-ID + temp, compute group = fix group + // id = fix-ID + temp int n = strlen(id) + 6; id_temp = new char[n]; strcpy(id_temp,id); strcat(id_temp,"_temp"); - + char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = group->names[igroup]; - if (strcmp(style,"nvt") == 0) - newarg[2] = (char *) "temp"; - else if (strcmp(style,"nvt/asphere") == 0) - 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"; + newarg[2] = (char *) "temp"; + modify->add_compute(3,newarg); delete [] newarg; tflag = 1; - - // Nose/Hoover temp init - - eta = eta_dot = 0.0; - eta2 = eta2_dot = 0.0; -} - -/* ---------------------------------------------------------------------- */ - -FixNVT::~FixNVT() -{ - // delete temperature if fix created it - - if (tflag) modify->delete_compute(id_temp); - delete [] id_temp; -} - -/* ---------------------------------------------------------------------- */ - -int FixNVT::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= FINAL_INTEGRATE; - mask |= THERMO_ENERGY; - mask |= INITIAL_INTEGRATE_RESPA; - mask |= FINAL_INTEGRATE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::init() -{ - int icompute = modify->find_compute(id_temp); - if (icompute < 0) error->all("Temperature 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; - dtf = 0.5 * update->dt * force->ftm2v; - dthalf = 0.5 * update->dt; - dt4 = 0.25 * update->dt; - dt8 = 0.125 * update->dt; - - drag_factor = 1.0 - (update->dt * t_freq * drag); - - if (strcmp(update->integrate_style,"respa") == 0) { - nlevels_respa = ((Respa *) update->integrate)->nlevels; - step_respa = ((Respa *) update->integrate)->step; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::setup(int vflag) -{ - t_target = t_start; // used by compute_scalar() - t_current = temperature->compute_scalar(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::initial_integrate(int vflag) -{ - double dtfm; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - t_target = t_start + delta * (t_stop-t_start); - - // update eta, eta_dot, eta2, eta2_dot - - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (chain) { - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - factor2 = exp(-dt8*eta2_dot); - eta_dot *= factor2; - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dt4; - eta_dot *= drag_factor; - eta_dot *= factor2; - eta += dthalf*eta_dot; - eta2 += dthalf*eta2_dot; - - } else { - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - } - - factor = exp(-dthalf*eta_dot); - - // update v and x of atoms in group - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - } - - } else { - 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 { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - } - } - - if (chain) { - eta_dot *= factor2; - f_eta = t_freq*t_freq * (factor * factor * t_current/t_target - 1.0); - eta_dot += f_eta * dt4; - eta_dot *= factor2; - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::final_integrate() -{ - double dtfm; - - // update v of atoms in group - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (chain) factor = 1.0; - - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - } - } - } - - } else { - 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 { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - } - } - } - } - - // compute current T - - t_current = temperature->compute_scalar(); - - // update eta, eta_dot, eta2, eta2_dot - // chains require additional velocity update and recompute of current T - - if (chain) { - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - factor2 = exp(-dt8*eta2_dot); - eta_dot *= factor2; - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dt4; - eta_dot *= drag_factor; - eta_dot *= factor2; - eta += dthalf*eta_dot; - eta2 += dthalf*eta2_dot; - - factor = exp(-dthalf*eta_dot); - - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] *= factor; - v[i][1] *= factor; - v[i][2] *= factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - 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] *= factor; - v[i][1] *= factor; - v[i][2] *= factor; - temperature->restore_bias(i,v[i]); - } - } - } - - } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] *= factor; - v[i][1] *= factor; - v[i][2] *= factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - 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] *= factor; - v[i][1] *= factor; - v[i][2] *= factor; - temperature->restore_bias(i,v[i]); - } - } - } - } - - t_current = temperature->compute_scalar(); - - eta_dot *= factor2; - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta * dt4; - eta_dot *= factor2; - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - - } else { - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::initial_integrate_respa(int vflag, int ilevel, int iloop) -{ - // set timesteps by level - - double dtfm; - dtv = step_respa[ilevel]; - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // atom quantities - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // 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; - delta /= update->endstep - update->beginstep; - t_target = t_start + delta * (t_stop-t_start); - - if (chain) { - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - factor2 = exp(-dt8*eta2_dot); - eta_dot *= factor2; - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dt4; - eta_dot *= drag_factor; - eta_dot *= factor2; - eta += dthalf*eta_dot; - eta2 += dthalf*eta2_dot; - - } else { - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - } - - factor = exp(-dthalf*eta_dot); - - } else factor = 1.0; - - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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 { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - } - } - } - - } else { - 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]; - } - } - - } else { - double tmp = temperature->compute_scalar(); - 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(i,v[i]); - } - } - } - } - - if (ilevel == 0) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - } - - if (ilevel == nlevels_respa-1 && chain) { - eta_dot *= factor2; - f_eta = t_freq*t_freq * (factor * factor * t_current/t_target - 1.0); - eta_dot += f_eta * dt4; - eta_dot *= factor2; - eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::final_integrate_respa(int ilevel, int iloop) -{ - // set timesteps by level - - double dtfm; - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // outermost level - update eta_dot and apply to v via final_integrate() - // all other levels - NVE update of v - - if (ilevel == nlevels_respa-1) final_integrate(); - else { - double **v = atom->v; - double **f = atom->f; - double *rmass = atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - if (rmass) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - - } else { - 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]; - } - } - } - } -} - -/* ---------------------------------------------------------------------- - pack entire state of Fix into one write -------------------------------------------------------------------------- */ - -void FixNVT::write_restart(FILE *fp) -{ - int n = 0; - double list[4]; - list[n++] = eta; - list[n++] = eta_dot; - list[n++] = eta2; - list[n++] = eta2_dot; - - if (comm->me == 0) { - int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); - } -} - -/* ---------------------------------------------------------------------- - use state info from restart file to restart the Fix -------------------------------------------------------------------------- */ - -void FixNVT::restart(char *buf) -{ - int n = 0; - double *list = (double *) buf; - eta = list[n++]; - eta_dot = list[n++]; - eta2 = list[n++]; - eta2_dot = list[n++]; -} - -/* ---------------------------------------------------------------------- */ - -int FixNVT::modify_param(int narg, char **arg) -{ - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all("Illegal fix_modify command"); - if (tflag) { - modify->delete_compute(id_temp); - tflag = 0; - } - 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 temperature ID"); - temperature = modify->compute[icompute]; - - if (temperature->tempflag == 0) - error->all("Fix_modify temperature 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; -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::reset_target(double t_new) -{ - t_start = t_stop = t_new; -} - -/* ---------------------------------------------------------------------- */ - -void FixNVT::reset_dt() -{ - dtv = update->dt; - dtf = 0.5 * update->dt * force->ftm2v; - dthalf = 0.5 * update->dt; - dt4 = 0.25 * update->dt; - dt8 = 0.125 * update->dt; - - drag_factor = 1.0 - (update->dt * t_freq * drag); -} - -/* ---------------------------------------------------------------------- */ - -double FixNVT::compute_scalar() -{ - double ke = temperature->dof * force->boltz * t_target; - double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq)); - return energy; } diff --git a/src/fix_nvt.h b/src/fix_nvt.h index e7bfbfe784..0e43604092 100644 --- a/src/fix_nvt.h +++ b/src/fix_nvt.h @@ -20,42 +20,14 @@ FixStyle(nvt,FixNVT) #ifndef LMP_FIX_NVT_H #define LMP_FIX_NVT_H -#include "fix.h" +#include "fix_nh.h" namespace LAMMPS_NS { -class FixNVT : public Fix { +class FixNVT : public FixNH { public: FixNVT(class LAMMPS *, int, char **); - virtual ~FixNVT(); - int setmask(); - virtual void init(); - void setup(int); - virtual void initial_integrate(int); - virtual void final_integrate(); - virtual void initial_integrate_respa(int, int, int); - void final_integrate_respa(int, int); - double compute_scalar(); - void write_restart(FILE *); - void restart(char *); - int modify_param(int, char **); - void reset_target(double); - void reset_dt(); - - protected: - int which,chain; - double t_start,t_stop; - double t_current,t_target; - double t_freq,drag,drag_factor; - double f_eta,eta_dot,eta,eta2_dot,eta2,factor,factor2; - double dtv,dtf,dthalf,dt4,dt8; - - int nlevels_respa; - double *step_respa; - - char *id_temp; - class Compute *temperature; - int tflag; + ~FixNVT() {} }; } diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp index caecc8fe47..2a483ddb39 100644 --- a/src/fix_nvt_sllod.cpp +++ b/src/fix_nvt_sllod.cpp @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains @@ -16,16 +16,11 @@ ------------------------------------------------------------------------- */ #include "math.h" -#include "string.h" -#include "stdlib.h" -#include "fix_nvt_sllod.h" #include "math_extra.h" +#include "fix_nvt_sllod.h" #include "atom.h" -#include "force.h" #include "domain.h" #include "group.h" -#include "update.h" -#include "respa.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" @@ -35,17 +30,41 @@ using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp +enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ -FixNVTSlodd::FixNVTSlodd(LAMMPS *lmp, int narg, char **arg) : - FixNVT(lmp, narg, arg) {} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSlodd::init() +FixNVTSllod::FixNVTSllod(LAMMPS *lmp, int narg, char **arg) : + FixNH(lmp, narg, arg) { - FixNVT::init(); + if (!tstat_flag) + error->all("Temperature control must be used with fix nvt/sllod"); + if (pstat_flag) + error->all("Pressure control can not be used with fix nvt/sllod"); + + // create a new compute temp style + // id = fix-ID + temp + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "temp/deform"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixNVTSllod::init() +{ + FixNH::init(); if (!temperature->tempbias) error->all("Temperature for fix nvt/sllod does not have a bias"); @@ -67,25 +86,12 @@ void FixNVTSlodd::init() error->all("Using fix nvt/sllod with no fix deform defined"); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + perform half-step scaling of velocities +-----------------------------------------------------------------------*/ -void FixNVTSlodd::initial_integrate(int vflag) +void FixNVTSllod::nh_v_temp() { - double dtfm; - - double delta = update->ntimestep - update->firststep; - delta /= update->nsteps; - t_target = t_start + delta * (t_stop-t_start); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - factor = exp(-dthalf*eta_dot); - - // update v and x of atoms in group // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo // thermostat thermal velocity only // vdelu = SLLOD correction = Hrate*Hinv*vthermal @@ -95,11 +101,7 @@ void FixNVTSlodd::initial_integrate(int vflag) if (nondeformbias) double tmp = temperature->compute_scalar(); - double **x = atom->x; double **v = atom->v; - double **f = atom->f; - double *mass = atom->mass; - int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; @@ -113,144 +115,10 @@ void FixNVTSlodd::initial_integrate(int vflag) 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(i,v[i]); - - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSlodd::final_integrate() -{ - double dtfm; - - // update v of atoms in group - // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo - // thermostat thermal velocity only - // vdelu = SLLOD correction = Hrate*Hinv*vthermal - // for non temp/deform BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - if (nondeformbias) double tmp = temperature->compute_scalar(); - - double **v = atom->v; - double **f = atom->f; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - 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]]; - 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]; + v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0]; + v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1]; + v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2]; temperature->restore_bias(i,v[i]); } } - - // compute current T - - t_current = temperature->compute_scalar(); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; -} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag) -{ - if (flag) return; // only used by NPT,NPH - - // set timesteps by level - - double dtfm; - dtv = step_respa[ilevel]; - dtf = 0.5 * step_respa[ilevel] * force->ftm2v; - dthalf = 0.5 * step_respa[ilevel]; - - // atom quantities - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double *mass = atom->mass; - int *type = atom->type; - int *mask = atom->mask; - 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 - - if (ilevel == nlevels_respa-1) { - double delta = update->ntimestep - update->firststep; - delta /= update->nsteps; - t_target = t_start + delta * (t_stop-t_start); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - factor = exp(-dthalf*eta_dot); - } else factor = 1.0; - - // update v of atoms in group - - if (nondeformbias) double tmp = temperature->compute_scalar(); - - 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(i,v[i]); - } - } - - // innermost level - also update x of only atoms in group - - if (ilevel == 0) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - } - } - } } diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h index 4661bb92ba..d4d5dd7b3a 100644 --- a/src/fix_nvt_sllod.h +++ b/src/fix_nvt_sllod.h @@ -1,7 +1,7 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov + www.cs.sandia.gov/~sjplimp/lammps.html + Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains @@ -13,27 +13,27 @@ #ifdef FIX_CLASS -FixStyle(nvt/sllod,FixNVTSlodd) +FixStyle(nvt/sllod,FixNVTSllod) #else -#ifndef LMP_FIX_NVT_SLODD_H -#define LMP_FIX_NVT_SLODD_H +#ifndef LMP_FIX_NVT_SLLOD_H +#define LMP_FIX_NVT_SLLOD_H -#include "fix_nvt.h" +#include "fix_nh.h" namespace LAMMPS_NS { -class FixNVTSlodd : public FixNVT { +class FixNVTSllod : public FixNH { public: - FixNVTSlodd(class LAMMPS *, int, char **); + FixNVTSllod(class LAMMPS *, int, char **); + ~FixNVTSllod() {} void init(); - void initial_integrate(int); - void final_integrate(); - void initial_integrate_respa(int, int, int); private: int nondeformbias; + + void nh_v_temp(); }; } diff --git a/src/fix_nvt_sphere.cpp b/src/fix_nvt_sphere.cpp index be4f7d7740..9873917c9f 100644 --- a/src/fix_nvt_sphere.cpp +++ b/src/fix_nvt_sphere.cpp @@ -12,489 +12,37 @@ ------------------------------------------------------------------------- */ #include "string.h" -#include "stdlib.h" -#include "math.h" #include "fix_nvt_sphere.h" -#include "atom.h" -#include "atom_vec.h" -#include "force.h" -#include "comm.h" #include "group.h" -#include "update.h" #include "modify.h" -#include "compute.h" #include "error.h" using namespace LAMMPS_NS; -#define INERTIA 0.4 // moment of inertia for sphere - -enum{NOBIAS,BIAS}; - /* ---------------------------------------------------------------------- */ FixNVTSphere::FixNVTSphere(LAMMPS *lmp, int narg, char **arg) : - FixNVT(lmp, narg, arg) + FixNHSphere(lmp, narg, arg) { - // error checks + if (!tstat_flag) + error->all("Temperature control must be used with fix nvt/sphere"); + if (pstat_flag) + error->all("Pressure control can not be used with fix nvt/sphere"); - if (!atom->omega_flag || !atom->torque_flag) - error->all("Fix nvt/sphere requires atom attributes omega, torque"); - if (!atom->radius_flag && !atom->avec->shape_type) - error->all("Fix nvt/sphere requires atom attribute radius or shape"); -} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSphere::init() -{ - int i,itype; - - // check that all particles are finite-size and spherical - // no point particles allowed - - if (atom->radius_flag) { - double *radius = atom->radius; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (radius[i] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - } - - } else { - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - itype = type[i]; - if (shape[itype][0] == 0.0) - error->one("Fix nvt/sphere requires extended particles"); - if (shape[itype][0] != shape[itype][1] || - shape[itype][0] != shape[itype][2]) - error->one("Fix nvt/sphere requires spherical particle shapes"); - } - } - - FixNVT::init(); -} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSphere::initial_integrate(int vflag) -{ - int itype; - double dtfm,dtirotate; - - double delta = update->ntimestep - update->beginstep; - delta /= update->endstep - update->beginstep; - t_target = t_start + delta * (t_stop-t_start); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; - eta += dtv*eta_dot; - factor = exp(-dthalf*eta_dot); - - // update v and x of only atoms in group - - double **x = atom->x; - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // set timestep here since dt may have changed or come via rRESPA - - double dtfrotate = dtf / INERTIA; - - // update v,x,omega of atoms in group - // d_omega/dt = torque / inertia - // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - if (radius) { - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } - - } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - 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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - 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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } - } - - } else { - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } - - } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype]; - 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]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype]; - 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(i,v[i]); - x[i][0] += dtv * v[i][0]; - x[i][1] += dtv * v[i][1]; - x[i][2] += dtv * v[i][2]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = omega[i][0]*factor + dtirotate*torque[i][0]; - omega[i][1] = omega[i][1]*factor + dtirotate*torque[i][1]; - omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; - } - } - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void FixNVTSphere::final_integrate() -{ - int itype; - double dtfm,dtirotate; - - // update v,omega of atoms in group - // d_omega/dt = torque / inertia - // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias - // for BIAS: - // calculate temperature since some computes require temp - // computed on current nlocal atoms to remove bias - // OK to not test returned v = 0, since factor is multiplied by v - - double **v = atom->v; - double **f = atom->f; - double **omega = atom->omega; - double **torque = atom->torque; - double *radius = atom->radius; - double *rmass = atom->rmass; - double *mass = atom->mass; - double **shape = atom->shape; - int *type = atom->type; - int *mask = atom->mask; - int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // set timestep here since dt may have changed or come via rRESPA - - double dtfrotate = dtf / INERTIA; - - if (radius) { - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - - dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } - - } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype] * 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]; - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype] * 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(i,v[i]); - - dtirotate = dtfrotate / (radius[i]*radius[i]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } - } - - } else { - if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / rmass[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]; - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[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(i,v[i]); - - dtirotate = dtfrotate / (shape[itype][0]*shape[itype][0]*rmass[i]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } - - } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - dtfm = dtf / mass[itype] * 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]; - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } else { - double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - itype = type[i]; - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[itype] * 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(i,v[i]); - - dtirotate = dtfrotate / - (shape[itype][0]*shape[itype][0]*mass[itype]); - omega[i][0] = (omega[i][0] + dtirotate*torque[i][0]) * factor; - omega[i][1] = (omega[i][1] + dtirotate*torque[i][1]) * factor; - omega[i][2] = (omega[i][2] + dtirotate*torque[i][2]) * factor; - } - } - } - } - } - - // compute current T - - t_current = temperature->compute_scalar(); - - // update eta_dot - - f_eta = t_freq*t_freq * (t_current/t_target - 1.0); - eta_dot += f_eta*dthalf; - eta_dot *= drag_factor; + // create a new compute temp style + // id = fix-ID + temp + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "temp/sphere"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; } diff --git a/src/fix_nvt_sphere.h b/src/fix_nvt_sphere.h index 67c397641f..60dd0a52e0 100644 --- a/src/fix_nvt_sphere.h +++ b/src/fix_nvt_sphere.h @@ -20,17 +20,14 @@ FixStyle(nvt/sphere,FixNVTSphere) #ifndef LMP_FIX_NVT_SPHERE_H #define LMP_FIX_NVT_SPHERE_H -#include "fix_nvt.h" +#include "fix_nh_sphere.h" namespace LAMMPS_NS { -class FixNVTSphere : public FixNVT { +class FixNVTSphere : public FixNHSphere { public: FixNVTSphere(class LAMMPS *, int, char **); ~FixNVTSphere() {} - void init(); - void initial_integrate(int); - void final_integrate(); }; } diff --git a/src/fix_press_berendsen.cpp b/src/fix_press_berendsen.cpp index 154d149628..3d912e3e30 100644 --- a/src/fix_press_berendsen.cpp +++ b/src/fix_press_berendsen.cpp @@ -33,7 +33,8 @@ using namespace LAMMPS_NS; #define MAX(A,B) ((A) > (B)) ? (A) : (B) enum{NOBIAS,BIAS}; -enum{XYZ,XY,YZ,XZ,ANISO}; +enum{NONE,XYZ,XY,YZ,XZ}; +enum{ISO,ANISO}; /* ---------------------------------------------------------------------- */ @@ -44,80 +45,89 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : box_change = 1; - // Berendsen barostat should be applied every step + // Berendsen barostat applied every step nevery = 1; - if (strcmp(arg[3],"xyz") == 0) { - if (narg < 7) error->all("Illegal fix press/berendsen command"); - - press_couple = XYZ; - p_start[0] = p_start[1] = p_start[2] = atof(arg[4]); - p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]); - p_period[0] = p_period[1] = p_period[2] = atof(arg[6]); - p_flag[0] = p_flag[1] = p_flag[2] = 1; - if (domain->dimension == 2) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } - - } else { - if (strcmp(arg[3],"xy") == 0) press_couple = XY; - else if (strcmp(arg[3],"yz") == 0) press_couple = YZ; - else if (strcmp(arg[3],"xz") == 0) press_couple = XZ; - else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO; - else error->all("Illegal fix press/berendsen command"); - - if (narg < 8) error->all("Illegal fix press/berendsen command"); - - if (domain->dimension == 2 && - (press_couple == XY || press_couple == YZ || press_couple == XZ)) - error->all("Invalid fix press/berendsen command for a 2d simulation"); - - if (strcmp(arg[4],"NULL") == 0) { - p_start[0] = p_stop[0] = p_period[0] = 0.0; - p_flag[0] = 0; - } else { - p_start[0] = atof(arg[4]); - p_stop[0] = atof(arg[5]); - p_flag[0] = 1; - } - if (strcmp(arg[6],"NULL") == 0) { - p_start[1] = p_stop[1] = p_period[1] = 0.0; - p_flag[1] = 0; - } else { - p_start[1] = atof(arg[6]); - p_stop[1] = atof(arg[7]); - p_flag[1] = 1; - } - if (strcmp(arg[8],"NULL") == 0) { - p_start[2] = p_stop[2] = p_period[2] = 0.0; - p_flag[2] = 0; - } else { - if (domain->dimension == 2) - error->all("Invalid fix press/berendsen command for a 2d simulation"); - p_start[2] = atof(arg[8]); - p_stop[2] = atof(arg[9]); - p_flag[2] = 1; - } - - double period = atof(arg[10]); - if (p_flag[0]) p_period[0] = period; - if (p_flag[1]) p_period[1] = period; - if (p_flag[2]) p_period[2] = period; - } - - // process extra keywords + // default values + pcouple = NONE; bulkmodulus = 10.0; allremap = 1; - int iarg; - if (press_couple == XYZ) iarg = 7; - else iarg = 11; + for (int i = 0; i < 3; i++) { + p_start[i] = p_stop[i] = p_period[i] = 0.0; + p_flag[i] = 0; + p_period[i] = 0.0; + } + + // process keywords + + dimension = domain->dimension; + + int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"modulus") == 0) { + if (strcmp(arg[iarg],"iso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix press/berendsen command"); + pcouple = XYZ; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + } else if (strcmp(arg[iarg],"aniso") == 0) { + if (iarg+4 > narg) error->all("Illegal fix press/berendsen command"); + pcouple = NONE; + p_start[0] = p_start[1] = p_start[2] = atof(arg[iarg+1]); + p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[iarg+2]); + p_period[0] = p_period[1] = p_period[2] = atof(arg[iarg+3]); + p_flag[0] = p_flag[1] = p_flag[2] = 1; + if (dimension == 2) { + p_start[2] = p_stop[2] = p_period[2] = 0.0; + p_flag[2] = 0; + } + iarg += 4; + + } else if (strcmp(arg[iarg],"x") == 0) { + if (iarg+4 > narg) error->all("Illegal fix press/berendsen command"); + p_start[0] = atof(arg[iarg+1]); + p_stop[0] = atof(arg[iarg+2]); + p_period[0] = atof(arg[iarg+3]); + p_flag[0] = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"y") == 0) { + if (iarg+4 > narg) error->all("Illegal fix press/berendsen command"); + p_start[1] = atof(arg[iarg+1]); + p_stop[1] = atof(arg[iarg+2]); + p_period[1] = atof(arg[iarg+3]); + p_flag[1] = 1; + iarg += 4; + } else if (strcmp(arg[iarg],"z") == 0) { + if (iarg+4 > narg) error->all("Illegal fix press/berendsen command"); + p_start[2] = atof(arg[iarg+1]); + p_stop[2] = atof(arg[iarg+2]); + p_period[2] = atof(arg[iarg+3]); + p_flag[2] = 1; + iarg += 4; + if (dimension == 2) + error->all("Invalid fix press/berendsen command for a 2d simulation"); + + } else if (strcmp(arg[iarg],"couple") == 0) { + if (iarg+2 > narg) error->all("Illegal fix press/berendsen command"); + if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ; + else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY; + else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ; + else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ; + else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE; + else error->all("Illegal fix press/berendsen command"); + iarg += 2; + + } else if (strcmp(arg[iarg],"modulus") == 0) { if (iarg+2 > narg) error->all("Illegal fix press/berendsen command"); bulkmodulus = atof(arg[iarg+1]); if (bulkmodulus <= 0.0) @@ -134,21 +144,18 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : // error checks - if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); - if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); - if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) - error->all("Invalid fix press/berendsen command pressure settings"); + if (dimension == 2 && p_flag[2]) + error->all("Invalid fix press/berendsen command for a 2d simulation"); + if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) + error->all("Invalid fix press/berendsen command for a 2d simulation"); - if (press_couple == XY && - (p_start[0] != p_start[1] || p_stop[0] != p_stop[1])) + if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0 || p_flag[2] == 0)) error->all("Invalid fix press/berendsen command pressure settings"); - if (press_couple == YZ && - (p_start[1] != p_start[2] || p_stop[1] != p_stop[2])) + if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) error->all("Invalid fix press/berendsen command pressure settings"); - if (press_couple == XZ && - (p_start[0] != p_start[2] || p_stop[0] != p_stop[2])) + if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) + error->all("Invalid fix press/berendsen command pressure settings"); + if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) error->all("Invalid fix press/berendsen command pressure settings"); if (p_flag[0] && domain->xperiodic == 0) @@ -158,12 +165,34 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : if (p_flag[2] && domain->zperiodic == 0) error->all("Cannot use fix press/berendsen on a non-periodic dimension"); - if (p_flag[0] && p_period[0] <= 0.0) - error->all("Fix press/berendsen period must be > 0.0"); - if (p_flag[1] && p_period[1] <= 0.0) - error->all("Fix press/berendsen period must be > 0.0"); - if (p_flag[2] && p_period[2] <= 0.0) - error->all("Fix press/berendsen period must be > 0.0"); + if (pcouple == XYZ && + (p_start[0] != p_start[1] || p_start[0] != p_start[2] || + p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[1] || p_period[0] != p_period[2])) + error->all("Invalid fix press/berendsen pressure settings"); + if (pcouple == XY && + (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || + p_period[0] != p_period[1])) + error->all("Invalid fix press/berendsen pressure settings"); + if (pcouple == YZ && + (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] || + p_period[1] != p_period[2])) + error->all("Invalid fix press/berendsen pressure settings"); + if (pcouple == XZ && + (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] || + p_period[0] != p_period[2])) + error->all("Invalid fix press/berendsen pressure settings"); + + if ((p_flag[0] && p_period[0] <= 0.0) || + (p_flag[1] && p_period[1] <= 0.0) || + (p_flag[2] && p_period[2] <= 0.0)) + error->all("Fix press/berendsen damping parameters must be > 0.0"); + + // pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof + // else pstyle = ANISO -> 3 dof + + if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO; + else pstyle = ANISO; // create a new compute temp style // id = fix-ID + temp @@ -235,13 +264,15 @@ void FixPressBerendsen::init() if (domain->triclinic) error->all("Cannot use fix press/berendsen with triclinic box"); + // insure no conflict with fix deform + for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { int *dimflag = ((FixDeform *) modify->fix[i])->dimflag; if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) || (p_flag[2] && dimflag[2])) error->all("Cannot use fix press/berendsen and " - "fix deform on same dimension"); + "fix deform on same component of stress tensor"); } // set temperature and pressure ptrs @@ -298,7 +329,7 @@ void FixPressBerendsen::end_of_step() { // compute new T,P - if (press_couple == XYZ) { + if (pstyle == ISO) { double tmp = temperature->compute_scalar(); tmp = pressure->compute_scalar(); } else { @@ -336,21 +367,24 @@ void FixPressBerendsen::couple() { double *tensor = pressure->vector; - if (press_couple == XYZ) + if (pstyle == ISO) p_current[0] = p_current[1] = p_current[2] = pressure->scalar; - else if (press_couple == XY) { + else if (pcouple == XYZ) { + double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]); + p_current[0] = p_current[1] = p_current[2] = ave; + } else if (pcouple == XY) { double ave = 0.5 * (tensor[0] + tensor[1]); p_current[0] = p_current[1] = ave; p_current[2] = tensor[2]; - } else if (press_couple == YZ) { + } else if (pcouple == YZ) { double ave = 0.5 * (tensor[1] + tensor[2]); p_current[1] = p_current[2] = ave; p_current[0] = tensor[0]; - } else if (press_couple == XZ) { + } else if (pcouple == XZ) { double ave = 0.5 * (tensor[0] + tensor[2]); p_current[0] = p_current[2] = ave; p_current[1] = tensor[1]; - } if (press_couple == ANISO) { + } else { p_current[0] = tensor[0]; p_current[1] = tensor[1]; p_current[2] = tensor[2]; diff --git a/src/fix_press_berendsen.h b/src/fix_press_berendsen.h index e30c04b713..cbe9a642a1 100644 --- a/src/fix_press_berendsen.h +++ b/src/fix_press_berendsen.h @@ -38,7 +38,7 @@ class FixPressBerendsen : public Fix { int dimension,which; double bulkmodulus; - int press_couple,allremap; + int pstyle,pcouple,allremap; int p_flag[3]; // 1 if control P on this dim, 0 if not double p_start[3],p_stop[3]; double p_period[3],p_target[3]; diff --git a/src/math_extra.h b/src/math_extra.h index e907ed78e1..8f6bb8ed60 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -44,6 +44,9 @@ namespace MathExtra { inline void transpose_times3(const double mat1[3][3], const double mat2[3][3], double ans[3][3]); + inline void times3_transpose(const double mat1[3][3], + const double mat2[3][3], + double ans[3][3]); inline void invert3(const double mat[3][3], double ans[3][3]); inline void times_column3(const double mat[3][3], const double*vec, double *ans); @@ -53,6 +56,7 @@ namespace MathExtra { double ans[3][3]); inline void row_times3(const double *v, const double m[3][3], double *ans); + inline void scalar_times3(const double f, double m[3][3]); inline void mldivide3(const double mat[3][3], const double *vec, double *ans, LAMMPS_NS::Error *error); inline void write3(const double mat[3][3]); @@ -209,6 +213,24 @@ void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3], ans[2][2] = m[0][2]*m2[0][2]+m[1][2]*m2[1][2]+m[2][2]*m2[2][2]; } +/* ---------------------------------------------------------------------- + multiply mat1 times transpose of mat2 +------------------------------------------------------------------------- */ + +void MathExtra::times3_transpose(const double m[3][3], const double m2[3][3], + double ans[3][3]) +{ + ans[0][0] = m[0][0]*m2[0][0]+m[0][1]*m2[0][1]+m[0][2]*m2[0][2]; + ans[0][1] = m[0][0]*m2[1][0]+m[0][1]*m2[1][1]+m[0][2]*m2[1][2]; + ans[0][2] = m[0][0]*m2[2][0]+m[0][1]*m2[2][1]+m[0][2]*m2[2][2]; + ans[1][0] = m[1][0]*m2[0][0]+m[1][1]*m2[0][1]+m[1][2]*m2[0][2]; + ans[1][1] = m[1][0]*m2[1][0]+m[1][1]*m2[1][1]+m[1][2]*m2[1][2]; + ans[1][2] = m[1][0]*m2[2][0]+m[1][1]*m2[2][1]+m[1][2]*m2[2][2]; + ans[2][0] = m[2][0]*m2[0][0]+m[2][1]*m2[0][1]+m[2][2]*m2[0][2]; + ans[2][1] = m[2][0]*m2[1][0]+m[2][1]*m2[1][1]+m[2][2]*m2[1][2]; + ans[2][2] = m[2][0]*m2[2][0]+m[2][1]*m2[2][1]+m[2][2]*m2[2][2]; +} + /* ---------------------------------------------------------------------- invert a matrix does NOT checks for singular or badly scaled matrix @@ -285,6 +307,23 @@ void MathExtra::row_times3(const double *v, const double m[3][3], ans[2] = v[0]*m[0][2]+v[1]*m[1][2]+m[2][2]*v[2]; } +/* ---------------------------------------------------------------------- + matrix times scalar, in place +------------------------------------------------------------------------- */ + +inline void MathExtra::scalar_times3(const double f, double m[3][3]) +{ + m[0][0] *= f; + m[0][1] *= f; + m[0][2] *= f; + m[1][0] *= f; + m[1][1] *= f; + m[1][2] *= f; + m[2][0] *= f; + m[2][1] *= f; + m[2][2] *= f; +} + /* ---------------------------------------------------------------------- solve Ax = b or M ans = v use gaussian elimination & partial pivoting on matrix diff --git a/src/min_linesearch.cpp b/src/min_linesearch.cpp index da0db9d21c..b94afa3337 100644 --- a/src/min_linesearch.cpp +++ b/src/min_linesearch.cpp @@ -45,8 +45,8 @@ using namespace LAMMPS_NS; #define ALPHA_MAX 1.0 #define ALPHA_REDUCE 0.5 #define BACKTRACK_SLOPE 0.4 -#define QUADRATIC_TOL 0.1 #define IDEAL_TOL 1.0e-8 +#define QUADRATIC_TOL 0.1 #define EPS_QUAD 1.0e-28 // same as in other min classes @@ -245,6 +245,15 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha, } if (nextra_global) modify->min_store(); + // Important diagnostic: test the gradient against energy + // double etmp; + // double alphatmp = alphamax*1.0e-4; + // etmp = alpha_step(alphatmp,1,nfunc); + // printf("alpha = %g dele = %g dele_force = %g err = %g\n", + // alphatmp,etmp-eoriginal,-alphatmp*fdothall, + // etmp-eoriginal+alphatmp*fdothall); + // alpha_step(0.0,1,nfunc); + // backtrack with alpha until energy decrease is sufficient while (1) { @@ -254,7 +263,13 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha, de_ideal = -BACKTRACK_SLOPE*alpha*fdothall; de = ecurrent - eoriginal; - if (de <= de_ideal) return 0; + if (de <= de_ideal) { + if (nextra_global) { + int itmp = modify->min_reset_ref(); + if (itmp) ecurrent = energy_force(1); + } + return 0; + } // reduce alpha @@ -317,6 +332,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0,fh0,ff0; double dot[2],dotall[2]; double *xatom,*x0atom,*fatom,*hatom; + double alphamax; // fdothall = projection of search dir along downhill gradient // if search direction is not downhill, exit with error @@ -336,18 +352,18 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, if (output->thermo->normflag) fdothall /= atom->natoms; if (fdothall <= 0.0) return DOWNHILL; - // set alpha so no dof is changed by more than max allowed amount + // set alphamax so no dof is changed by more than max allowed amount // for atom coords, max amount = dmax // for extra per-atom dof, max amount = extra_max[] // for extra global dof, max amount is set by fix - // also insure alpha <= ALPHA_MAX + // also insure alphamax <= ALPHA_MAX // else will have to backtrack from huge value when forces are tiny // if all search dir components are already 0.0, exit with error hme = 0.0; for (i = 0; i < n3; i++) hme = MAX(hme,fabs(h[i])); MPI_Allreduce(&hme,&hmaxall,1,MPI_DOUBLE,MPI_MAX,world); - alpha = MIN(ALPHA_MAX,dmax/hmaxall); + alphamax = MIN(ALPHA_MAX,dmax/hmaxall); if (nextra_atom) for (m = 0; m < nextra_atom; m++) { hme = 0.0; @@ -355,15 +371,16 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, n = extra_nlen[m]; for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i])); MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world); - alpha = MIN(alpha,extra_max[m]/hmax); + alphamax = MIN(alphamax,extra_max[m]/hmax); hmaxall = MAX(hmaxall,hmax); } if (nextra_global) { double alpha_extra = modify->max_alpha(hextra); - alpha = MIN(alpha,alpha_extra); + alphamax = MIN(alphamax,alpha_extra); for (i = 0; i < nextra_global; i++) hmaxall = MAX(hmaxall,fabs(hextra[i])); } + if (hmaxall == 0.0) return ZEROFORCE; // store box and values of all dof at start of linesearch @@ -382,10 +399,20 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, // backtrack with alpha until energy decrease is sufficient // or until get to small energy change, then perform quadratic projection + alpha = alphamax; fhprev = fdothall; engprev = eoriginal; alphaprev = 0.0; + // Important diagnostic: test the gradient against energy + // double etmp; + // double alphatmp = alphamax*1.0e-4; + // etmp = alpha_step(alphatmp,1,nfunc); + // printf("alpha = %g dele = %g dele_force = %g err = %g\n", + // alphatmp,etmp-eoriginal,-alphatmp*fdothall, + // etmp-eoriginal+alphatmp*fdothall); + // alpha_step(0.0,1,nfunc); + while (1) { ecurrent = alpha_step(alpha,1,nfunc); @@ -437,6 +464,10 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) { ecurrent = alpha_step(alpha0,1,nfunc); + if (nextra_global) { + int itmp = modify->min_reset_ref(); + if (itmp) ecurrent = energy_force(1); + } return 0; } @@ -444,7 +475,14 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, de_ideal = -BACKTRACK_SLOPE*alpha*fdothall; de = ecurrent - eoriginal; - if (de <= de_ideal) return 0; + + if (de <= de_ideal) { + if (nextra_global) { + int itmp = modify->min_reset_ref(); + if (itmp) ecurrent = energy_force(1); + } + return 0; + } // save previous state diff --git a/src/modify.cpp b/src/modify.cpp index 92357b9eaa..265547319d 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -513,6 +513,21 @@ int Modify::min_dof() return ndof; } +/* ---------------------------------------------------------------------- + reset reference state of fix, only for relevant fixes +------------------------------------------------------------------------- */ + +int Modify::min_reset_ref() +{ + int itmp,itmpall; + itmpall = 0; + for (int i = 0; i < n_min_energy; i++) { + itmp = fix[list_min_energy[i]]->min_reset_ref(); + if (itmp) itmpall = 1; + } + return itmpall; +} + /* ---------------------------------------------------------------------- add a new fix or replace one with same ID ------------------------------------------------------------------------- */ diff --git a/src/modify.h b/src/modify.h index 03ad635996..904d8d41a1 100644 --- a/src/modify.h +++ b/src/modify.h @@ -68,6 +68,7 @@ class Modify : protected Pointers { void min_clearstore(); void min_pushstore(); void min_popstore(); + int min_reset_ref(); double max_alpha(double *); int min_dof();