diff --git a/src/USER-UEF/fix_nh_uef.cpp b/src/USER-UEF/fix_nh_uef.cpp
index 1bc9386ad..fe342f1ef 100644
--- a/src/USER-UEF/fix_nh_uef.cpp
+++ b/src/USER-UEF/fix_nh_uef.cpp
@@ -1,861 +1,861 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
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
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: David Nicholson (MIT)
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_nh_uef.h"
#include "atom.h"
#include "force.h"
#include "group.h"
#include "comm.h"
#include "citeme.h"
#include "irregular.h"
#include "modify.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "domain.h"
#include "error.h"
#include "output.h"
#include "timer.h"
#include "neighbor.h"
#include "compute_pressure_uef.h"
#include "compute_temp_uef.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{ISO,ANISO,TRICLINIC};
/* ----------------------------------------------------------------------
Put all of the uef-only keywords at the back of arg and make narg smaller
so FixNH::FixNH() only sees the keywords it knows. Save the numer of
remaining keywords in rem.
---------------------------------------------------------------------- */
-char ** FixNHUef::arg_kludge(int &narg, char **arg, int &rem)
+char ** LAMMPS_NS::uef_arg_kludge(int &narg, char **arg, int &rem)
{
int iarg = 3;
bool flags[3]= {false,false,false};
rem=0;
char *tmp[3];
while (iarg < narg)
{
if (strcmp(arg[iarg],"erate" ) == 0 && !flags[0])
{
tmp[0] = arg[iarg];
tmp[1] = arg[iarg+1];
tmp[2] = arg[iarg+2];
for (int k=iarg+3; k<narg; k++)
arg[k-3] = arg[k];
arg[narg-1] = tmp[2];
arg[narg-2] = tmp[1];
arg[narg-3] = tmp[0];
rem += 3;
flags[0] = true;
}
else if (strcmp(arg[iarg],"strain" ) == 0 && !flags[1])
{
tmp[0] = arg[iarg];
tmp[1] = arg[iarg+1];
tmp[2] = arg[iarg+2];
for (int k=iarg+3; k<narg; k++)
arg[k-3] = arg[k];
arg[narg-1] = tmp[2];
arg[narg-2] = tmp[1];
arg[narg-3] = tmp[0];
rem += 3;
flags[1] = true;
}
else if(strcmp(arg[iarg],"ext") == 0 && !flags[2])
{
tmp[0] = arg[iarg];
tmp[1] = arg[iarg+1];
for (int k=iarg+2; k<narg; k++)
arg[k-2] = arg[k];
arg[narg-1] = tmp[1];
arg[narg-2] = tmp[0];
rem += 2;
flags[2] = true;
}
else
iarg++;
}
narg -= rem;
return arg;
}
//citation info
static const char cite_user_uef_package[] =
"USER-UEF package:\n\n"
"@Article{NicholsonRutledge16,\n"
"author = {David A. Nicholson and Gregory C. Rutledge},\n"
"title = {Molecular simulation of flow-enhanced nucleation in n-eicosane melts under steady shear and uniaxial extension},\n"
"journal = {The Journal of Chemical Physics},\n"
"volume = {145},\n"
"number = {24},\n"
"pages = {244903},\n"
"year = {2016}\n"
"}\n\n";
/* ----------------------------------------------------------------------
* Parse the remaing keywords, do some error checking, and initalize
* temp/pressure fixes
---------------------------------------------------------------------- */
FixNHUef::FixNHUef(LAMMPS *lmp, int narg, char **arg) :
- FixNH(lmp, narg, arg_kludge(narg,arg,rem))
+ FixNH(lmp, narg, uef_arg_kludge(narg,arg,rem))
{
if (lmp->citeme) lmp->citeme->add(cite_user_uef_package);
//initialization
erate[0] = erate[1] = 0;
//default values
strain[0]=strain[1]= 0;
ext_flags[0]=ext_flags[1]=ext_flags[2] = true;
// need to initialize these
omega_dot[0]=omega_dot[1]=omega_dot[2]=0;
// parse remaining input
bool erate_flag = false;
int iarg = narg;
narg += rem;
while (iarg <narg)
{
if (strcmp(arg[iarg],"erate")==0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix nvt/npt/uef command");
erate[0] = force->numeric(FLERR,arg[iarg+1]);
erate[1] = force->numeric(FLERR,arg[iarg+2]);
erate_flag = true;
iarg += 3;
}
else if (strcmp(arg[iarg],"strain")==0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix nvt/npt/uef command");
strain[0] = force->numeric(FLERR,arg[iarg+1]);
strain[1] = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
}
else if (strcmp(arg[iarg],"ext")==0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/uef command");
if (strcmp(arg[iarg+1],"x")==0)
ext_flags[1] = ext_flags[2] = false;
else if (strcmp(arg[iarg+1],"y")==0)
ext_flags[0] = ext_flags[2] = false;
else if (strcmp(arg[iarg+1],"z")==0)
ext_flags[0] = ext_flags[1] = false;
else if (strcmp(arg[iarg+1],"xy")==0)
ext_flags[2] = false;
else if (strcmp(arg[iarg+1],"xz")==0)
ext_flags[1] = false;
else if (strcmp(arg[iarg+1],"yz")==0)
ext_flags[0] = false;
else if (strcmp(arg[iarg+1],"xyz")!=0)
error->all(FLERR,"Illegal fix nvt/npt/uef command");
iarg += 2;
}
else
error->all(FLERR,"Illegal fix nvt/npt/uef command");
}
if (!erate_flag)
error->all(FLERR,"Keyword erate must be set for fix npt/npt/uef command");
if (mtchain_default_flag) mtchain=1;
if (!domain->triclinic)
error->all(FLERR,"Simulation box must be triclinic for fix/nvt/npt/uef");
//check for conditions that impose a deviatoric stress
if (pstyle == TRICLINIC)
error->all(FLERR,"Only normal stresses can be controlled with fix/nvt/npt/uef");
double erate_tmp[3];
erate_tmp[0]=erate[0];
erate_tmp[1]=erate[1];
erate_tmp[2]=-erate[0]-erate[1];
if (pstyle == ANISO)
{
if (!(ext_flags[0] & ext_flags[1] & ext_flags[2]))
error->all(FLERR,"The ext keyword may only be used with iso pressure control");
for (int k=0;k<3;k++)
for (int j=0;j<3;j++)
if (p_flag[k] && p_flag[j])
{
double tol = 1e-6;
if ( !nearly_equal(p_start[k],p_start[j],tol) || !nearly_equal(p_stop[k],p_stop[j],tol))
error->all(FLERR,"All controlled stresses must have the same value in fix/nvt/npt/uef");
if ( !nearly_equal(erate_tmp[k],erate_tmp[j],tol) || !nearly_equal(erate_tmp[k],erate_tmp[j],tol))
error->all(FLERR,"Dimensions with controlled stresses must have same strain rate in fix/nvt/npt/uef");
}
}
// conditions that produce a deviatoric stress have already
// been eliminated.
deviatoric_flag=0;
// need pre_exchange and irregular migration
pre_exchange_flag = 1;
irregular = new Irregular(lmp);
// flag that I change the box here (in case of nvt)
box_change_shape = 1;
// initialize the UEFBox class which computes the box at each step
uefbox = new UEF_utils::UEFBox();
uefbox->set_strain(strain[0],strain[1]);
// reset fixedpoint to the stagnation point. I don't allow fixedpoint
// to be set by the user.
fixedpoint[0] = domain->boxlo[0];
fixedpoint[1] = domain->boxlo[1];
fixedpoint[2] = domain->boxlo[2];
// Create temp and pressure computes for uef
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/uef";
modify->add_compute(3,newarg);
delete [] newarg;
tcomputeflag = 1;
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/uef";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
pcomputeflag = 1;
nevery = 1;
}
/* ----------------------------------------------------------------------
* Erase the UEFBox object and get rid of the pressure compute if the nvt
* version is being used. Everything else will be done in base destructor
* ---------------------------------------------------------------------- */
FixNHUef::~FixNHUef()
{
delete uefbox;
if (pcomputeflag && !pstat_flag)
{
modify->delete_compute(id_press);
delete [] id_press;
}
}
/* ----------------------------------------------------------------------
* Make the end_of_step() routine callable
* ---------------------------------------------------------------------- */
int FixNHUef::setmask()
{
int mask = FixNH::setmask();
mask |= END_OF_STEP;
return mask;
}
/* ----------------------------------------------------------------------
* Run FixNH::init() and do more error checking. Set the pressure
* pointer in the case that the nvt version is used
* ---------------------------------------------------------------------- */
void FixNHUef::init()
{
FixNH::init();
// find conflict with fix/deform or other box chaging fixes
for (int i=0; i < modify->nfix; i++)
{
if (strcmp(modify->fix[i]->id,id) != 0)
if (modify->fix[i]->box_change_shape != 0)
error->all(FLERR,"Can't use another fix which changes box shape with fix/nvt/npt/uef");
}
// this will make the pressure compute for nvt
if (!pstat_flag)
if (pcomputeflag)
{
int icomp = modify->find_compute(id_press);
if (icomp<0)
error->all(FLERR,"Pressure ID for fix/nvt/uef doesn't exist");
pressure = modify->compute[icomp];
}
if (strcmp(pressure->style,"pressure/uef") != 0)
error->all(FLERR,"Using fix nvt/npt/uef without a compute pressure/uef");
if (strcmp(temperature->style,"temp/uef") != 0)
error->all(FLERR,"Using fix nvt/npt/uef without a compute temp/uef");
}
/* ----------------------------------------------------------------------
* Run FixNH::setup() make sure the box is OK and set the rotation matrix
* for the first step
* ---------------------------------------------------------------------- */
void FixNHUef::setup(int j)
{
double box[3][3];
double vol = domain->xprd * domain->yprd * domain->zprd;
uefbox->get_box(box,vol);
double tol = 1e-4;
// ensure the box is ok for uef
bool isok = true;
isok &= nearly_equal(domain->h[0],box[0][0],tol);
isok &= nearly_equal(domain->h[1],box[1][1],tol);
isok &= nearly_equal(domain->h[2],box[2][2],tol);
isok &= nearly_equal(domain->xy,box[0][1],tol);
isok &= nearly_equal(domain->yz,box[1][2],tol);
isok &= nearly_equal(domain->xz,box[0][2],tol);
if (!isok)
error->all(FLERR,"Initial box is not close enough to the expected uef box");
uefbox->get_rot(rot);
((ComputeTempUef*) temperature)->yes_rot();
((ComputePressureUef*) pressure)->in_fix = true;
((ComputePressureUef*) pressure)->update_rot();
FixNH::setup(j);
}
/* ----------------------------------------------------------------------
* rotate -> initial integration step -> rotate back
* ---------------------------------------------------------------------- */
void FixNHUef::initial_integrate(int vflag)
{
inv_rotate_x(rot);
inv_rotate_v(rot);
inv_rotate_f(rot);
((ComputeTempUef*) temperature)->no_rot();
FixNH::initial_integrate(vflag);
rotate_x(rot);
rotate_v(rot);
rotate_f(rot);
((ComputeTempUef*) temperature)->yes_rot();
}
/* ----------------------------------------------------------------------
* rotate -> initial integration step -> rotate back (RESPA)
* ---------------------------------------------------------------------- */
void FixNHUef::initial_integrate_respa(int vflag, int ilevel, int iloop)
{
inv_rotate_x(rot);
inv_rotate_v(rot);
inv_rotate_f(rot);
((ComputeTempUef*) temperature)->no_rot();
FixNH::initial_integrate_respa(vflag,ilevel,iloop);
rotate_x(rot);
rotate_v(rot);
rotate_f(rot);
((ComputeTempUef*) temperature)->yes_rot();
}
/* ----------------------------------------------------------------------
* rotate -> final integration step -> rotate back
* ---------------------------------------------------------------------- */
void FixNHUef::final_integrate()
{
// update rot here since it must directly follow the virial calculation
((ComputePressureUef*) pressure)->update_rot();
inv_rotate_v(rot);
inv_rotate_f(rot);
((ComputeTempUef*) temperature)->no_rot();
FixNH::final_integrate();
rotate_v(rot);
rotate_f(rot);
((ComputeTempUef*) temperature)->yes_rot();
}
/* ----------------------------------------------------------------------
* at outer level: call this->final_integrate()
* at other levels: rotate -> 2nd verlet step -> rotate back
* ---------------------------------------------------------------------- */
void FixNHUef::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
{
inv_rotate_v(rot);
inv_rotate_f(rot);
nve_v();
rotate_v(rot);
rotate_f(rot);
}
}
/* ----------------------------------------------------------------------
SLLOD velocity update in time-reversible (i think) increments
v -> exp(-edot*dt/2)*v
v -> v +f/m*dt
v -> exp(-edot*dt/2)*v
-----------------------------------------------------------------------*/
void FixNHUef::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;
double ex = erate[0]*dtf/2;
double ey = erate[1]*dtf/2;
double ez = -ex-ey;
double e0 = exp(-ex);
double e1 = exp(-ey);
double e2 = exp(-ez);
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] *= e0;
v[i][1] *= e1;
v[i][2] *= e2;
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
v[i][0] *= e0;
v[i][1] *= e1;
v[i][2] *= e2;
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] *= e0;
v[i][1] *= e1;
v[i][2] *= e2;
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
v[i][0] *= e0;
v[i][1] *= e1;
v[i][2] *= e2;
}
}
}
}
/* ----------------------------------------------------------------------
Don't actually move atoms in remap(), just change the box
-----------------------------------------------------------------------*/
void FixNHUef::remap()
{
double vol = domain->xprd * domain->yprd * domain->zprd;
double domega = dto*(omega_dot[0]+omega_dot[1]+omega_dot[2])/3.;
// constant volume strain associated with barostat
// box scaling
double ex = dto*omega_dot[0]-domega;
double ey = dto*omega_dot[1]-domega;
uefbox->step_deform(ex,ey);
strain[0] += ex;
strain[1] += ey;
// volume change
vol = vol*exp(3*domega);
double box[3][3];
uefbox->get_box(box,vol);
domain->boxhi[0] = domain->boxlo[0]+box[0][0];
domain->boxhi[1] = domain->boxlo[1]+box[1][1];
domain->boxhi[2] = domain->boxlo[2]+box[2][2];
domain->xy = box[0][1];
domain->xz = box[0][2];
domain->yz = box[1][2];
domain->set_global_box();
domain->set_local_box();
uefbox->get_rot(rot);
}
/* ----------------------------------------------------------------------
SLLOD position update in time-reversible (i think) increments
x -> exp(edot*dt/2)*x
x -> x + v*dt
x -> exp(edot*dt/2)*x
-----------------------------------------------------------------------*/
void FixNHUef::nve_x()
{
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double ex = erate[0]*dtv;
strain[0] += ex;
double e0 = exp((ex+omega_dot[0]*dtv)/2);
double ey = erate[1]*dtv;
strain[1] += ey;
double e1 = exp((ey+omega_dot[1]*dtv)/2.);
double ez = -ex -ey;
double e2 = exp((ez+omega_dot[2]*dtv)/2.);
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] *= e0;
x[i][1] *= e1;
x[i][2] *= e2;
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
x[i][0] *= e0;
x[i][1] *= e1;
x[i][2] *= e2;
}
}
uefbox->step_deform(ex,ey);
double box[3][3];
double vol = domain->xprd * domain->yprd * domain->zprd;
uefbox->get_box(box,vol);
domain->boxhi[0] = domain->boxlo[0]+box[0][0];
domain->boxhi[1] = domain->boxlo[1]+box[1][1];
domain->boxhi[2] = domain->boxlo[2]+box[2][2];
domain->xy = box[0][1];
domain->xz = box[0][2];
domain->yz = box[1][2];
domain->set_global_box();
domain->set_local_box();
uefbox->get_rot(rot);
}
/* ----------------------------------------------------------------------
* Do the lattice reduction if necessary.
-----------------------------------------------------------------------*/
void FixNHUef::pre_exchange()
{
// only need to reset things if the lattice needs to be reduced
if (uefbox->reduce())
{
// go to lab frame
inv_rotate_x(rot);
inv_rotate_v(rot);
inv_rotate_f(rot);
// get & set the new box and rotation matrix
double vol = domain->xprd * domain->yprd * domain->zprd;
double box[3][3];
uefbox->get_box(box,vol);
domain->boxhi[0] = domain->boxlo[0]+box[0][0];
domain->boxhi[1] = domain->boxlo[1]+box[1][1];
domain->boxhi[2] = domain->boxlo[2]+box[2][2];
domain->xy = box[0][1];
domain->xz = box[0][2];
domain->yz = box[1][2];
domain->set_global_box();
domain->set_local_box();
uefbox->get_rot(rot);
// rotate to the new upper triangular frame
rotate_v(rot);
rotate_x(rot);
rotate_f(rot);
// put all atoms in the new box
double **x = atom->x;
imageint *image = atom->image;
int nlocal = atom->nlocal;
for (int i=0; i<nlocal; i++) domain->remap(x[i],image[i]);
// move atoms to the right processors
domain->x2lamda(atom->nlocal);
irregular->migrate_atoms();
domain->lamda2x(atom->nlocal);
}
}
/* ----------------------------------------------------------------------
* The following are routines to rotate between the lab and upper triangular
* (UT) frames. For most of the time the simulation is in the UT frame.
* To get to the lab frame, apply the inv_rotate_[..](rot) and to
* get back to the UT frame apply rotate_[..](rot).
*
* Note: the rotate_x() functions also apply a shift to/from the fixedpoint
* to make the integration a little simpler.
* ---------------------------------------------------------------------- */
void FixNHUef::rotate_x(double r[3][3])
{
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double xn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
xn[0]=r[0][0]*x[i][0]+r[0][1]*x[i][1]+r[0][2]*x[i][2];
xn[1]=r[1][0]*x[i][0]+r[1][1]*x[i][1]+r[1][2]*x[i][2];
xn[2]=r[2][0]*x[i][0]+r[2][1]*x[i][1]+r[2][2]*x[i][2];
x[i][0]=xn[0]+domain->boxlo[0];
x[i][1]=xn[1]+domain->boxlo[1];
x[i][2]=xn[2]+domain->boxlo[2];
}
}
}
void FixNHUef::inv_rotate_x(double r[3][3])
{
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double xn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
x[i][0] -= domain->boxlo[0];
x[i][1] -= domain->boxlo[1];
x[i][2] -= domain->boxlo[2];
xn[0]=r[0][0]*x[i][0]+r[1][0]*x[i][1]+r[2][0]*x[i][2];
xn[1]=r[0][1]*x[i][0]+r[1][1]*x[i][1]+r[2][1]*x[i][2];
xn[2]=r[0][2]*x[i][0]+r[1][2]*x[i][1]+r[2][2]*x[i][2];
x[i][0]=xn[0];
x[i][1]=xn[1];
x[i][2]=xn[2];
}
}
}
void FixNHUef::rotate_v(double r[3][3])
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double vn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
vn[0]=r[0][0]*v[i][0]+r[0][1]*v[i][1]+r[0][2]*v[i][2];
vn[1]=r[1][0]*v[i][0]+r[1][1]*v[i][1]+r[1][2]*v[i][2];
vn[2]=r[2][0]*v[i][0]+r[2][1]*v[i][1]+r[2][2]*v[i][2];
v[i][0]=vn[0]; v[i][1]=vn[1]; v[i][2]=vn[2];
}
}
}
void FixNHUef::inv_rotate_v(double r[3][3])
{
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double vn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
vn[0]=r[0][0]*v[i][0]+r[1][0]*v[i][1]+r[2][0]*v[i][2];
vn[1]=r[0][1]*v[i][0]+r[1][1]*v[i][1]+r[2][1]*v[i][2];
vn[2]=r[0][2]*v[i][0]+r[1][2]*v[i][1]+r[2][2]*v[i][2];
v[i][0]=vn[0]; v[i][1]=vn[1]; v[i][2]=vn[2];
}
}
}
void FixNHUef::rotate_f(double r[3][3])
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double fn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
fn[0]=r[0][0]*f[i][0]+r[0][1]*f[i][1]+r[0][2]*f[i][2];
fn[1]=r[1][0]*f[i][0]+r[1][1]*f[i][1]+r[1][2]*f[i][2];
fn[2]=r[2][0]*f[i][0]+r[2][1]*f[i][1]+r[2][2]*f[i][2];
f[i][0]=fn[0]; f[i][1]=fn[1]; f[i][2]=fn[2];
}
}
}
void FixNHUef::inv_rotate_f(double r[3][3])
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
double fn[3];
for (int i=0;i<nlocal;i++)
{
if (mask[i] & groupbit)
{
fn[0]=r[0][0]*f[i][0]+r[1][0]*f[i][1]+r[2][0]*f[i][2];
fn[1]=r[0][1]*f[i][0]+r[1][1]*f[i][1]+r[2][1]*f[i][2];
fn[2]=r[0][2]*f[i][0]+r[1][2]*f[i][1]+r[2][2]*f[i][2];
f[i][0]=fn[0]; f[i][1]=fn[1]; f[i][2]=fn[2];
}
}
}
/* ----------------------------------------------------------------------
* Increase the size of the restart list to add in the strains
* ---------------------------------------------------------------------- */
int FixNHUef::size_restart_global()
{
return FixNH::size_restart_global() +2;
}
/* ----------------------------------------------------------------------
* Pack the strains after packing the default FixNH values
* ---------------------------------------------------------------------- */
int FixNHUef::pack_restart_data(double *list)
{
int n = FixNH::pack_restart_data(list);
list[n++] = strain[0];
list[n++] = strain[1];
return n;
}
/* ----------------------------------------------------------------------
* read and set the strains after the default FixNH values
* ---------------------------------------------------------------------- */
void FixNHUef::restart(char *buf)
{
int n = size_restart_global();
FixNH::restart(buf);
double *list = (double *) buf;
strain[0] = list[n-2];
strain[1] = list[n-1];
uefbox->set_strain(strain[0],strain[1]);
}
/* ----------------------------------------------------------------------
* If the step writes a restart, reduce the box beforehand. This makes sure
* the unique box shape can be found once the restart is read and that
* all of the atoms lie within the box.
* This may only be necessary for RESPA runs, but I'm leaving it in anyway.
* ---------------------------------------------------------------------- */
void FixNHUef::end_of_step()
{
if (update->ntimestep==output->next_restart)
{
pre_exchange();
domain->x2lamda(atom->nlocal);
domain->pbc();
timer->stamp();
comm->exchange();
comm->borders();
domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
neighbor->build();
timer->stamp(Timer::NEIGH);
}
}
/* ----------------------------------------------------------------------
* reduce the simulation box after a run is complete. otherwise it won't
* be possible to resume from a write_restart since the initialization of
* the simulation box requires reduced simulation box
* ---------------------------------------------------------------------- */
void FixNHUef::post_run()
{
pre_exchange();
domain->x2lamda(atom->nlocal);
domain->pbc();
timer->stamp();
comm->exchange();
comm->borders();
domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
neighbor->build();
timer->stamp(Timer::NEIGH);
}
/* ----------------------------------------------------------------------
* public read for rotation matrix
* ---------------------------------------------------------------------- */
void FixNHUef::get_rot(double r[3][3])
{
r[0][0] = rot[0][0];
r[0][1] = rot[0][1];
r[0][2] = rot[0][2];
r[1][0] = rot[1][0];
r[1][1] = rot[1][1];
r[1][2] = rot[1][2];
r[2][0] = rot[2][0];
r[2][1] = rot[2][1];
r[2][2] = rot[2][2];
}
/* ----------------------------------------------------------------------
* public read for ext flags
* ---------------------------------------------------------------------- */
void FixNHUef::get_ext_flags(bool* e)
{
e[0] = ext_flags[0];
e[1] = ext_flags[1];
e[2] = ext_flags[2];
}
/* ----------------------------------------------------------------------
* public read for simulation box
* ---------------------------------------------------------------------- */
void FixNHUef::get_box(double b[3][3])
{
double box[3][3];
double vol = domain->xprd * domain->yprd * domain->zprd;
uefbox->get_box(box,vol);
b[0][0] = box[0][0];
b[0][1] = box[0][1];
b[0][2] = box[0][2];
b[1][0] = box[1][0];
b[1][1] = box[1][1];
b[1][2] = box[1][2];
b[2][0] = box[2][0];
b[2][1] = box[2][1];
b[2][2] = box[2][2];
}
/* ----------------------------------------------------------------------
* comparing floats
* it's imperfect, but should work provided no infinities
* ---------------------------------------------------------------------- */
bool FixNHUef::nearly_equal(double a, double b, double epsilon)
{
double absa = fabs(a);
double absb = fabs(b);
double diff = fabs(a-b);
if (a == b) return true;
else if ( (absa+absb) < epsilon)
return diff < epsilon*epsilon;
else
return diff/(absa+absb) < epsilon;
}
diff --git a/src/USER-UEF/fix_nh_uef.h b/src/USER-UEF/fix_nh_uef.h
index 46916005a..8b0f07a34 100644
--- a/src/USER-UEF/fix_nh_uef.h
+++ b/src/USER-UEF/fix_nh_uef.h
@@ -1,125 +1,126 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
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
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: David Nicholson (MIT)
------------------------------------------------------------------------- */
#ifndef LMP_FIX_NH_UEF_H
#define LMP_FIX_NH_UEF_H
#include "uef_utils.h"
#include "fix_nh.h"
namespace LAMMPS_NS {
+
+char **uef_arg_kludge(int&, char**, int&);
class FixNHUef : public FixNH {
public:
FixNHUef(class LAMMPS *, int, char **);
virtual ~FixNHUef();
virtual int setmask();
virtual void init();
virtual void setup(int);
virtual void pre_exchange();
virtual int pack_restart_data(double*);
virtual void restart(char *);
virtual void end_of_step();
virtual void initial_integrate(int);
virtual void final_integrate();
virtual void initial_integrate_respa(int, int, int);
virtual void final_integrate_respa(int, int);
virtual void post_run();
void get_rot(double[3][3]);
void get_ext_flags(bool*);
void get_box(double[3][3]);
protected:
virtual void remap();
virtual int size_restart_global();
virtual void nve_x();
virtual void nve_v();
- char **arg_kludge(int&, char**, int&);
void rotate_x(double [3][3]);
void inv_rotate_x(double[3][3]);
void rotate_v(double[3][3]);
void inv_rotate_v(double[3][3]);
void rotate_f(double[3][3]);
void inv_rotate_f(double[3][3]);
double strain[2],erate[2]; // strain/strain rate : [e_x, e_y]
// always assume traceless e_z = -e_x-e_y
int rem; //this is for the narg kluge
UEF_utils::UEFBox *uefbox; // interface for the special simulation box
double rot[3][3]; // rotation matrix
bool ext_flags[3]; // flags for external "free surfaces"
bool nearly_equal(double,double,double);
//bool rotate_output; // experimental feature. Too many issues for now
};
}
#endif
/* ERROR/WARNING messages:
This is a base class for FixNH so it will inherit most of its error/warning messages along with the following:
E: Illegal fix nvt/npt/uef command
Self-explanatory
E: Keyword erate must be set for fix nvt/npt/uef command
Self-explanatory.
E: Simulation box must be triclinic for fix/nvt/npt/uef
Self-explanatory.
E: Only normal stresses can be controlled with fix/nvt/npt/uef
The keywords xy xz and yz cannot be used for pressure control
E: The ext keyword may only be used with iso pressure control
Self-explanatory
E: All controlled stresses must have the same value in fix/nvt/npt/uef
Stress control is only possible when the stress specified for each dimension is the same
E: Dimensions with controlled stresses must have same strain rate in fix/nvt/npt/uef
Stress-controlled dimensions with the same strain rate must have the same target stress.
E: Can't use another fix which changes box shape with fix/nvt/npt/uef
The fix npt/nvt/uef command must have full control over the box shape. You cannot use a simultaneous fix deform command, for example.
E: Pressure ID for fix/nvt/uef doesn't exist
The compute pressure introduced via fix_modify does not exist
E: Using fix nvt/npt/uef without a compute pressure/uef
Self-explanatory.
E: Using fix nvt/npt/uef without a compute temp/uef
Self-explanatory.
E: Initial box is not close enough to the expected uef box
The initial box does not correspond to the shape required by the value of the strain keyword. If the default strain value of zero was used, the initial box is not cubic.
*/