diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index 12bc17852..61431d54d 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -1,394 +1,396 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "string.h" #include "compute_temp_asphere.h" #include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "group.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{ROTATE,ALL}; #define INERTIA 0.2 // moment of inertia prefactor for ellipsoid /* ---------------------------------------------------------------------- */ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 3) error->all(FLERR,"Illegal compute temp/asphere command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 0; id_bias = NULL; mode = ALL; int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"bias") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command"); tempbias = 1; int n = strlen(arg[iarg+1]) + 1; id_bias = new char[n]; strcpy(id_bias,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"dof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/asphere command"); if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE; else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL; else error->all(FLERR,"Illegal compute temp/asphere command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/asphere command"); } vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempAsphere::~ComputeTempAsphere() { delete [] id_bias; delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::init() { // error check avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); if (!avec) error->all(FLERR,"Compute temp/asphere requires atom style ellipsoid"); // check that all particles are finite-size, no point particles allowed int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (ellipsoid[i] < 0) error->one(FLERR,"Compute temp/asphere requires extended particles"); if (tempbias) { int i = modify->find_compute(id_bias); if (i < 0) error->all(FLERR,"Could not find compute ID for temperature bias"); tbias = modify->compute[i]; if (tbias->tempflag == 0) error->all(FLERR,"Bias compute does not calculate temperature"); if (tbias->tempbias == 0) error->all(FLERR,"Bias compute does not calculate a velocity bias"); if (tbias->igroup != igroup) error->all(FLERR,"Bias compute group does not match compute group"); tbias->init(); tbias->setup(); if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; else tempbias = 1; } } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::dof_compute() { if (fix_dof) adjust_dof_fix(); // 6 dof for 3d, 3 dof for 2d // which dof are included also depends on mode // assume full rotation of extended particles // user should correct this via compute_modify if needed double natoms = group->count(igroup); int nper; if (domain->dimension == 3) { if (mode == ALL) nper = 6; else nper = 3; } else { if (mode == ALL) nper = 3; else nper = 1; } dof = nper*natoms; // additional adjustments to dof if (tempbias == 1) { if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms; } else if (tempbias == 2) { int *mask = atom->mask; int nlocal = atom->nlocal; tbias->dof_remove_pre(); int count = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (tbias->dof_remove(i)) count++; int count_all; MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof -= nper*count_all; } dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempAsphere::compute_scalar() { invoked_scalar = update->ntimestep; if (tempbias) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } AtomVecEllipsoid::Bonus *bonus = avec->bonus; double **v = atom->v; double **angmom = atom->angmom; double *rmass = atom->rmass; int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; double *shape,*quat; double wbody[3],inertia[3]; double rot[3][3]; // sum translational and rotational energy for each particle // no point particles since divide by inertia double t = 0.0; if (mode == ALL) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; t += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; t += inertia[0]*wbody[0]*wbody[0] + inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic || tempbias == 2) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempAsphere::compute_vector() { int i; invoked_vector = update->ntimestep; if (tempbias) { if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } AtomVecEllipsoid::Bonus *bonus = avec->bonus; double **v = atom->v; double **angmom = atom->angmom; double *rmass = atom->rmass; int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; double *shape,*quat; double wbody[3],inertia[3],t[6]; double rot[3][3]; double massone; // sum translational and rotational energy for each particle // no point particles since divide by inertia for (i = 0; i < 6; i++) t[i] = 0.0; if (mode == ALL) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = rmass[i]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; // rotational kinetic energy t[0] += inertia[0]*wbody[0]*wbody[0]; t[1] += inertia[1]*wbody[1]*wbody[1]; t[2] += inertia[2]*wbody[2]*wbody[2]; t[3] += inertia[0]*wbody[0]*wbody[1]; t[4] += inertia[1]*wbody[0]*wbody[2]; t[5] += inertia[2]*wbody[1]*wbody[2]; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { // principal moments of inertia shape = bonus[ellipsoid[i]].shape; quat = bonus[ellipsoid[i]].quat; massone = rmass[i]; inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]); inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]); inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]); // wbody = angular velocity in body frame MathExtra::quat_to_mat(quat,rot); MathExtra::transpose_matvec(rot,angmom[i],wbody); wbody[0] /= inertia[0]; wbody[1] /= inertia[1]; wbody[2] /= inertia[2]; // rotational kinetic energy t[0] += inertia[0]*wbody[0]*wbody[0]; t[1] += inertia[1]*wbody[1]*wbody[1]; t[2] += inertia[2]*wbody[2]*wbody[2]; t[3] += inertia[0]*wbody[0]*wbody[1]; t[4] += inertia[1]*wbody[0]*wbody[2]; t[5] += inertia[2]*wbody[1]*wbody[2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempAsphere::remove_bias(int i, double *v) { if (tbias) tbias->remove_bias(i,v); } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempAsphere::restore_bias(int i, double *v) { if (tbias) tbias->restore_bias(i,v); } diff --git a/src/USER-CUDA/compute_temp_cuda.cpp b/src/USER-CUDA/compute_temp_cuda.cpp index b042de453..1db65294c 100644 --- a/src/USER-CUDA/compute_temp_cuda.cpp +++ b/src/USER-CUDA/compute_temp_cuda.cpp @@ -1,212 +1,215 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator Original Version: http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov See the README file in the top-level LAMMPS directory. ----------------------------------------------------------------------- USER-CUDA Package and associated modifications: https://sourceforge.net/projects/lammpscuda/ Christian Trott, christian.trott@tu-ilmenau.de Lars Winterfeld, lars.winterfeld@tu-ilmenau.de Theoretical Physics II, University of Technology Ilmenau, Germany See the README file in the USER-CUDA directory. This software is distributed under the GNU General Public License. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- 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 "mpi.h" #include #include #include #include "compute_temp_cuda.h" #include "compute_temp_cuda_cu.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" #include "error.h" #include "user_cuda.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempCuda::ComputeTempCuda(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { cuda = lmp->cuda; if(cuda == NULL) error->all(FLERR,"You cannot use a /cuda class, without activating 'cuda' acceleration. Provide '-c on' as command-line argument to LAMMPS.."); if (narg != 3) error->all(FLERR,"Illegal compute temp/cuda command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[6]; cu_t_vector = 0; cu_t_scalar = 0; cudable=true; } /* ---------------------------------------------------------------------- */ ComputeTempCuda::~ComputeTempCuda() { delete [] vector; delete cu_t_vector; delete cu_t_scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempCuda::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; + fix_dof = 0; for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempCuda::dof_compute() { double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempCuda::compute_scalar() { if(cuda->begin_setup) { if(not cu_t_vector) cu_t_vector = new cCudaData (t_vector,6); if(not cu_t_scalar) cu_t_scalar = new cCudaData (&t_scalar,1); invoked_scalar = update->ntimestep; Cuda_ComputeTempCuda_Scalar(&cuda->shared_data,groupbit,(ENERGY_CFLOAT*) cu_t_scalar->dev_data()); cu_t_scalar->download(); } else { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; } t_scalar=t; } MPI_Allreduce(&t_scalar,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; if(scalar>1e15) { cuda->cu_v->download(); cuda->cu_x->download(); cuda->cu_type->download(); double **v = atom->v; double **x = atom->x; printf("Out of v-range atoms: \n"); for(int i=0;inlocal;i++) if((v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2])>1e5) printf("%i %i // %lf %lf %lf // %lf %lf %lf\n",atom->tag[i],atom->type[i],x[i][0], x[i][1], x[i][2],v[i][0], v[i][1], v[i][2]); error->all(FLERR,"Temperature out of range. Simulations will be abortet.\n"); } return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempCuda::compute_vector() { int i; if(cuda->begin_setup) { if(not cu_t_vector) cu_t_vector = new cCudaData (t_vector,6); if(not cu_t_scalar) cu_t_scalar = new cCudaData (&t_scalar,1); invoked_vector = update->ntimestep; Cuda_ComputeTempCuda_Vector(&cuda->shared_data,groupbit,(ENERGY_CFLOAT*) cu_t_vector->dev_data()); cu_t_vector->download(); } else { invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; } for (i = 0; i < 6; i++) t_vector[i]=t[i]; } MPI_Allreduce(t_vector,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/USER-CUDA/compute_temp_partial_cuda.cpp b/src/USER-CUDA/compute_temp_partial_cuda.cpp index 6df69fa1b..3de4fbb35 100644 --- a/src/USER-CUDA/compute_temp_partial_cuda.cpp +++ b/src/USER-CUDA/compute_temp_partial_cuda.cpp @@ -1,357 +1,360 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator Original Version: http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov See the README file in the top-level LAMMPS directory. ----------------------------------------------------------------------- USER-CUDA Package and associated modifications: https://sourceforge.net/projects/lammpscuda/ Christian Trott, christian.trott@tu-ilmenau.de Lars Winterfeld, lars.winterfeld@tu-ilmenau.de Theoretical Physics II, University of Technology Ilmenau, Germany See the README file in the USER-CUDA directory. This software is distributed under the GNU General Public License. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- 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 "mpi.h" #include #include #include #include "compute_temp_partial_cuda.h" #include "compute_temp_partial_cuda_cu.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" #include "memory.h" #include "error.h" #include "user_cuda.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempPartialCuda::ComputeTempPartialCuda(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { cuda = lmp->cuda; if(cuda == NULL) error->all(FLERR,"You cannot use a /cuda class, without activating 'cuda' acceleration. Provide '-c on' as command-line argument to LAMMPS.."); if (narg != 6) error->all(FLERR,"Illegal compute temp/partial command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; xflag = force->inumeric(FLERR,arg[3]); yflag = force->inumeric(FLERR,arg[4]); zflag = force->inumeric(FLERR,arg[5]); if (zflag && domain->dimension == 2) error->all(FLERR,"Compute temp/partial cannot use vz for 2d systemx"); maxbias = 0; vbiasall = NULL; vector = new double[6]; cu_t_vector = 0; cu_t_scalar = 0; cu_vbiasall=NULL; cudable=true; } /* ---------------------------------------------------------------------- */ ComputeTempPartialCuda::~ComputeTempPartialCuda() { memory->destroy(vbiasall); delete [] vector; delete cu_t_vector; delete cu_t_scalar; delete cu_vbiasall; } /* ---------------------------------------------------------------------- */ void ComputeTempPartialCuda::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; + fix_dof = 0; for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempPartialCuda::dof_compute() { double natoms = group->count(igroup); int nper = xflag+yflag+zflag; dof = nper * natoms; dof -= (1.0*nper/domain->dimension)*fix_dof + extra_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ int ComputeTempPartialCuda::dof_remove(int i) { int nper = xflag+yflag+zflag; return (domain->dimension - nper); } /* ---------------------------------------------------------------------- */ double ComputeTempPartialCuda::compute_scalar() { if(cuda->begin_setup) { if(not cu_t_vector) cu_t_vector = new cCudaData (t_vector,6); if(not cu_t_scalar) cu_t_scalar = new cCudaData (&t_scalar,1); invoked_scalar = update->ntimestep; Cuda_ComputeTempPartialCuda_Scalar(&cuda->shared_data,groupbit,(ENERGY_CFLOAT*) cu_t_scalar->dev_data(),xflag,yflag,zflag); cu_t_scalar->download(); } else { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * rmass[i]; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * mass[type[i]]; } t_scalar=t; } MPI_Allreduce(&t_scalar,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; if(scalar>1e15) { cuda->cu_v->download(); cuda->cu_x->download(); cuda->cu_type->download(); double **v = atom->v; double **x = atom->x; printf("Out of v-range atoms: \n"); for(int i=0;inlocal;i++) if((v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2])>1e5) printf("%i %i // %lf %lf %lf // %lf %lf %lf\n",atom->tag[i],atom->type[i],x[i][0], x[i][1], x[i][2],v[i][0], v[i][1], v[i][2]); error->all(FLERR,"Temperature out of range. Simulations will be abortet.\n"); } return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempPartialCuda::compute_vector() { int i; if(cuda->begin_setup) { if(not cu_t_vector) cu_t_vector = new cCudaData (t_vector,6); if(not cu_t_scalar) cu_t_scalar = new cCudaData (&t_scalar,1); invoked_vector = update->ntimestep; Cuda_ComputeTempPartialCuda_Vector(&cuda->shared_data,groupbit,(ENERGY_CFLOAT*) cu_t_vector->dev_data(),xflag,yflag,zflag); cu_t_vector->download(); } else { invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * xflag*v[i][0]*v[i][0]; t[1] += massone * yflag*v[i][1]*v[i][1]; t[2] += massone * zflag*v[i][2]*v[i][2]; t[3] += massone * xflag*yflag*v[i][0]*v[i][1]; t[4] += massone * xflag*zflag*v[i][0]*v[i][2]; t[5] += massone * yflag*zflag*v[i][1]*v[i][2]; } for (i = 0; i < 6; i++) t_vector[i]=t[i]; } MPI_Allreduce(t_vector,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempPartialCuda::remove_bias(int i, double *v) { if (!xflag) { vbias[0] = v[0]; v[0] = 0.0; } if (!yflag) { vbias[1] = v[1]; v[1] = 0.0; } if (!zflag) { vbias[2] = v[2]; v[2] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempPartialCuda::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/partial:vbiasall"); delete cu_vbiasall; cu_vbiasall = new cCudaData ((double*)vbiasall, atom->nmax, 3); } if(cuda->begin_setup) { Cuda_ComputeTempPartialCuda_RemoveBiasAll(&cuda->shared_data,groupbit,xflag,yflag,zflag,cu_vbiasall->dev_data()); } else { if (!xflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][0] = v[i][0]; v[i][0] = 0.0; } } if (!yflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][1] = v[i][1]; v[i][1] = 0.0; } } if (!zflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][2] = v[i][2]; v[i][2] = 0.0; } } } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempPartialCuda::restore_bias(int i, double *v) { if (!xflag) v[0] += vbias[0]; if (!yflag) v[1] += vbias[1]; if (!zflag) v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempPartialCuda::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if(cuda->begin_setup) { Cuda_ComputeTempPartialCuda_RestoreBiasAll(&cuda->shared_data,groupbit,xflag,yflag,zflag,cu_vbiasall->dev_data()); } else { if (!xflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][0] += vbiasall[i][0]; } if (!yflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][1] += vbiasall[i][1]; } if (!zflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][2] += vbiasall[i][2]; } } } /* ---------------------------------------------------------------------- */ double ComputeTempPartialCuda::memory_usage() { double bytes = maxbias * sizeof(double); return bytes; } diff --git a/src/USER-EFF/compute_temp_deform_eff.cpp b/src/USER-EFF/compute_temp_deform_eff.cpp index 692372591..45639145c 100644 --- a/src/USER-EFF/compute_temp_deform_eff.cpp +++ b/src/USER-EFF/compute_temp_deform_eff.cpp @@ -1,318 +1,323 @@ /* ---------------------------------------------------------------------- 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: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_deform_eff.h" #include "domain.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "group.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp /* ---------------------------------------------------------------------- */ ComputeTempDeformEff::ComputeTempDeformEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/deform/eff command"); if (!atom->electron_flag) error->all(FLERR,"Compute temp/deform/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempDeformEff::~ComputeTempDeformEff() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::init() { - int i; - - fix_dof = 0; - for (i = 0; i < modify->nfix; i++) - fix_dof += modify->fix[i]->dof(igroup); - dof_compute(); - // check fix deform remap settings + int i; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP && comm->me == 0) error->warning(FLERR,"Using compute temp/deform/eff with inconsistent " "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) - error->warning(FLERR,"Using compute temp/deform/eff with no fix deform defined"); + error->warning(FLERR, + "Using compute temp/deform/eff with no fix deform defined"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempDeformEff::setup() +{ + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; + fix_dof = -1; + dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::dof_compute() { double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; // just include nuclear dof int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (fabs(spin[i]) == 1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempDeformEff::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo // vthermal = thermal velocity = v - vstream double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (mass) { t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2])* mass[type[i]]; if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempDeformEff::compute_vector() { double lamda[3],vstream[3],vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; if (fabs(spin[i])==1) { t[0] += mefactor * massone * ervel[i]*ervel[i]; t[1] += mefactor * massone * ervel[i]*ervel[i]; t[2] += mefactor * massone * ervel[i]*ervel[i]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeformEff::remove_bias(int i, double *v) { double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; domain->x2lamda(atom->x[i],lamda); vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity NOTE: only removes translational velocity bias from electrons ------------------------------------------------------------------------- */ void ComputeTempDeformEff::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/deform/eff:vbiasall"); } double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(atom->x[i],lamda); vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[i][0] -= vbiasall[i][0]; v[i][1] -= vbiasall[i][1]; v[i][2] -= vbiasall[i][2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeformEff::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeformEff::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempDeformEff::memory_usage() { double bytes = maxbias * sizeof(double); return bytes; } diff --git a/src/USER-EFF/compute_temp_eff.cpp b/src/USER-EFF/compute_temp_eff.cpp index 9b06f05a2..4dcf7720c 100644 --- a/src/USER-EFF/compute_temp_eff.cpp +++ b/src/USER-EFF/compute_temp_eff.cpp @@ -1,163 +1,165 @@ /* ---------------------------------------------------------------------- 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: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_eff.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempEff::ComputeTempEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Compute temp/eff requires atom style electron"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[size_vector]; } /* ---------------------------------------------------------------------- */ ComputeTempEff::~ComputeTempEff() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempEff::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempEff::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; int *spin = atom->spin; int *mask = atom->mask; int nlocal = atom->nlocal; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (fabs(spin[i])==1) one++; } int nelectrons; MPI_Allreduce(&one,&nelectrons,1,MPI_INT,MPI_SUM,world); // Assume 3/2 k T per nucleus dof -= domain->dimension * nelectrons; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempEff::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double t = 0.0; if (mass) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; if (fabs(spin[i])==1) t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; } } } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempEff::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; if (fabs(spin[i])==1) { t[0] += mefactor*massone*ervel[i]*ervel[i]; t[1] += mefactor*massone*ervel[i]*ervel[i]; t[2] += mefactor*massone*ervel[i]*ervel[i]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/USER-EFF/compute_temp_region_eff.cpp b/src/USER-EFF/compute_temp_region_eff.cpp index d0ee66e1b..e5a7a5c40 100644 --- a/src/USER-EFF/compute_temp_region_eff.cpp +++ b/src/USER-EFF/compute_temp_region_eff.cpp @@ -1,292 +1,294 @@ /* ---------------------------------------------------------------------- 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: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ #include "mpi.h" #include "math.h" #include "string.h" #include "compute_temp_region_eff.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "domain.h" #include "region.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (!atom->electron_flag) error->all(FLERR,"Compute temp/region/eff requires atom style electron"); if (narg != 4) error->all(FLERR,"Illegal compute temp/region/eff command"); iregion = domain->find_region(arg[3]); if (iregion == -1) error->all(FLERR,"Region ID for compute temp/region/eff does not exist"); int n = strlen(arg[3]) + 1; idregion = new char[n]; strcpy(idregion,arg[3]); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempRegionEff::~ComputeTempRegionEff() { delete [] idregion; memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::init() { // set index and check validity of region iregion = domain->find_region(idregion); if (iregion == -1) error->all(FLERR,"Region ID for compute temp/region/eff does not exist"); } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; dof = 0.0; } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::dof_remove_pre() { domain->regions[iregion]->prematch(); } /* ---------------------------------------------------------------------- */ int ComputeTempRegionEff::dof_remove(int i) { double *x = atom->x[i]; if (domain->regions[iregion]->match(x[0],x[1],x[2])) return 0; return 1; } /* ---------------------------------------------------------------------- */ double ComputeTempRegionEff::compute_scalar() { invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; region->prematch(); int count = 0; int ecount = 0; double t = 0.0; if (mass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { count++; t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; if (fabs(spin[i])==1) { t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; ecount++; } } } double tarray[2],tarray_all[2]; // Assume 3/2 k T per nucleus tarray[0] = count-ecount; tarray[1] = t; MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); dof = domain->dimension * tarray_all[0] - extra_dof; int one = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (fabs(spin[i])==1) one++; } if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); else scalar = 0.0; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempRegionEff::compute_vector() { int i; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *ervel = atom->ervel; double *mass = atom->mass; int *spin = atom->spin; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double mefactor = domain->dimension/4.0; Region *region = domain->regions[iregion]; region->prematch(); double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; if (fabs(spin[i])==1) { t[0] += mefactor * massone * ervel[i]*ervel[i]; t[1] += mefactor * massone * ervel[i]*ervel[i]; t[2] += mefactor * massone * ervel[i]*ervel[i]; } } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity NOTE: the following commands do not bias the radial electron velocities ------------------------------------------------------------------------- */ void ComputeTempRegionEff::remove_bias(int i, double *v) { double *x = atom->x[i]; if (domain->regions[iregion]->match(x[0],x[1],x[2])) vbias[0] = vbias[1] = vbias[2] = 0.0; else { vbias[0] = v[0]; vbias[1] = v[1]; vbias[2] = v[2]; v[0] = v[1] = v[2] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRegionEff::remove_bias_all() { double **x = atom->x; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/region:vbiasall"); } Region *region = domain->regions[iregion]; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (region->match(x[i][0],x[i][1],x[i][2])) vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0; else { vbiasall[i][0] = v[i][0]; vbiasall[i][1] = v[i][1]; vbiasall[i][2] = v[i][2]; v[i][0] = v[i][1] = v[i][2] = 0.0; } } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempRegionEff::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempRegionEff::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempRegionEff::memory_usage() { double bytes = maxbias * sizeof(double); return bytes; } diff --git a/src/USER-MISC/compute_temp_rotate.cpp b/src/USER-MISC/compute_temp_rotate.cpp index b34255be5..7c85fc1d1 100644 --- a/src/USER-MISC/compute_temp_rotate.cpp +++ b/src/USER-MISC/compute_temp_rotate.cpp @@ -1,277 +1,279 @@ /* ---------------------------------------------------------------------- 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: Laurent Joly (U Lyon, France), ljoly.ulyon@gmail.com ------------------------------------------------------------------------- */ #include "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_rotate.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "domain.h" #include "lattice.h" #include "error.h" #include "memory.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempRotate::ComputeTempRotate(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/rotate command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempRotate::~ComputeTempRotate() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempRotate::init() { masstotal = group->mass(igroup); } /* ---------------------------------------------------------------------- */ void ComputeTempRotate::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempRotate::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempRotate::compute_scalar() { double vthermal[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; double dx,dy,dz; double unwrap[3]; invoked_scalar = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vcm); group->xcm(igroup,masstotal,xcm); group->inertia(igroup,xcm,inertia); group->angmom(igroup,xcm,angmom); group->omega(angmom,inertia,omega); double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; imageint *image = atom->image; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/rotate:vbiasall"); } double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[2] = v[i][2] - vbiasall[i][2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempRotate::compute_vector() { double vthermal[3]; double vcm[3],xcm[3],inertia[3][3],angmom[3],omega[3]; double dx,dy,dz; double unwrap[3]; invoked_vector = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vcm); group->xcm(igroup,masstotal,xcm); group->inertia(igroup,xcm,inertia); group->angmom(igroup,xcm,angmom); group->omega(angmom,inertia,omega); double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; imageint *image = atom->image; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/rotate:vbiasall"); } double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->unmap(x[i],image[i],unwrap); dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; vbiasall[i][0] = vcm[0] + dz*omega[1]-dy*omega[2]; vbiasall[i][1] = vcm[1] + dx*omega[2]-dz*omega[0]; vbiasall[i][2] = vcm[2] + dy*omega[0]-dx*omega[1]; vthermal[0] = v[i][0] - vbiasall[i][0]; vthermal[1] = v[i][1] - vbiasall[i][1]; vthermal[2] = v[i][2] - vbiasall[i][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRotate::remove_bias(int i, double *v) { v[0] -= vbiasall[i][0]; v[1] -= vbiasall[i][1]; v[2] -= vbiasall[i][2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRotate::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] -= vbiasall[i][0]; v[i][1] -= vbiasall[i][1]; v[i][2] -= vbiasall[i][2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempRotate::restore_bias(int i, double *v) { v[0] += vbiasall[i][0]; v[1] += vbiasall[i][1]; v[2] += vbiasall[i][2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempRotate::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempRotate::memory_usage() { double bytes = maxbias * sizeof(double); return bytes; } diff --git a/src/compute.cpp b/src/compute.cpp index ca4251927..207823527 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -1,316 +1,317 @@ /* ---------------------------------------------------------------------- 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 "lmptype.h" #include "mpi.h" #include "stdlib.h" #include "string.h" #include "ctype.h" #include "compute.h" #include "atom.h" #include "domain.h" #include "comm.h" #include "group.h" #include "modify.h" #include "fix.h" #include "atom_masks.h" #include "memory.h" #include "error.h" #include "force.h" using namespace LAMMPS_NS; #define DELTA 4 #define BIG MAXTAGINT /* ---------------------------------------------------------------------- */ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) { if (narg < 3) error->all(FLERR,"Illegal compute command"); // compute ID, group, and style // ID must be all alphanumeric chars or underscores int n = strlen(arg[0]) + 1; id = new char[n]; strcpy(id,arg[0]); for (int i = 0; i < n-1; i++) if (!isalnum(id[i]) && id[i] != '_') error->all(FLERR, "Compute ID must be alphanumeric or underscore characters"); igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find compute group ID"); groupbit = group->bitmask[igroup]; n = strlen(arg[2]) + 1; style = new char[n]; strcpy(style,arg[2]); // set child class defaults scalar_flag = vector_flag = array_flag = 0; peratom_flag = local_flag = 0; tempflag = pressflag = peflag = 0; pressatomflag = peatomflag = 0; tempbias = 0; timeflag = 0; comm_forward = comm_reverse = 0; + dynamic = 0; dynamic_group_allow = 1; cudable = 0; invoked_scalar = invoked_vector = invoked_array = -1; invoked_peratom = invoked_local = -1; invoked_flag = 0; // set modify defaults extra_dof = domain->dimension; - dynamic = 0; + dynamic_user = 0; // setup list of timesteps ntime = maxtime = 0; tlist = NULL; // setup map for molecule IDs molmap = NULL; datamask = ALL_MASK; datamask_ext = ALL_MASK; } /* ---------------------------------------------------------------------- */ Compute::~Compute() { delete [] id; delete [] style; memory->destroy(tlist); memory->destroy(molmap); } /* ---------------------------------------------------------------------- */ void Compute::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal compute_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"extra") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); extra_dof = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"dynamic") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); - if (strcmp(arg[iarg+1],"no") == 0) dynamic = 0; - else if (strcmp(arg[iarg+1],"yes") == 0) dynamic = 1; + if (strcmp(arg[iarg+1],"no") == 0) dynamic_user = 0; + else if (strcmp(arg[iarg+1],"yes") == 0) dynamic_user = 1; else error->all(FLERR,"Illegal compute_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"thermo") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute_modify command"); if (strcmp(arg[iarg+1],"no") == 0) thermoflag = 0; else if (strcmp(arg[iarg+1],"yes") == 0) thermoflag = 1; else error->all(FLERR,"Illegal compute_modify command"); iarg += 2; } else error->all(FLERR,"Illegal compute_modify command"); } } /* ---------------------------------------------------------------------- calculate adjustment in DOF due to fixes ------------------------------------------------------------------------- */ void Compute::adjust_dof_fix() { fix_dof = 0; for (int i = 0; i < modify->nfix; i++) fix_dof += modify->fix[i]->dof(igroup); } /* ---------------------------------------------------------------------- reset extra_dof to its default value ------------------------------------------------------------------------- */ void Compute::reset_extra_dof() { extra_dof = domain->dimension; } /* ---------------------------------------------------------------------- */ void Compute::reset_extra_compute_fix(const char *) { error->all(FLERR, "Compute does not allow an extra compute or fix to be reset"); } /* ---------------------------------------------------------------------- add ntimestep to list of timesteps the compute will be called on do not add if already in list search from top downward, since list of times is in decreasing order ------------------------------------------------------------------------- */ void Compute::addstep(bigint ntimestep) { // i = location in list to insert ntimestep int i; for (i = ntime-1; i >= 0; i--) { if (ntimestep == tlist[i]) return; if (ntimestep < tlist[i]) break; } i++; // extend list as needed if (ntime == maxtime) { maxtime += DELTA; memory->grow(tlist,maxtime,"compute:tlist"); } // move remainder of list upward and insert ntimestep for (int j = ntime-1; j >= i; j--) tlist[j+1] = tlist[j]; tlist[i] = ntimestep; ntime++; } /* ---------------------------------------------------------------------- return 1/0 if ntimestep is or is not in list of calling timesteps if value(s) on top of list are less than ntimestep, delete them search from top downward, since list of times is in decreasing order ------------------------------------------------------------------------- */ int Compute::matchstep(bigint ntimestep) { for (int i = ntime-1; i >= 0; i--) { if (ntimestep < tlist[i]) return 0; if (ntimestep == tlist[i]) return 1; if (ntimestep > tlist[i]) ntime--; } return 0; } /* ---------------------------------------------------------------------- clean out list of timesteps to call the compute on ------------------------------------------------------------------------- */ void Compute::clearstep() { ntime = 0; } /* ---------------------------------------------------------------------- identify molecule IDs with atoms in group warn if any atom in group has molecule ID = 0 warn if any molecule has only some atoms in group return Ncount = # of molecules with atoms in group set molmap to NULL if molecule IDs include all in range from 1 to Ncount else: molecule IDs range from idlo to idhi set molmap to vector of length idhi-idlo+1 molmap[id-idlo] = index from 0 to Ncount-1 return idlo and idhi ------------------------------------------------------------------------- */ int Compute::molecules_in_group(tagint &idlo, tagint &idhi) { int i; memory->destroy(molmap); molmap = NULL; // find lo/hi molecule ID for any atom in group // warn if atom in group has ID = 0 tagint *molecule = atom->molecule; int *mask = atom->mask; int nlocal = atom->nlocal; tagint lo = BIG; tagint hi = -BIG; int flag = 0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (molecule[i] == 0) flag = 1; lo = MIN(lo,molecule[i]); hi = MAX(hi,molecule[i]); } int flagall; MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall && comm->me == 0) error->warning(FLERR,"Atom with molecule ID = 0 included in " "compute molecule group"); MPI_Allreduce(&lo,&idlo,1,MPI_LMP_TAGINT,MPI_MIN,world); MPI_Allreduce(&hi,&idhi,1,MPI_LMP_TAGINT,MPI_MAX,world); if (idlo == BIG) return 0; // molmap = vector of length nlen // set to 1 for IDs that appear in group across all procs, else 0 tagint nlen_tag = idhi-idlo+1; if (nlen_tag > MAXSMALLINT) error->all(FLERR,"Too many molecules for compute"); int nlen = (int) nlen_tag; memory->create(molmap,nlen,"compute:molmap"); for (i = 0; i < nlen; i++) molmap[i] = 0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) molmap[molecule[i]-idlo] = 1; int *molmapall; memory->create(molmapall,nlen,"compute:molmapall"); MPI_Allreduce(molmap,molmapall,nlen,MPI_INT,MPI_MAX,world); // nmolecules = # of non-zero IDs in molmap // molmap[i] = index of molecule, skipping molecules not in group with -1 int nmolecules = 0; for (i = 0; i < nlen; i++) if (molmapall[i]) molmap[i] = nmolecules++; else molmap[i] = -1; memory->destroy(molmapall); // warn if any molecule has some atoms in group and some not in group flag = 0; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) continue; if (molecule[i] < idlo || molecule[i] > idhi) continue; if (molmap[molecule[i]-idlo] >= 0) flag = 1; } MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); if (flagall && comm->me == 0) error->warning(FLERR, "One or more compute molecules has atoms not in group"); // if molmap simply stores 1 to Nmolecules, then free it if (idlo == 1 && idhi == nmolecules && nlen == nmolecules) { memory->destroy(molmap); molmap = NULL; } return nmolecules; } diff --git a/src/compute.h b/src/compute.h index b456453bb..9344eb92f 100644 --- a/src/compute.h +++ b/src/compute.h @@ -1,195 +1,196 @@ /* -*- c++ -*- ---------------------------------------------------------- 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_COMPUTE_H #define LMP_COMPUTE_H #include "pointers.h" namespace LAMMPS_NS { class Compute : protected Pointers { public: char *id,*style; int igroup,groupbit; double scalar; // computed global scalar double *vector; // computed global vector double **array; // computed global array double *vector_atom; // computed per-atom vector double **array_atom; // computed per-atom array double *vector_local; // computed local vector double **array_local; // computed local array int scalar_flag; // 0/1 if compute_scalar() function exists int vector_flag; // 0/1 if compute_vector() function exists int array_flag; // 0/1 if compute_array() function exists int size_vector; // length of global vector int size_array_rows; // rows in global array int size_array_cols; // columns in global array int peratom_flag; // 0/1 if compute_peratom() function exists int size_peratom_cols; // 0 = vector, N = columns in peratom array int local_flag; // 0/1 if compute_local() function exists int size_local_rows; // rows in local vector or array int size_local_cols; // 0 = vector, N = columns in local array int extscalar; // 0/1 if global scalar is intensive/extensive int extvector; // 0/1/-1 if global vector is all int/ext/extlist int *extlist; // list of 0/1 int/ext for each vec component int extarray; // 0/1 if global array is all intensive/extensive int tempflag; // 1 if Compute can be used as temperature // must have both compute_scalar, compute_vector int pressflag; // 1 if Compute can be used as pressure (uses virial) // must have both compute_scalar, compute_vector int pressatomflag; // 1 if Compute calculates per-atom virial int peflag; // 1 if Compute calculates PE (uses Force energies) int peatomflag; // 1 if Compute calculates per-atom PE int tempbias; // 0/1 if Compute temp includes self/extra bias int timeflag; // 1 if Compute stores list of timesteps it's called on int ntime; // # of entries in time list int maxtime; // max # of entries time list can hold bigint *tlist; // list of timesteps the Compute is called on int invoked_flag; // non-zero if invoked or accessed this step, 0 if not bigint invoked_scalar; // last timestep on which compute_scalar() was invoked bigint invoked_vector; // ditto for compute_vector() bigint invoked_array; // ditto for compute_array() bigint invoked_peratom; // ditto for compute_peratom() bigint invoked_local; // ditto for compute_local() double dof; // degrees-of-freedom for temperature int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int dynamic_group_allow; // 1 if can be used with dynamic group, else 0 unsigned int datamask; unsigned int datamask_ext; int cudable; // 1 if compute is CUDA-enabled Compute(class LAMMPS *, int, char **); virtual ~Compute(); void modify_params(int, char **); void adjust_dof_fix(); void reset_extra_dof(); virtual void init() = 0; virtual void init_list(int, class NeighList *) {} virtual void setup() {} virtual double compute_scalar() {return 0.0;} virtual void compute_vector() {} virtual void compute_array() {} virtual void compute_peratom() {} virtual void compute_local() {} virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} virtual void unpack_forward_comm(int, int, double *) {} virtual int pack_reverse_comm(int, int, double *) {return 0;} virtual void unpack_reverse_comm(int, int *, double *) {} virtual void dof_remove_pre() {} virtual int dof_remove(int) {return 0;} virtual void remove_bias(int, double *) {} virtual void remove_bias_all() {} virtual void restore_bias(int, double *) {} virtual void restore_bias_all() {} virtual void reset_extra_compute_fix(const char *); void addstep(bigint); int matchstep(bigint); void clearstep(); virtual double memory_usage() {return 0.0;} virtual int unsigned data_mask() {return datamask;} virtual int unsigned data_mask_ext() {return datamask_ext;} protected: int extra_dof; // extra DOF for temperature computes int fix_dof; // DOF due to fixes int dynamic; // recount atoms for temperature computes + int dynamic_user; // user request for temp compute to be dynamic int thermoflag; // 1 if include fix PE for PE computes double vbias[3]; // stored velocity bias for one atom double **vbiasall; // stored velocity bias for all atoms int maxbias; // size of vbiasall array int *molmap; // convert molecule ID to local index int molecules_in_group(tagint &, tagint &); inline int sbmask(int j) { return j >> SBBITS & 3; } // union data struct for packing 32-bit and 64-bit ints into double bufs // see atom_vec.h for documentation union ubuf { double d; int64_t i; ubuf(double arg) : d(arg) {} ubuf(int64_t arg) : i(arg) {} ubuf(int arg) : i(arg) {} }; }; } #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Compute ID must be alphanumeric or underscore characters Self-explanatory. E: Could not find compute group ID Self-explanatory. E: Compute does not allow an extra compute or fix to be reset This is an internal LAMMPS error. Please report it to the developers. W: Atom with molecule ID = 0 included in compute molecule group The group used in a compute command that operates on moleclues includes atoms with no molecule ID. This is probably not what you want. E: Too many molecules for compute The limit is 2^31 = ~2 billion molecules. W: One or more compute molecules has atoms not in group The group used in a compute command that operates on moleclues does not include all the atoms in some molecules. This is probably not what you want. */ diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index c54a9ba9b..6c0fb09c0 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -1,133 +1,135 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "string.h" #include "compute_temp.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTemp::~ComputeTemp() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTemp::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTemp::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0.0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTemp::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTemp::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } diff --git a/src/compute_temp_com.cpp b/src/compute_temp_com.cpp index 0c94d2a41..64d23efcd 100644 --- a/src/compute_temp_com.cpp +++ b/src/compute_temp_com.cpp @@ -1,218 +1,220 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_com.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "domain.h" #include "lattice.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempCOM::~ComputeTempCOM() { delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::init() { masstotal = group->mass(igroup); } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempCOM::compute_scalar() { double vthermal[3]; invoked_scalar = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vbias[0]; vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempCOM::compute_vector() { int i; double vthermal[3]; invoked_vector = update->ntimestep; if (dynamic) masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vbias); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vthermal[0] = v[i][0] - vbias[0]; vthermal[1] = v[i][1] - vbias[1]; vthermal[2] = v[i][2] - vbias[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempCOM::remove_bias(int, double *v) { v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempCOM::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] -= vbias[0]; v[i][1] -= vbias[1]; v[i][2] -= vbias[2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempCOM::restore_bias(int, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempCOM::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbias[0]; v[i][1] += vbias[1]; v[i][2] += vbias[2]; } } diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 7c72e08f3..7acce8334 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -1,291 +1,293 @@ /* ---------------------------------------------------------------------- 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: Pieter in 't Veld (SNL) ------------------------------------------------------------------------- */ #include "mpi.h" #include "string.h" #include "compute_temp_deform.h" #include "domain.h" #include "atom.h" #include "update.h" #include "force.h" #include "modify.h" #include "fix.h" #include "fix_deform.h" #include "group.h" #include "comm.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp /* ---------------------------------------------------------------------- */ ComputeTempDeform::ComputeTempDeform(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal compute temp/deform command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempDeform::~ComputeTempDeform() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::init() { int i; // check fix deform remap settings for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP && comm->me == 0) error->warning(FLERR,"Using compute temp/deform with inconsistent " "fix deform remap option"); break; } if (i == modify->nfix && comm->me == 0) error->warning(FLERR, "Using compute temp/deform with no fix deform defined"); } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempDeform::compute_scalar() { double lamda[3],vstream[3],vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; // lamda = 0-1 triclinic lamda coords // vstream = streaming velocity = Hrate*lamda + Hratelo // vthermal = thermal velocity = v - vstream double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempDeform::compute_vector() { double lamda[3],vstream[3],vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; double massone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(x[i],lamda); vstream[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vstream[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vstream[2] = h_rate[2]*lamda[2] + h_ratelo[2]; vthermal[0] = v[i][0] - vstream[0]; vthermal[1] = v[i][1] - vstream[1]; vthermal[2] = v[i][2] - vstream[2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeform::remove_bias(int i, double *v) { double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; domain->x2lamda(atom->x[i],lamda); vbias[0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbias[1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbias[2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[0] -= vbias[0]; v[1] -= vbias[1]; v[2] -= vbias[2]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempDeform::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/deform:vbiasall"); } double lamda[3]; double *h_rate = domain->h_rate; double *h_ratelo = domain->h_ratelo; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->x2lamda(atom->x[i],lamda); vbiasall[i][0] = h_rate[0]*lamda[0] + h_rate[5]*lamda[1] + h_rate[4]*lamda[2] + h_ratelo[0]; vbiasall[i][1] = h_rate[1]*lamda[1] + h_rate[3]*lamda[2] + h_ratelo[1]; vbiasall[i][2] = h_rate[2]*lamda[2] + h_ratelo[2]; v[i][0] -= vbiasall[i][0]; v[i][1] -= vbiasall[i][1]; v[i][2] -= vbiasall[i][2]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeform::restore_bias(int i, double *v) { v[0] += vbias[0]; v[1] += vbias[1]; v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempDeform::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { v[i][0] += vbiasall[i][0]; v[i][1] += vbiasall[i][1]; v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempDeform::memory_usage() { double bytes = 3*maxbias * sizeof(double); return bytes; } diff --git a/src/compute_temp_partial.cpp b/src/compute_temp_partial.cpp index 4fcb27dd4..19aab68e5 100644 --- a/src/compute_temp_partial.cpp +++ b/src/compute_temp_partial.cpp @@ -1,264 +1,266 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "compute_temp_partial.h" #include "atom.h" #include "update.h" #include "force.h" #include "domain.h" #include "group.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg != 6) error->all(FLERR,"Illegal compute temp/partial command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; xflag = force->inumeric(FLERR,arg[3]); yflag = force->inumeric(FLERR,arg[4]); zflag = force->inumeric(FLERR,arg[5]); if (zflag && domain->dimension == 2) error->all(FLERR,"Compute temp/partial cannot use vz for 2d systemx"); maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempPartial::~ComputeTempPartial() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempPartial::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- DOF for a body of N atoms with S constraints (e.g. from SHAKE) DOF = nper/dim (dim*N - S), where dim = dimensionality = 2 or 3 ------------------------------------------------------------------------- */ void ComputeTempPartial::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = xflag+yflag+zflag; dof = nper * natoms; dof -= (1.0*nper/domain->dimension)*fix_dof + extra_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ int ComputeTempPartial::dof_remove(int i) { int nper = xflag+yflag+zflag; return (domain->dimension - nper); } /* ---------------------------------------------------------------------- */ double ComputeTempPartial::compute_scalar() { invoked_scalar = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * rmass[i]; } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] + zflag*v[i][2]*v[i][2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempPartial::compute_vector() { int i; invoked_vector = update->ntimestep; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * xflag*v[i][0]*v[i][0]; t[1] += massone * yflag*v[i][1]*v[i][1]; t[2] += massone * zflag*v[i][2]*v[i][2]; t[3] += massone * xflag*yflag*v[i][0]*v[i][1]; t[4] += massone * xflag*zflag*v[i][0]*v[i][2]; t[5] += massone * yflag*zflag*v[i][1]*v[i][2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempPartial::remove_bias(int i, double *v) { if (!xflag) { vbias[0] = v[0]; v[0] = 0.0; } if (!yflag) { vbias[1] = v[1]; v[1] = 0.0; } if (!zflag) { vbias[2] = v[2]; v[2] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempPartial::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/partial:vbiasall"); } if (!xflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][0] = v[i][0]; v[i][0] = 0.0; } } if (!yflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][1] = v[i][1]; v[i][1] = 0.0; } } if (!zflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { vbiasall[i][2] = v[i][2]; v[i][2] = 0.0; } } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempPartial::restore_bias(int i, double *v) { if (!xflag) v[0] += vbias[0]; if (!yflag) v[1] += vbias[1]; if (!zflag) v[2] += vbias[2]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempPartial::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (!xflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][0] += vbiasall[i][0]; } if (!yflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][1] += vbiasall[i][1]; } if (!zflag) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][2] += vbiasall[i][2]; } } /* ---------------------------------------------------------------------- */ double ComputeTempPartial::memory_usage() { double bytes = 3*maxbias * sizeof(double); return bytes; } diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index d1e107a82..58afbb3a5 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -1,563 +1,565 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_profile.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "modify.h" #include "fix.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; enum{TENSOR,BIN}; /* ---------------------------------------------------------------------- */ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 7) error->all(FLERR,"Illegal compute temp/profile command"); scalar_flag = 1; extscalar = 0; tempflag = 1; tempbias = 1; xflag = force->inumeric(FLERR,arg[3]); yflag = force->inumeric(FLERR,arg[4]); zflag = force->inumeric(FLERR,arg[5]); if (zflag && domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot use vz for 2d systemx"); ncount = 0; ivx = ivy = ivz = 0; if (xflag) ivx = ncount++; if (yflag) ivy = ncount++; if (zflag) ivz = ncount++; ncount += 2; nbinx = nbiny = nbinz = 1; int iarg = 6; if (strcmp(arg[iarg],"x") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbinx = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"y") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbiny = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"z") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinz = force->inumeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"xy") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbiny = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"yz") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbiny = force->inumeric(FLERR,arg[iarg+1]); nbinz = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"xz") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbinz = force->inumeric(FLERR,arg[iarg+2]); iarg += 3; } else if (strcmp(arg[iarg],"xyz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (domain->dimension == 2) error->all(FLERR,"Compute temp/profile cannot bin z for 2d systems"); nbinx = force->inumeric(FLERR,arg[iarg+1]); nbiny = force->inumeric(FLERR,arg[iarg+2]); nbinz = force->inumeric(FLERR,arg[iarg+3]); iarg += 4; } else error->all(FLERR,"Illegal compute temp/profile command"); // optional keywords outflag = TENSOR; while (iarg < narg) { if (strcmp(arg[iarg],"out") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/profile command"); if (strcmp(arg[iarg+1],"tensor") == 0) outflag = TENSOR; else if (strcmp(arg[iarg+1],"bin") == 0) outflag = BIN; else error->all(FLERR,"Illegal compute temp/profile command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/profile command"); } // setup nbins = nbinx*nbiny*nbinz; if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command"); memory->create(vbin,nbins,ncount,"temp/profile:vbin"); memory->create(binave,nbins,ncount,"temp/profile:binave"); if (outflag == TENSOR) { vector_flag = 1; size_vector = 6; extvector = 1; vector = new double[size_vector]; } else { array_flag = 1; size_array_rows = nbins; size_array_cols = 2; extarray = 0; memory->create(tbin,nbins,"temp/profile:tbin"); memory->create(tbinall,nbins,"temp/profile:tbinall"); memory->create(array,nbins,2,"temp/profile:array"); } maxatom = 0; bin = NULL; } /* ---------------------------------------------------------------------- */ ComputeTempProfile::~ComputeTempProfile() { memory->destroy(vbin); memory->destroy(binave); memory->destroy(bin); if (outflag == TENSOR) delete [] vector; else { memory->destroy(tbin); memory->destroy(tbinall); memory->destroy(array); } } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::init() { fix_dof = -1; dof_compute(); // ptrs to domain data box_change = domain->box_change; triclinic = domain->triclinic; periodicity = domain->periodicity; if (triclinic) { boxlo = domain->boxlo_lamda; boxhi = domain->boxhi_lamda; prd = domain->prd_lamda; } else { boxlo = domain->boxlo; boxhi = domain->boxhi; prd = domain->prd; } if (!box_change) bin_setup(); } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempProfile::compute_scalar() { int ibin; double vthermal[3]; invoked_scalar = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::compute_vector() { int i,ibin; double vthermal[3]; invoked_vector = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- */ void ComputeTempProfile::compute_array() { int i,ibin; double vthermal[3]; invoked_array = update->ntimestep; bin_average(); double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; for (i = 0; i < nbins; i++) tbin[i] = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vthermal[0] = v[i][0] - binave[ibin][ivx]; else vthermal[0] = v[i][0]; if (yflag) vthermal[1] = v[i][1] - binave[ibin][ivy]; else vthermal[1] = v[i][1]; if (zflag) vthermal[2] = v[i][2] - binave[ibin][ivz]; else vthermal[2] = v[i][2]; if (rmass) tbin[ibin] += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else tbin[ibin] += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(tbin,tbinall,nbins,MPI_DOUBLE,MPI_SUM,world); int nper = domain->dimension; for (i = 0; i < nbins; i++) { array[i][0] = binave[i][ncount-1]; if (array[i][0] > 0.0) { dof = nper*array[i][0] - extra_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; array[i][1] = tfactor*tbinall[i]; } else array[i][1] = 0.0; } } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempProfile::remove_bias(int i, double *v) { int ibin = bin[i]; if (xflag) v[0] -= binave[ibin][ivx]; if (yflag) v[1] -= binave[ibin][ivy]; if (zflag) v[2] -= binave[ibin][ivz]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempProfile::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; int ibin; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) v[i][0] -= binave[ibin][ivx]; if (yflag) v[i][1] -= binave[ibin][ivy]; if (zflag) v[i][2] -= binave[ibin][ivz]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempProfile::restore_bias(int i, double *v) { int ibin = bin[i]; if (xflag) v[0] += binave[ibin][ivx]; if (yflag) v[1] += binave[ibin][ivy]; if (zflag) v[2] += binave[ibin][ivz]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempProfile::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; int ibin; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) v[i][0] += binave[ibin][ivx]; if (yflag) v[i][1] += binave[ibin][ivy]; if (zflag) v[i][2] += binave[ibin][ivz]; } } /* ---------------------------------------------------------------------- compute average COM velocity in each bin ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_average() { int i,j,ibin; if (box_change) bin_setup(); bin_assign(); // clear bins, including particle mass and count for (i = 0; i < nbins; i++) for (j = 0; j < ncount; j++) vbin[i][j] = 0.0; // sum each particle's mass-weighted velocity, mass, count to appropriate bin double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *mask = atom->mask; int *type = atom->type; int nlocal = atom->nlocal; int nc2 = ncount-2; int nc1 = ncount-1; if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; if (xflag) vbin[ibin][ivx] += rmass[i]*v[i][0]; if (yflag) vbin[ibin][ivy] += rmass[i]*v[i][1]; if (zflag) vbin[ibin][ivz] += rmass[i]*v[i][2]; vbin[ibin][nc2] += rmass[i]; vbin[ibin][nc1] += 1.0; } } else { double onemass; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { ibin = bin[i]; onemass = mass[type[i]]; if (xflag) vbin[ibin][ivx] += onemass*v[i][0]; if (yflag) vbin[ibin][ivy] += onemass*v[i][1]; if (zflag) vbin[ibin][ivz] += onemass*v[i][2]; vbin[ibin][nc2] += onemass; vbin[ibin][nc1] += 1.0; } } // sum bins across processors MPI_Allreduce(vbin[0],binave[0],nbins*ncount,MPI_DOUBLE,MPI_SUM,world); // compute ave COM velocity in each bin, checking for no particles for (i = 0; i < nbins; i++) if (binave[i][nc1] > 0.0) for (j = 0; j < nc2; j++) binave[i][j] /= binave[i][nc2]; } /* ---------------------------------------------------------------------- set bin sizes, redo if box size changes ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_setup() { invdelta[0] = nbinx / prd[0]; invdelta[1] = nbiny / prd[1]; invdelta[2] = nbinz / prd[2]; } /* ---------------------------------------------------------------------- assign all atoms to bins ------------------------------------------------------------------------- */ void ComputeTempProfile::bin_assign() { // reallocate bin array if necessary if (atom->nlocal > maxatom) { maxatom = atom->nmax; memory->destroy(bin); memory->create(bin,maxatom,"temp/profile:bin"); } // assign each atom to a bin, accounting for PBC // if triclinic, do this in lamda space double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; int ibinx,ibiny,ibinz; double coord; if (triclinic) domain->x2lamda(nlocal); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (nbinx > 1) { coord = x[i][0]; if (periodicity[0]) { if (coord < boxlo[0]) coord += prd[0]; if (coord >= boxhi[0]) coord -= prd[0]; } ibinx = static_cast ((coord - boxlo[0]) * invdelta[0]); ibinx = MAX(ibinx,0); ibinx = MIN(ibinx,nbinx-1); } else ibinx = 0; if (nbiny > 1) { coord = x[i][1]; if (periodicity[1]) { if (coord < boxlo[1]) coord += prd[1]; if (coord >= boxhi[1]) coord -= prd[1]; } ibiny = static_cast ((coord - boxlo[1]) * invdelta[1]); ibiny = MAX(ibiny,0); ibiny = MIN(ibiny,nbiny-1); } else ibiny = 0; if (nbinz > 1) { coord = x[i][2]; if (periodicity[2]) { if (coord < boxlo[2]) coord += prd[2]; if (coord >= boxhi[2]) coord -= prd[2]; } ibinz = static_cast ((coord - boxlo[2]) * invdelta[2]); ibinz = MAX(ibinz,0); ibinz = MIN(ibinz,nbinz-1); } else ibinz = 0; bin[i] = nbinx*nbiny*ibinz + nbinx*ibiny + ibinx; } if (triclinic) domain->lamda2x(nlocal); } /* ---------------------------------------------------------------------- */ double ComputeTempProfile::memory_usage() { double bytes = maxatom * sizeof(int); bytes += nbins*ncount * sizeof(double); return bytes; } diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index 33a89229f..c4b1ee4dc 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -1,295 +1,297 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "stdlib.h" #include "string.h" #include "compute_temp_ramp.h" #include "atom.h" #include "update.h" #include "force.h" #include "group.h" #include "modify.h" #include "fix.h" #include "domain.h" #include "lattice.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 9) error->all(FLERR,"Illegal compute temp command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 1; // parse optional args scaleflag = 1; int iarg = 9; while (iarg < narg) { if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/ramp command"); if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; else error->all(FLERR,"Illegal compute temp/ramp command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/ramp command"); } // setup scaling if (scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; // read standard args and apply scaling if (strcmp(arg[3],"vx") == 0) v_dim = 0; else if (strcmp(arg[3],"vy") == 0) v_dim = 1; else if (strcmp(arg[3],"vz") == 0) v_dim = 2; else error->all(FLERR,"Illegal compute temp/ramp command"); if (v_dim == 0) { v_lo = xscale*force->numeric(FLERR,arg[4]); v_hi = xscale*force->numeric(FLERR,arg[5]); } else if (v_dim == 1) { v_lo = yscale*force->numeric(FLERR,arg[4]); v_hi = yscale*force->numeric(FLERR,arg[5]); } else if (v_dim == 2) { v_lo = zscale*force->numeric(FLERR,arg[4]); v_hi = zscale*force->numeric(FLERR,arg[5]); } if (strcmp(arg[6],"x") == 0) coord_dim = 0; else if (strcmp(arg[6],"y") == 0) coord_dim = 1; else if (strcmp(arg[6],"z") == 0) coord_dim = 2; else error->all(FLERR,"Illegal compute temp/ramp command"); if (coord_dim == 0) { coord_lo = xscale*force->numeric(FLERR,arg[7]); coord_hi = xscale*force->numeric(FLERR,arg[8]); } else if (coord_dim == 1) { coord_lo = yscale*force->numeric(FLERR,arg[7]); coord_hi = yscale*force->numeric(FLERR,arg[8]); } else if (coord_dim == 2) { coord_lo = zscale*force->numeric(FLERR,arg[7]); coord_hi = zscale*force->numeric(FLERR,arg[8]); } maxbias = 0; vbiasall = NULL; vector = new double[6]; } /* ---------------------------------------------------------------------- */ ComputeTempRamp::~ComputeTempRamp() { memory->destroy(vbiasall); delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::dof_compute() { if (fix_dof) adjust_dof_fix(); double natoms = group->count(igroup); int nper = domain->dimension; dof = nper * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempRamp::compute_scalar() { double fraction,vramp,vthermal[3]; invoked_scalar = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double t = 0.0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vramp = v_lo + fraction*(v_hi - v_lo); vthermal[0] = v[i][0]; vthermal[1] = v[i][1]; vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; if (rmass) t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * rmass[i]; else t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + vthermal[2]*vthermal[2]) * mass[type[i]]; } MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempRamp::compute_vector() { int i; double fraction,vramp,vthermal[3]; invoked_vector = update->ntimestep; double **x = atom->x; double **v = atom->v; double *mass = atom->mass; double *rmass = atom->rmass; int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; double massone,t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vramp = v_lo + fraction*(v_hi - v_lo); vthermal[0] = v[i][0]; vthermal[1] = v[i][1]; vthermal[2] = v[i][2]; vthermal[v_dim] -= vramp; if (rmass) massone = rmass[i]; else massone = mass[type[i]]; t[0] += massone * vthermal[0]*vthermal[0]; t[1] += massone * vthermal[1]*vthermal[1]; t[2] += massone * vthermal[2]*vthermal[2]; t[3] += massone * vthermal[0]*vthermal[1]; t[4] += massone * vthermal[0]*vthermal[2]; t[5] += massone * vthermal[1]*vthermal[2]; } MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRamp::remove_bias(int i, double *v) { double fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vbias[v_dim] = v_lo + fraction*(v_hi - v_lo); v[v_dim] -= vbias[v_dim]; } /* ---------------------------------------------------------------------- remove velocity bias from all atoms to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempRamp::remove_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (nlocal > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; memory->create(vbiasall,maxbias,3,"temp/ramp:vbiasall"); } double fraction; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { fraction = (atom->x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo); fraction = MAX(fraction,0.0); fraction = MIN(fraction,1.0); vbiasall[i][v_dim] = v_lo + fraction*(v_hi - v_lo); v[i][v_dim] -= vbiasall[i][v_dim]; } } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempRamp::restore_bias(int i, double *v) { v[v_dim] += vbias[v_dim]; } /* ---------------------------------------------------------------------- add back in velocity bias to all atoms removed by remove_bias_all() assume remove_bias_all() was previously called ------------------------------------------------------------------------- */ void ComputeTempRamp::restore_bias_all() { double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) v[i][v_dim] += vbiasall[i][v_dim]; } /* ---------------------------------------------------------------------- */ double ComputeTempRamp::memory_usage() { double bytes = 3*maxbias * sizeof(double); return bytes; } diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 8ba73f3c5..5e5a9c88e 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -1,331 +1,333 @@ /* ---------------------------------------------------------------------- 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 "mpi.h" #include "string.h" #include "compute_temp_sphere.h" #include "atom.h" #include "atom_vec.h" #include "update.h" #include "force.h" #include "domain.h" #include "modify.h" #include "group.h" #include "error.h" using namespace LAMMPS_NS; enum{ROTATE,ALL}; #define INERTIA 0.4 // moment of inertia prefactor for sphere /* ---------------------------------------------------------------------- */ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { if (narg < 3) error->all(FLERR,"Illegal compute temp/sphere command"); scalar_flag = vector_flag = 1; size_vector = 6; extscalar = 0; extvector = 1; tempflag = 1; tempbias = 0; id_bias = NULL; mode = ALL; int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"bias") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/sphere command"); tempbias = 1; int n = strlen(arg[iarg+1]) + 1; id_bias = new char[n]; strcpy(id_bias,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"dof") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute temp/sphere command"); if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE; else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL; else error->all(FLERR,"Illegal compute temp/sphere command"); iarg += 2; } else error->all(FLERR,"Illegal compute temp/sphere command"); } vector = new double[6]; // error checks if (!atom->sphere_flag) error->all(FLERR,"Compute temp/sphere requires atom style sphere"); } /* ---------------------------------------------------------------------- */ ComputeTempSphere::~ComputeTempSphere() { delete [] id_bias; delete [] vector; } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::init() { if (tempbias) { int i = modify->find_compute(id_bias); if (i < 0) error->all(FLERR,"Could not find compute ID for temperature bias"); tbias = modify->compute[i]; if (tbias->tempflag == 0) error->all(FLERR,"Bias compute does not calculate temperature"); if (tbias->tempbias == 0) error->all(FLERR,"Bias compute does not calculate a velocity bias"); if (tbias->igroup != igroup) error->all(FLERR,"Bias compute group does not match compute group"); tbias->init(); tbias->setup(); if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; else tempbias = 1; } } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::setup() { + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; fix_dof = -1; dof_compute(); } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::dof_compute() { int count,count_all; if (fix_dof) adjust_dof_fix(); // 6 or 3 dof for extended/point particles for 3d // 3 or 2 dof for extended/point particles for 2d // which dof are included also depends on mode // assume full rotation of extended particles // user should correct this via compute_modify if needed double *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; count = 0; if (domain->dimension == 3) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (radius[i] == 0.0) { if (mode == ALL) count += 3; } else { if (mode == ALL) count += 6; else count += 3; } } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (radius[i] == 0.0) { if (mode == ALL) count += 2; } else { if (mode == ALL) count += 3; else count += 1; } } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof = count_all; // additional adjustments to dof if (tempbias == 1) { if (mode == ALL) { double natoms = group->count(igroup); dof -= tbias->dof_remove(-1) * natoms; } } else if (tempbias == 2) { int *mask = atom->mask; int nlocal = atom->nlocal; tbias->dof_remove_pre(); count = 0; if (domain->dimension == 3) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (tbias->dof_remove(i)) { if (radius[i] == 0.0) { if (mode == ALL) count += 3; } else { if (mode == ALL) count += 6; else count += 3; } } } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { if (tbias->dof_remove(i)) { if (radius[i] == 0.0) { if (mode == ALL) count += 2; } else { if (mode == ALL) count += 3; else count += 1; } } } } MPI_Allreduce(&count,&count_all,1,MPI_INT,MPI_SUM,world); dof -= count_all; } dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } /* ---------------------------------------------------------------------- */ double ComputeTempSphere::compute_scalar() { invoked_scalar = update->ntimestep; if (tempbias) { if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); tbias->remove_bias_all(); } double **v = atom->v; double **omega = atom->omega; double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; // point particles will not contribute rotation due to radius = 0 double t = 0.0; if (mode == ALL) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i]; t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i]; } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world); if (dynamic || tempbias == 2) dof_compute(); scalar *= tfactor; return scalar; } /* ---------------------------------------------------------------------- */ void ComputeTempSphere::compute_vector() { invoked_vector = update->ntimestep; if (tempbias) { if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); tbias->remove_bias_all(); } double **v = atom->v; double **omega = atom->omega; double *rmass = atom->rmass; double *radius = atom->radius; int *mask = atom->mask; int nlocal = atom->nlocal; // point particles will not contribute rotation due to radius = 0 double massone,inertiaone,t[6]; for (int i = 0; i < 6; i++) t[i] = 0.0; if (mode == ALL) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { massone = rmass[i]; t[0] += massone * v[i][0]*v[i][0]; t[1] += massone * v[i][1]*v[i][1]; t[2] += massone * v[i][2]*v[i][2]; t[3] += massone * v[i][0]*v[i][1]; t[4] += massone * v[i][0]*v[i][2]; t[5] += massone * v[i][1]*v[i][2]; inertiaone = INERTIA*rmass[i]*radius[i]*radius[i]; t[0] += inertiaone * omega[i][0]*omega[i][0]; t[1] += inertiaone * omega[i][1]*omega[i][1]; t[2] += inertiaone * omega[i][2]*omega[i][2]; t[3] += inertiaone * omega[i][0]*omega[i][1]; t[4] += inertiaone * omega[i][0]*omega[i][2]; t[5] += inertiaone * omega[i][1]*omega[i][2]; } } else { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { inertiaone = INERTIA*rmass[i]*radius[i]*radius[i]; t[0] += inertiaone * omega[i][0]*omega[i][0]; t[1] += inertiaone * omega[i][1]*omega[i][1]; t[2] += inertiaone * omega[i][2]*omega[i][2]; t[3] += inertiaone * omega[i][0]*omega[i][1]; t[4] += inertiaone * omega[i][0]*omega[i][2]; t[5] += inertiaone * omega[i][1]*omega[i][2]; } } if (tempbias) tbias->restore_bias_all(); MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); for (int i = 0; i < 6; i++) vector[i] *= force->mvv2e; } /* ---------------------------------------------------------------------- remove velocity bias from atom I to leave thermal velocity ------------------------------------------------------------------------- */ void ComputeTempSphere::remove_bias(int i, double *v) { tbias->remove_bias(i,v); } /* ---------------------------------------------------------------------- add back in velocity bias to atom I removed by remove_bias() assume remove_bias() was previously called ------------------------------------------------------------------------- */ void ComputeTempSphere::restore_bias(int i, double *v) { tbias->restore_bias(i,v); } diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index f7a8d22a3..e9b23a1c9 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -1,2294 +1,2298 @@ /* ---------------------------------------------------------------------- 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 "math_extra.h" #include "atom.h" #include "force.h" #include "group.h" #include "comm.h" #include "neighbor.h" #include "irregular.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" using namespace LAMMPS_NS; using namespace FixConst; #define DELTAFLIP 0.1 #define TILTMAX 1.5 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(FLERR,"Illegal fix nvt/npt/nph command"); restart_global = 1; dynamic_group_allow = 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; id_dilate = NULL; mtchain = mpchain = 3; nc_tchain = nc_pchain = 1; mtk_flag = 1; deviatoric_flag = 0; nreset_h0 = 0; eta_mass_flag = 1; omega_mass_flag = 0; etap_mass_flag = 0; flipflag = 1; // turn on tilt factor scaling, whenever applicable dimension = domain->dimension; scaleyz = scalexz = scalexy = 0; if (domain->yperiodic && domain->xy != 0.0) scalexy = 1; if (domain->zperiodic && dimension == 3) { if (domain->yz != 0.0) scaleyz = 1; if (domain->xz != 0.0) scalexz = 1; } // set fixed-point to default = center of cell fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]); fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]); fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]); // used by FixNVTSllod to preserve non-default value mtchain_default_flag = 1; 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] = p_target[i] = 0.0; p_flag[i] = 0; } // process keywords int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"temp") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); tstat_flag = 1; t_start = force->numeric(FLERR,arg[iarg+1]); t_target = t_start; t_stop = force->numeric(FLERR,arg[iarg+2]); t_period = force->numeric(FLERR,arg[iarg+3]); if (t_start < 0.0 || t_stop <= 0.0) error->all(FLERR, "Target temperature for fix nvt/npt/nph cannot be 0.0"); iarg += 4; } else if (strcmp(arg[iarg],"iso") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = XYZ; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,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(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,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(FLERR,"Illegal fix nvt/npt/nph command"); pcouple = NONE; scalexy = scalexz = scaleyz = 0; p_start[0] = p_start[1] = p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = p_stop[1] = p_stop[2] = force->numeric(FLERR,arg[iarg+2]); - p_period[0] = p_period[1] = p_period[2] = force->numeric(FLERR,arg[iarg+3]); + p_period[0] = p_period[1] = p_period[2] = + force->numeric(FLERR,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] = force->numeric(FLERR,arg[iarg+3]); + p_period[3] = p_period[4] = p_period[5] = + force->numeric(FLERR,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(FLERR,"Illegal fix nvt/npt/nph command"); p_start[0] = force->numeric(FLERR,arg[iarg+1]); p_stop[0] = force->numeric(FLERR,arg[iarg+2]); p_period[0] = force->numeric(FLERR,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(FLERR,"Illegal fix nvt/npt/nph command"); p_start[1] = force->numeric(FLERR,arg[iarg+1]); p_stop[1] = force->numeric(FLERR,arg[iarg+2]); p_period[1] = force->numeric(FLERR,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(FLERR,"Illegal fix nvt/npt/nph command"); p_start[2] = force->numeric(FLERR,arg[iarg+1]); p_stop[2] = force->numeric(FLERR,arg[iarg+2]); p_period[2] = force->numeric(FLERR,arg[iarg+3]); p_flag[2] = 1; deviatoric_flag = 1; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"yz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[3] = force->numeric(FLERR,arg[iarg+1]); p_stop[3] = force->numeric(FLERR,arg[iarg+2]); p_period[3] = force->numeric(FLERR,arg[iarg+3]); p_flag[3] = 1; deviatoric_flag = 1; scaleyz = 0; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xz") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[4] = force->numeric(FLERR,arg[iarg+1]); p_stop[4] = force->numeric(FLERR,arg[iarg+2]); p_period[4] = force->numeric(FLERR,arg[iarg+3]); p_flag[4] = 1; deviatoric_flag = 1; scalexz = 0; iarg += 4; if (dimension == 2) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); } else if (strcmp(arg[iarg],"xy") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); p_start[5] = force->numeric(FLERR,arg[iarg+1]); p_stop[5] = force->numeric(FLERR,arg[iarg+2]); p_period[5] = force->numeric(FLERR,arg[iarg+3]); p_flag[5] = 1; deviatoric_flag = 1; scalexy = 0; iarg += 4; } else if (strcmp(arg[iarg],"couple") == 0) { if (iarg+2 > narg) error->all(FLERR,"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(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"drag") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); drag = force->numeric(FLERR,arg[iarg+1]); if (drag < 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"dilate") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"all") == 0) allremap = 1; else { allremap = 0; delete [] id_dilate; int n = strlen(arg[iarg+1]) + 1; id_dilate = new char[n]; strcpy(id_dilate,arg[iarg+1]); int idilate = group->find(id_dilate); if (idilate == -1) error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist"); } iarg += 2; } else if (strcmp(arg[iarg],"tchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); mtchain = force->inumeric(FLERR,arg[iarg+1]); // used by FixNVTSllod to preserve non-default value mtchain_default_flag = 0; if (mtchain < 1) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"pchain") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); mpchain = force->inumeric(FLERR,arg[iarg+1]); if (mpchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"mtk") == 0) { if (iarg+2 > narg) error->all(FLERR,"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(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"tloop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nc_tchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_tchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"ploop") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nc_pchain = force->inumeric(FLERR,arg[iarg+1]); if (nc_pchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"nreset") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); nreset_h0 = force->inumeric(FLERR,arg[iarg+1]); if (nreset_h0 < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scalexy") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1; else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scalexz") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1; else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"scaleyz") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1; else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"flip") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1; else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0; else error->all(FLERR,"Illegal fix nvt/npt/nph command"); iarg += 2; } else if (strcmp(arg[iarg],"fixedpoint") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command"); fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]); fixedpoint[1] = force->numeric(FLERR,arg[iarg+2]); fixedpoint[2] = force->numeric(FLERR,arg[iarg+3]); iarg += 4; } else error->all(FLERR,"Illegal fix nvt/npt/nph command"); } // error checks if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4])) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (dimension == 2 && (pcouple == YZ || pcouple == XZ)) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (dimension == 2 && (scalexz == 1 || scaleyz == 1 )) error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation"); if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0)) error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings"); // require periodicity in tensile dimension if (p_flag[0] && domain->xperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); if (p_flag[1] && domain->yperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); if (p_flag[2] && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension"); // require periodicity in 2nd dim of off-diagonal tilt component if (p_flag[3] && domain->zperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (p_flag[4] && domain->zperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (p_flag[5] && domain->yperiodic == 0) error->all(FLERR, "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension"); if (scaleyz == 1 && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with yz scaling when z is non-periodic dimension"); if (scalexz == 1 && domain->zperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with xz scaling when z is non-periodic dimension"); if (scalexy == 1 && domain->yperiodic == 0) error->all(FLERR,"Cannot use fix nvt/npt/nph " "with xy scaling when y is non-periodic dimension"); if (p_flag[3] && scaleyz == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both yz dynamics and yz scaling"); if (p_flag[4] && scalexz == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both xz dynamics and xz scaling"); if (p_flag[5] && scalexy == 1) error->all(FLERR,"Cannot use fix nvt/npt/nph with " "both xy dynamics and xy scaling"); if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5])) error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in " "fix nvt/npt/nph with non-triclinic box"); if (pcouple == XYZ && dimension == 3 && (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(FLERR,"Invalid fix nvt/npt/nph pressure settings"); if (pcouple == XYZ && dimension == 2 && (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] || p_period[0] != p_period[1])) error->all(FLERR,"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(FLERR,"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(FLERR,"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(FLERR,"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(FLERR,"Fix nvt/npt/nph damping parameters must be > 0.0"); // set pstat_flag and box change and restart_pbc variables pstat_flag = 0; for (int i = 0; i < 6; i++) if (p_flag[i]) pstat_flag = 1; if (pstat_flag) { if (p_flag[0] || p_flag[1] || p_flag[2]) box_change_size = 1; if (p_flag[3] || p_flag[4] || p_flag[5]) box_change_shape = 1; no_change_box = 1; if (allremap == 0) restart_pbc = 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; // pre_exchange only required if flips can occur due to shape changes pre_exchange_flag = 0; if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) pre_exchange_flag = 1; if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0)) pre_exchange_flag = 1; // 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; if (pre_exchange_flag) irregular = new Irregular(lmp); else irregular = NULL; // initialize vol0,t0 to zero to signal uninitialized // values then assigned in init(), if necessary vol0 = t0 = 0.0; } /* ---------------------------------------------------------------------- */ FixNH::~FixNH() { delete [] id_dilate; delete [] rfix; delete irregular; // 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; if (pre_exchange_flag) mask |= PRE_EXCHANGE; return mask; } /* ---------------------------------------------------------------------- */ void FixNH::init() { // recheck that dilate group has not been deleted if (allremap == 0) { int idilate = group->find(id_dilate); if (idilate == -1) error->all(FLERR,"Fix nvt/npt/nph dilate group ID does not exist"); dilate_group_bit = group->bitmask[idilate]; } // ensure 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(FLERR,"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(FLERR,"Temperature ID for fix nvt/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(FLERR,"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; dto = dthalf; 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]); } pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); } if (tstat_flag) tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); // tally the number of dimensions that are barostatted // 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 (strstr(update->integrate_style,"respa")) { nlevels_respa = ((Respa *) update->integrate)->nlevels; step_respa = ((Respa *) update->integrate)->step; dto = 0.5*step_respa[0]; } // 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 - - tdof = temperature->dof; - // t_target is needed by NPH and NPT in compute_scalar() // If no thermostat or using fix nphug, // t_target must be defined by other means. if (tstat_flag && strcmp(style,"nphug") != 0) { compute_temp_target(); } else if (pstat_flag) { // t0 = reference temperature for masses // 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(); + tdof = temperature->dof; + if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); else pressure->compute_vector(); couple(); pressure->addstep(update->ntimestep+1); } // masses and 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]; } } // masses and initial forces on barostat variables 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]); } // masses and 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]; } } } /* ---------------------------------------------------------------------- 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) { compute_temp_target(); 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(); pressure->compute_scalar(); } else { temperature->compute_vector(); pressure->compute_vector(); } couple(); pressure->addstep(update->ntimestep+1); } 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(); // re-compute temp before nh_v_press() // only needed for temperature computes with BIAS on reneighboring steps: // b/c some biases store per-atom values (e.g. temp/profile) // per-atom values are invalid if reneigh/comm occurred // since temp->compute() in initial_integrate() if (which == BIAS && neighbor->ago == 0) t_current = temperature->compute_scalar(); if (pstat_flag) nh_v_press(); // compute new T,P after velocities rescaled by nh_v_press() // compute appropriately coupled elements of mvv_current t_current = temperature->compute_scalar(); + tdof = temperature->dof; + if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); else pressure->compute_vector(); couple(); pressure->addstep(update->ntimestep+1); } if (pstat_flag) nh_omega_dot(); // update eta_dot // update eta_press_dot if (tstat_flag) nhc_temp_integrate(); if (pstat_flag && mpchain) nhc_press_integrate(); } /* ---------------------------------------------------------------------- */ void FixNH::initial_integrate_respa(int vflag, int ilevel, int iloop) { // 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 // 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) { compute_temp_target(); 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(); pressure->compute_scalar(); } else { - temperature->compute_vector(); + temperature->compute_vector(); pressure->compute_vector(); } couple(); pressure->addstep(update->ntimestep+1); } if (pstat_flag) { compute_press_target(); nh_omega_dot(); nh_v_press(); } nve_v(); } else nve_v(); // innermost level - also update x only for atoms in group // if barostat, perform 1/2 step remap before and after if (ilevel == 0) { if (pstat_flag) remap(); nve_x(); if (pstat_flag) remap(); } // if barostat, redo KSpace coeffs at outermost level, // since volume has changed if (ilevel == nlevels_respa-1 && kspace_flag && pstat_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]; } } /* ---------------------------------------------------------------------- change box size remap all atoms or dilate group atoms depending on allremap flag if rigid bodies exist, scale rigid body centers-of-mass ------------------------------------------------------------------------- */ void FixNH::remap() { int i; double oldlo,oldhi; double expfac; double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; double *h = domain->h; // omega is not used, except for book-keeping for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i]; // convert pertinent atoms and rigid bodies to lamda coords if (allremap) domain->x2lamda(nlocal); else { for (i = 0; i < nlocal; i++) if (mask[i] & dilate_group_bit) 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 // this operation corresponds to applying the // translate and scale operations // corresponding to the solution of the following ODE: // // h_dot = omega_dot * h // // where h_dot, omega_dot and h are all upper-triangular // 3x3 tensors. In Voigt notation, the elements of the // RHS product tensor are: // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1] // // Ordering of operations preserves time symmetry. double dto2 = dto/2.0; double dto4 = dto/4.0; double dto8 = dto/8.0; // off-diagonal components, first half if (pstyle == TRICLINIC) { if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } if (p_flag[3]) { expfac = exp(dto4*omega_dot[1]); h[3] *= expfac; h[3] += dto2*(omega_dot[3]*h[2]); h[3] *= expfac; } if (p_flag[5]) { expfac = exp(dto4*omega_dot[0]); h[5] *= expfac; h[5] += dto2*(omega_dot[5]*h[1]); h[5] *= expfac; } if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } } // scale diagonal components // scale tilt factors with cell, if set if (p_flag[0]) { oldlo = domain->boxlo[0]; oldhi = domain->boxhi[0]; expfac = exp(dto*omega_dot[0]); domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0]; domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0]; } if (p_flag[1]) { oldlo = domain->boxlo[1]; oldhi = domain->boxhi[1]; expfac = exp(dto*omega_dot[1]); domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1]; domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1]; if (scalexy) h[5] *= expfac; } if (p_flag[2]) { oldlo = domain->boxlo[2]; oldhi = domain->boxhi[2]; expfac = exp(dto*omega_dot[2]); domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2]; domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2]; if (scalexz) h[4] *= expfac; if (scaleyz) h[3] *= expfac; } // off-diagonal components, second half if (pstyle == TRICLINIC) { if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } if (p_flag[3]) { expfac = exp(dto4*omega_dot[1]); h[3] *= expfac; h[3] += dto2*(omega_dot[3]*h[2]); h[3] *= expfac; } if (p_flag[5]) { expfac = exp(dto4*omega_dot[0]); h[5] *= expfac; h[5] += dto2*(omega_dot[5]*h[1]); h[5] *= expfac; } if (p_flag[4]) { expfac = exp(dto8*omega_dot[0]); h[4] *= expfac; h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]); h[4] *= expfac; } } domain->yz = h[3]; domain->xz = h[4]; domain->xy = h[5]; // tilt factor to cell length ratio can not exceed TILTMAX in one step if (domain->yz < -TILTMAX*domain->yprd || domain->yz > TILTMAX*domain->yprd || domain->xz < -TILTMAX*domain->xprd || domain->xz > TILTMAX*domain->xprd || domain->xy < -TILTMAX*domain->xprd || domain->xy > TILTMAX*domain->xprd) error->all(FLERR,"Fix npt/nph has tilted box too far in one step - " "periodic cell is too far from equilibrium state"); 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] & dilate_group_bit) 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 = size_restart_global(); double *list; memory->create(list,nsize,"nh:list"); pack_restart_data(list); if (comm->me == 0) { int size = nsize * sizeof(double); fwrite(&size,sizeof(int),1,fp); fwrite(list,sizeof(double),nsize,fp); } memory->destroy(list); } /* ---------------------------------------------------------------------- calculate the number of data to be packed ------------------------------------------------------------------------- */ int FixNH::size_restart_global() { int nsize = 2; if (tstat_flag) nsize += 1 + 2*mtchain; if (pstat_flag) { nsize += 16 + 2*mpchain; if (deviatoric_flag) nsize += 6; } return nsize; } /* ---------------------------------------------------------------------- pack restart data ------------------------------------------------------------------------- */ int FixNH::pack_restart_data(double *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]; } } return n; } /* ---------------------------------------------------------------------- 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(FLERR,"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(FLERR,"Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) error->all(FLERR, "Fix_modify temperature ID does not compute temperature"); if (temperature->igroup != 0 && comm->me == 0) error->warning(FLERR,"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(FLERR,"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(FLERR,"Illegal fix_modify command"); if (!pstat_flag) error->all(FLERR,"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(FLERR,"Could not find fix_modify pressure ID"); pressure = modify->compute[icompute]; if (pressure->pressflag == 0) error->all(FLERR,"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_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 += ke_target * 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]; } // 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 if (pstat_flag) { 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; } } double volume; double kt = boltz * t_target; 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 ke_target * 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 pdim*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_target(double t_new) { t_target = t_start = t_stop = t_new; } /* ---------------------------------------------------------------------- */ 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; dto = dthalf; // If using respa, then remap is performed in innermost level if (strstr(update->integrate_style,"respa")) dto = 0.5*step_respa[0]; if (pstat_flag) pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain); if (tstat_flag) tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain); } /* ---------------------------------------------------------------------- extract thermostat properties ------------------------------------------------------------------------- */ void *FixNH::extract(const char *str, int &dim) { dim=0; if (strcmp(str,"t_target") == 0) { return &t_target; } else if (strcmp(str,"mtchain") == 0) { return &mtchain; } dim=1; if (strcmp(str,"eta") == 0) { return η } return NULL; } /* ---------------------------------------------------------------------- perform half-step update of chain thermostat variables ------------------------------------------------------------------------- */ void FixNH::nhc_temp_integrate() { int ich; double expfac; double kecurrent = tdof * boltz * t_current; // Update masses, to preserve initial freq, if flag set if (eta_mass_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); } if (eta_mass[0] > 0.0) eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; else eta_dotdot[0] = 0.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] *= tdrag_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] *= tdrag_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; if (eta_mass[0] > 0.0) eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; else eta_dotdot[0] = 0.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,kecurrent; double kt = boltz * t_target; double lkt_press = kt; // Update masses, to preserve initial freq, if flag set if (omega_mass_flag) { 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]); } } if (etap_mass_flag) { 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]; } } 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] *= pdrag_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] *= pdrag_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 factor[3]; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2)); factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2)); factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2)); 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] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); v[i][1] += -dthalf*v[i][2]*omega_dot[3]; } v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; } } } 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] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]); v[i][1] += -dthalf*v[i][2]*omega_dot[3]; } v[i][0] *= factor[0]; v[i][1] *= factor[1]; v[i][2] *= factor[2]; 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 target temperature and kinetic energy -----------------------------------------------------------------------*/ void FixNH::compute_temp_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop-t_start); ke_target = tdof * boltz * t_target; } /* ---------------------------------------------------------------------- compute hydrostatic target pressure -----------------------------------------------------------------------*/ void FixNH::compute_press_target() { double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; 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 -----------------------------------------------------------------------*/ void FixNH::nh_omega_dot() { 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(); mtk_term1 = 0.0; if (mtk_flag) { if (pstyle == ISO) { mtk_term1 = tdof * boltz * t_current; mtk_term1 /= pdim * atom->natoms; } else { double *mvv_current = temperature->vector; for (int i = 0; i < 3; i++) if (p_flag[i]) mtk_term1 += mvv_current[i]; mtk_term1 /= pdim * atom->natoms; } } for (int i = 0; i < 3; i++) if (p_flag[i]) { f_omega = (p_current[i]-p_hydro)*volume / (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i]; if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p); omega_dot[i] += f_omega*dthalf; omega_dot[i] *= pdrag_factor; } mtk_term2 = 0.0; if (mtk_flag) { for (int i = 0; i < 3; i++) if (p_flag[i]) mtk_term2 += omega_dot[i]; mtk_term2 /= pdim * atom->natoms; } 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] *= pdrag_factor; } } } } /* ---------------------------------------------------------------------- if any tilt ratios exceed limits, set flip = 1 and compute new tilt values do not flip in x or y if non-periodic (can tilt but not flip) this is b/c the box length would be changed (dramatically) by flip if yz tilt exceeded, adjust C vector by one B vector if xz tilt exceeded, adjust C vector by one A vector if xy tilt exceeded, adjust B vector by one A vector check yz first since it may change xz, then xz check comes after if any flip occurs, create new box in domain image_flip() adjusts image flags due to box shape change induced by flip remap() puts atoms outside the new box back into the new box perform irregular on atoms in lamda coords to migrate atoms to new procs important that image_flip comes before remap, since remap may change image flags to new values, making eqs in doc of Domain:image_flip incorrect ------------------------------------------------------------------------- */ void FixNH::pre_exchange() { double xprd = domain->xprd; double yprd = domain->yprd; // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP // this avoids immediate re-flipping due to tilt oscillations double xtiltmax = (0.5+DELTAFLIP)*xprd; double ytiltmax = (0.5+DELTAFLIP)*yprd; int flipxy,flipxz,flipyz; flipxy = flipxz = flipyz = 0; if (domain->yperiodic) { if (domain->yz < -ytiltmax) { domain->yz += yprd; domain->xz += domain->xy; flipyz = 1; } else if (domain->yz >= ytiltmax) { domain->yz -= yprd; domain->xz -= domain->xy; flipyz = -1; } } if (domain->xperiodic) { if (domain->xz < -xtiltmax) { domain->xz += xprd; flipxz = 1; } else if (domain->xz >= xtiltmax) { domain->xz -= xprd; flipxz = -1; } if (domain->xy < -xtiltmax) { domain->xy += xprd; flipxy = 1; } else if (domain->xy >= xtiltmax) { domain->xy -= xprd; flipxy = -1; } } int flip = 0; if (flipxy || flipxz || flipyz) flip = 1; if (flip) { domain->set_global_box(); domain->set_local_box(); domain->image_flip(flipxy,flipxz,flipyz); 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]); domain->x2lamda(atom->nlocal); irregular->migrate_atoms(); domain->lamda2x(atom->nlocal); } } /* ---------------------------------------------------------------------- memory usage of Irregular ------------------------------------------------------------------------- */ double FixNH::memory_usage() { double bytes = 0.0; if (irregular) bytes += irregular->memory_usage(); return bytes; } diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp index 5aac69073..a01f1ded3 100644 --- a/src/fix_temp_berendsen.cpp +++ b/src/fix_temp_berendsen.cpp @@ -1,248 +1,250 @@ /* ---------------------------------------------------------------------- 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 "stdlib.h" #include "math.h" #include "fix_temp_berendsen.h" #include "atom.h" #include "force.h" #include "comm.h" #include "input.h" #include "variable.h" #include "group.h" #include "update.h" #include "modify.h" #include "compute.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; enum{NOBIAS,BIAS}; enum{CONSTANT,EQUAL}; /* ---------------------------------------------------------------------- */ FixTempBerendsen::FixTempBerendsen(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (narg != 6) error->all(FLERR,"Illegal fix temp/berendsen command"); // Berendsen thermostat should be applied every step nevery = 1; scalar_flag = 1; global_freq = nevery; extscalar = 1; tstr = NULL; if (strstr(arg[3],"v_") == arg[3]) { int n = strlen(&arg[3][2]) + 1; tstr = new char[n]; strcpy(tstr,&arg[3][2]); tstyle = EQUAL; } else { t_start = force->numeric(FLERR,arg[3]); t_target = t_start; tstyle = CONSTANT; } t_stop = force->numeric(FLERR,arg[4]); t_period = force->numeric(FLERR,arg[5]); // error checks if (t_period <= 0.0) error->all(FLERR,"Fix temp/berendsen period must be > 0.0"); // create a new compute temp style // id = fix-ID + temp, compute group = fix group 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"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; energy = 0; } /* ---------------------------------------------------------------------- */ FixTempBerendsen::~FixTempBerendsen() { delete [] tstr; // delete temperature if fix created it if (tflag) modify->delete_compute(id_temp); delete [] id_temp; } /* ---------------------------------------------------------------------- */ int FixTempBerendsen::setmask() { int mask = 0; mask |= END_OF_STEP; mask |= THERMO_ENERGY; return mask; } /* ---------------------------------------------------------------------- */ void FixTempBerendsen::init() { // check variable if (tstr) { tvar = input->variable->find(tstr); if (tvar < 0) error->all(FLERR,"Variable name for fix temp/berendsen does not exist"); if (input->variable->equalstyle(tvar)) tstyle = EQUAL; else error->all(FLERR,"Variable for fix temp/berendsen is invalid style"); } int icompute = modify->find_compute(id_temp); if (icompute < 0) error->all(FLERR,"Temperature ID for fix temp/berendsen does not exist"); temperature = modify->compute[icompute]; if (temperature->tempbias) which = BIAS; else which = NOBIAS; } /* ---------------------------------------------------------------------- */ void FixTempBerendsen::end_of_step() { double t_current = temperature->compute_scalar(); + double tdof = temperature->dof; + if (t_current == 0.0) error->all(FLERR, "Computed temperature for fix temp/berendsen cannot be 0.0"); double delta = update->ntimestep - update->beginstep; if (delta != 0.0) delta /= update->endstep - update->beginstep; // set current t_target // if variable temp, evaluate variable, wrap with clear/add if (tstyle == CONSTANT) t_target = t_start + delta * (t_stop-t_start); else { modify->clearstep_compute(); t_target = input->variable->compute_equal(tvar); if (t_target < 0.0) error->one(FLERR, "Fix temp/berendsen variable returned negative temperature"); modify->addstep_compute(update->ntimestep + nevery); } // rescale velocities by lamda // for BIAS: // temperature is current, so do not need to re-compute // OK to not test returned v = 0, since lamda is multiplied by v double lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current - 1.0)); - double efactor = 0.5 * force->boltz * temperature->dof; + double efactor = 0.5 * force->boltz * tdof; energy += t_current * (1.0-lamda*lamda) * efactor; double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { v[i][0] *= lamda; v[i][1] *= lamda; v[i][2] *= lamda; } } } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); v[i][0] *= lamda; v[i][1] *= lamda; v[i][2] *= lamda; temperature->restore_bias(i,v[i]); } } } } /* ---------------------------------------------------------------------- */ int FixTempBerendsen::modify_param(int narg, char **arg) { if (strcmp(arg[0],"temp") == 0) { if (narg < 2) error->all(FLERR,"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(FLERR,"Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) error->all(FLERR, "Fix_modify temperature ID does not compute temperature"); if (temperature->igroup != igroup && comm->me == 0) error->warning(FLERR,"Group for fix_modify temp != fix group"); return 2; } return 0; } /* ---------------------------------------------------------------------- */ void FixTempBerendsen::reset_target(double t_new) { t_target = t_start = t_stop = t_new; } /* ---------------------------------------------------------------------- */ double FixTempBerendsen::compute_scalar() { return energy; } /* ---------------------------------------------------------------------- extract thermostat properties ------------------------------------------------------------------------- */ void *FixTempBerendsen::extract(const char *str, int &dim) { dim=0; if (strcmp(str,"t_target") == 0) { return &t_target; } return NULL; } diff --git a/src/modify.cpp b/src/modify.cpp index 3abb15dd4..865ff48e4 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -1,1290 +1,1295 @@ /* ---------------------------------------------------------------------- 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 "stdio.h" #include "string.h" #include "modify.h" #include "style_compute.h" #include "style_fix.h" #include "atom.h" #include "comm.h" #include "fix.h" #include "compute.h" #include "group.h" #include "update.h" #include "domain.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; using namespace FixConst; #define DELTA 4 #define BIG 1.0e20 #define NEXCEPT 5 // change when add to exceptions in add_fix() /* ---------------------------------------------------------------------- */ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) { nfix = maxfix = 0; n_initial_integrate = n_post_integrate = 0; n_pre_exchange = n_pre_neighbor = 0; n_pre_force = n_post_force = 0; n_final_integrate = n_end_of_step = n_thermo_energy = 0; n_initial_integrate_respa = n_post_integrate_respa = 0; n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0; n_min_pre_exchange = n_min_pre_force = n_min_post_force = n_min_energy = 0; fix = NULL; fmask = NULL; list_initial_integrate = list_post_integrate = NULL; list_pre_exchange = list_pre_neighbor = NULL; list_pre_force = list_post_force = NULL; list_final_integrate = list_end_of_step = NULL; list_thermo_energy = NULL; list_initial_integrate_respa = list_post_integrate_respa = NULL; list_pre_force_respa = list_post_force_respa = NULL; list_final_integrate_respa = NULL; list_min_pre_exchange = list_min_pre_neighbor = NULL; list_min_pre_force = list_min_post_force = NULL; list_min_energy = NULL; end_of_step_every = NULL; list_timeflag = NULL; nfix_restart_global = 0; id_restart_global = style_restart_global = state_restart_global = NULL; nfix_restart_peratom = 0; id_restart_peratom = style_restart_peratom = NULL; index_restart_peratom = NULL; ncompute = maxcompute = 0; compute = NULL; // fill map with fixes listed in style_fix.h fix_map = new std::map(); #define FIX_CLASS #define FixStyle(key,Class) \ (*fix_map)[#key] = &fix_creator; #include "style_fix.h" #undef FixStyle #undef FIX_CLASS // fill map with computes listed in style_compute.h compute_map = new std::map(); #define COMPUTE_CLASS #define ComputeStyle(key,Class) \ (*compute_map)[#key] = &compute_creator; #include "style_compute.h" #undef ComputeStyle #undef COMPUTE_CLASS } /* ---------------------------------------------------------------------- */ Modify::~Modify() { // delete all fixes // do it via delete_fix() so callbacks in Atom are also updated correctly while (nfix) delete_fix(fix[0]->id); memory->sfree(fix); memory->destroy(fmask); // delete all computes for (int i = 0; i < ncompute; i++) delete compute[i]; memory->sfree(compute); delete [] list_initial_integrate; delete [] list_post_integrate; delete [] list_pre_exchange; delete [] list_pre_neighbor; delete [] list_pre_force; delete [] list_post_force; delete [] list_final_integrate; delete [] list_end_of_step; delete [] list_thermo_energy; delete [] list_initial_integrate_respa; delete [] list_post_integrate_respa; delete [] list_pre_force_respa; delete [] list_post_force_respa; delete [] list_final_integrate_respa; delete [] list_min_pre_exchange; delete [] list_min_pre_neighbor; delete [] list_min_pre_force; delete [] list_min_post_force; delete [] list_min_energy; delete [] end_of_step_every; delete [] list_timeflag; restart_deallocate(); delete compute_map; delete fix_map; } /* ---------------------------------------------------------------------- initialize all fixes and computes ------------------------------------------------------------------------- */ void Modify::init() { int i,j; // delete storage of restart info since it is not valid after 1st run restart_deallocate(); // create lists of fixes to call at each stage of run list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate); list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate); list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange); list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor); list_init(PRE_FORCE,n_pre_force,list_pre_force); list_init(POST_FORCE,n_post_force,list_post_force); list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate); list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step); list_init_thermo_energy(THERMO_ENERGY,n_thermo_energy,list_thermo_energy); list_init(INITIAL_INTEGRATE_RESPA, n_initial_integrate_respa,list_initial_integrate_respa); list_init(POST_INTEGRATE_RESPA, n_post_integrate_respa,list_post_integrate_respa); list_init(POST_FORCE_RESPA, n_post_force_respa,list_post_force_respa); list_init(PRE_FORCE_RESPA, n_pre_force_respa,list_pre_force_respa); list_init(FINAL_INTEGRATE_RESPA, n_final_integrate_respa,list_final_integrate_respa); list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange); list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor); list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force); list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force); list_init(MIN_ENERGY,n_min_energy,list_min_energy); // init each fix // not sure if now needs to come before compute init // used to b/c temperature computes called fix->dof() in their init, // and fix rigid required its own init before its dof() could be called, // but computes now do their DOF in setup() for (i = 0; i < nfix; i++) fix[i]->init(); // set global flag if any fix has its restart_pbc flag set restart_pbc_any = 0; for (i = 0; i < nfix; i++) if (fix[i]->restart_pbc) restart_pbc_any = 1; // create list of computes that store invocation times list_init_compute(); // init each compute // set invoked_scalar,vector,etc to -1 to force new run to re-compute them // add initial timestep to all computes that store invocation times // since any of them may be invoked by initial thermo // do not clear out invocation times stored within a compute, // b/c some may be holdovers from previous run, like for ave fixes for (i = 0; i < ncompute; i++) { compute[i]->init(); compute[i]->invoked_scalar = -1; compute[i]->invoked_vector = -1; compute[i]->invoked_array = -1; compute[i]->invoked_peratom = -1; compute[i]->invoked_local = -1; } addstep_compute_all(update->ntimestep); // error if any fix or compute is using a dynamic group when not allowed for (i = 0; i < nfix; i++) if (!fix[i]->dynamic_group_allow && group->dynamic[fix[i]->igroup]) { char str[128]; sprintf(str,"Fix %s does not allow use of dynamic group",fix[i]->id); error->all(FLERR,str); } for (i = 0; i < ncompute; i++) if (!compute[i]->dynamic_group_allow && group->dynamic[compute[i]->igroup]) { char str[128]; sprintf(str,"Compute %s does not allow use of dynamic group",fix[i]->id); error->all(FLERR,str); } // warn if any particle is time integrated more than once int nlocal = atom->nlocal; int *mask = atom->mask; int *flag = new int[nlocal]; for (i = 0; i < nlocal; i++) flag[i] = 0; int groupbit; for (i = 0; i < nfix; i++) { if (fix[i]->time_integrate == 0) continue; groupbit = fix[i]->groupbit; for (j = 0; j < nlocal; j++) if (mask[j] & groupbit) flag[j]++; } int check = 0; for (i = 0; i < nlocal; i++) if (flag[i] > 1) check = 1; delete [] flag; int checkall; MPI_Allreduce(&check,&checkall,1,MPI_INT,MPI_SUM,world); if (comm->me == 0 && checkall) error->warning(FLERR, "One or more atoms are time integrated more than once"); } /* ---------------------------------------------------------------------- setup for run, calls setup() of all fixes and computes called from Verlet, RESPA, Min ------------------------------------------------------------------------- */ void Modify::setup(int vflag) { // compute setup needs to come before fix setup - // b/c NH fixes need DOF of temperature computes + // b/c NH fixes need DOF of temperature computes + // fix group setup() is special case since populates a dynamic group + // needs to be done before temperature compute setup + + for (int i = 0; i < nfix; i++) + if (strcmp(fix[i]->style,"GROUP") == 0) fix[i]->setup(vflag); for (int i = 0; i < ncompute; i++) compute[i]->setup(); if (update->whichflag == 1) for (int i = 0; i < nfix; i++) fix[i]->setup(vflag); else if (update->whichflag == 2) for (int i = 0; i < nfix; i++) fix[i]->min_setup(vflag); } /* ---------------------------------------------------------------------- setup pre_exchange call, only for fixes that define pre_exchange called from Verlet, RESPA, Min, and WriteRestart with whichflag = 0 ------------------------------------------------------------------------- */ void Modify::setup_pre_exchange() { if (update->whichflag <= 1) for (int i = 0; i < n_pre_exchange; i++) fix[list_pre_exchange[i]]->setup_pre_exchange(); else if (update->whichflag == 2) for (int i = 0; i < n_min_pre_exchange; i++) fix[list_min_pre_exchange[i]]->min_setup_pre_exchange(); } /* ---------------------------------------------------------------------- setup pre_neighbor call, only for fixes that define pre_neighbor called from Verlet, RESPA ------------------------------------------------------------------------- */ void Modify::setup_pre_neighbor() { if (update->whichflag == 1) for (int i = 0; i < n_pre_neighbor; i++) fix[list_pre_neighbor[i]]->setup_pre_neighbor(); else if (update->whichflag == 2) for (int i = 0; i < n_min_pre_neighbor; i++) fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor(); } /* ---------------------------------------------------------------------- setup pre_force call, only for fixes that define pre_force called from Verlet, RESPA, Min ------------------------------------------------------------------------- */ void Modify::setup_pre_force(int vflag) { if (update->whichflag == 1) for (int i = 0; i < n_pre_force; i++) fix[list_pre_force[i]]->setup_pre_force(vflag); else if (update->whichflag == 2) for (int i = 0; i < n_min_pre_force; i++) fix[list_min_pre_force[i]]->min_setup_pre_force(vflag); } /* ---------------------------------------------------------------------- 1st half of integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::initial_integrate(int vflag) { for (int i = 0; i < n_initial_integrate; i++) fix[list_initial_integrate[i]]->initial_integrate(vflag); } /* ---------------------------------------------------------------------- post_integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_integrate() { for (int i = 0; i < n_post_integrate; i++) fix[list_post_integrate[i]]->post_integrate(); } /* ---------------------------------------------------------------------- pre_exchange call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::pre_exchange() { for (int i = 0; i < n_pre_exchange; i++) fix[list_pre_exchange[i]]->pre_exchange(); } /* ---------------------------------------------------------------------- pre_neighbor call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::pre_neighbor() { for (int i = 0; i < n_pre_neighbor; i++) fix[list_pre_neighbor[i]]->pre_neighbor(); } /* ---------------------------------------------------------------------- pre_force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::pre_force(int vflag) { for (int i = 0; i < n_pre_force; i++) fix[list_pre_force[i]]->pre_force(vflag); } /* ---------------------------------------------------------------------- post_force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_force(int vflag) { for (int i = 0; i < n_post_force; i++) fix[list_post_force[i]]->post_force(vflag); } /* ---------------------------------------------------------------------- 2nd half of integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::final_integrate() { for (int i = 0; i < n_final_integrate; i++) fix[list_final_integrate[i]]->final_integrate(); } /* ---------------------------------------------------------------------- end-of-timestep call, only for relevant fixes only call fix->end_of_step() on timesteps that are multiples of nevery ------------------------------------------------------------------------- */ void Modify::end_of_step() { for (int i = 0; i < n_end_of_step; i++) if (update->ntimestep % end_of_step_every[i] == 0) fix[list_end_of_step[i]]->end_of_step(); } /* ---------------------------------------------------------------------- thermo energy call, only for relevant fixes called by Thermo class compute_scalar() is fix call to return energy ------------------------------------------------------------------------- */ double Modify::thermo_energy() { double energy = 0.0; for (int i = 0; i < n_thermo_energy; i++) energy += fix[list_thermo_energy[i]]->compute_scalar(); return energy; } /* ---------------------------------------------------------------------- post_run call ------------------------------------------------------------------------- */ void Modify::post_run() { for (int i = 0; i < nfix; i++) fix[i]->post_run(); } /* ---------------------------------------------------------------------- setup rRESPA pre_force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::setup_pre_force_respa(int vflag, int ilevel) { for (int i = 0; i < n_pre_force_respa; i++) fix[list_pre_force_respa[i]]->setup_pre_force_respa(vflag,ilevel); } /* ---------------------------------------------------------------------- 1st half of rRESPA integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::initial_integrate_respa(int vflag, int ilevel, int iloop) { for (int i = 0; i < n_initial_integrate_respa; i++) fix[list_initial_integrate_respa[i]]-> initial_integrate_respa(vflag,ilevel,iloop); } /* ---------------------------------------------------------------------- rRESPA post_integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_integrate_respa(int ilevel, int iloop) { for (int i = 0; i < n_post_integrate_respa; i++) fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,iloop); } /* ---------------------------------------------------------------------- rRESPA pre_force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::pre_force_respa(int vflag, int ilevel, int iloop) { for (int i = 0; i < n_pre_force_respa; i++) fix[list_pre_force_respa[i]]->pre_force_respa(vflag,ilevel,iloop); } /* ---------------------------------------------------------------------- rRESPA post_force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::post_force_respa(int vflag, int ilevel, int iloop) { for (int i = 0; i < n_post_force_respa; i++) fix[list_post_force_respa[i]]->post_force_respa(vflag,ilevel,iloop); } /* ---------------------------------------------------------------------- 2nd half of rRESPA integrate call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::final_integrate_respa(int ilevel, int iloop) { for (int i = 0; i < n_final_integrate_respa; i++) fix[list_final_integrate_respa[i]]->final_integrate_respa(ilevel,iloop); } /* ---------------------------------------------------------------------- minimizer pre-exchange call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_pre_exchange() { for (int i = 0; i < n_min_pre_exchange; i++) fix[list_min_pre_exchange[i]]->min_pre_exchange(); } /* ---------------------------------------------------------------------- minimizer pre-neighbor call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_pre_neighbor() { for (int i = 0; i < n_min_pre_neighbor; i++) fix[list_min_pre_neighbor[i]]->min_pre_neighbor(); } /* ---------------------------------------------------------------------- minimizer pre-force call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_pre_force(int vflag) { for (int i = 0; i < n_min_pre_force; i++) fix[list_min_pre_force[i]]->min_pre_force(vflag); } /* ---------------------------------------------------------------------- minimizer force adjustment call, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_post_force(int vflag) { for (int i = 0; i < n_min_post_force; i++) fix[list_min_post_force[i]]->min_post_force(vflag); } /* ---------------------------------------------------------------------- minimizer energy/force evaluation, only for relevant fixes return energy and forces on extra degrees of freedom ------------------------------------------------------------------------- */ double Modify::min_energy(double *fextra) { int ifix,index; index = 0; double eng = 0.0; for (int i = 0; i < n_min_energy; i++) { ifix = list_min_energy[i]; eng += fix[ifix]->min_energy(&fextra[index]); index += fix[ifix]->min_dof(); } return eng; } /* ---------------------------------------------------------------------- store current state of extra minimizer dof, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_store() { for (int i = 0; i < n_min_energy; i++) fix[list_min_energy[i]]->min_store(); } /* ---------------------------------------------------------------------- manage state of extra minimizer dof on a stack, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_clearstore() { for (int i = 0; i < n_min_energy; i++) fix[list_min_energy[i]]->min_clearstore(); } void Modify::min_pushstore() { for (int i = 0; i < n_min_energy; i++) fix[list_min_energy[i]]->min_pushstore(); } void Modify::min_popstore() { for (int i = 0; i < n_min_energy; i++) fix[list_min_energy[i]]->min_popstore(); } /* ---------------------------------------------------------------------- displace extra minimizer dof along vector hextra, only for relevant fixes ------------------------------------------------------------------------- */ void Modify::min_step(double alpha, double *hextra) { int ifix,index; index = 0; for (int i = 0; i < n_min_energy; i++) { ifix = list_min_energy[i]; fix[ifix]->min_step(alpha,&hextra[index]); index += fix[ifix]->min_dof(); } } /* ---------------------------------------------------------------------- compute max allowed step size along vector hextra, only for relevant fixes ------------------------------------------------------------------------- */ double Modify::max_alpha(double *hextra) { int ifix,index; double alpha = BIG; index = 0; for (int i = 0; i < n_min_energy; i++) { ifix = list_min_energy[i]; double alpha_one = fix[ifix]->max_alpha(&hextra[index]); alpha = MIN(alpha,alpha_one); index += fix[ifix]->min_dof(); } return alpha; } /* ---------------------------------------------------------------------- extract extra minimizer dof, only for relevant fixes ------------------------------------------------------------------------- */ int Modify::min_dof() { int ndof = 0; for (int i = 0; i < n_min_energy; i++) ndof += fix[list_min_energy[i]]->min_dof(); return ndof; } /* ---------------------------------------------------------------------- reset minimizer 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 ------------------------------------------------------------------------- */ void Modify::add_fix(int narg, char **arg, int trysuffix) { if (narg < 3) error->all(FLERR,"Illegal fix command"); // cannot define fix before box exists unless style is in exception list // don't like this way of checking for exceptions by adding fixes to list, // but can't think of better way // too late if instantiate fix, then check flag set in fix constructor, // since some fixes access domain settings in their constructor // MUST change NEXCEPT above when add new fix to this list const char *exceptions[NEXCEPT] = {"GPU","OMP","INTEL","property/atom","cmap"}; if (domain->box_exist == 0) { int m; for (m = 0; m < NEXCEPT; m++) if (strcmp(arg[2],exceptions[m]) == 0) break; if (m == NEXCEPT) error->all(FLERR,"Fix command before simulation box is defined"); } // check group ID int igroup = group->find(arg[1]); if (igroup == -1) error->all(FLERR,"Could not find fix group ID"); // if fix ID exists: // set newflag = 0 so create new fix in same location in fix list // error if new style does not match old style // since can't replace it (all when-to-invoke ptrs would be invalid) // warn if new group != old group // delete old fix, but do not call update_callback(), // since will replace this fix and thus other fix locs will not change // set ptr to NULL in case new fix scans list of fixes, // e.g. scan will occur in add_callback() if called by new fix // if fix ID does not exist: // set newflag = 1 so create new fix // extend fix and fmask lists as necessary int ifix,newflag; for (ifix = 0; ifix < nfix; ifix++) if (strcmp(arg[0],fix[ifix]->id) == 0) break; if (ifix < nfix) { newflag = 0; int match = 0; if (strcmp(arg[2],fix[ifix]->style) == 0) match = 1; if (!match && trysuffix && lmp->suffix_enable) { char estyle[256]; if (lmp->suffix) { sprintf(estyle,"%s/%s",arg[2],lmp->suffix); if (strcmp(estyle,fix[ifix]->style) == 0) match = 1; } if (lmp->suffix2) { sprintf(estyle,"%s/%s",arg[2],lmp->suffix2); if (strcmp(estyle,fix[ifix]->style) == 0) match = 1; } } if (!match) error->all(FLERR,"Replacing a fix, but new style != old style"); if (fix[ifix]->igroup != igroup && comm->me == 0) error->warning(FLERR,"Replacing a fix, but new group != old group"); delete fix[ifix]; fix[ifix] = NULL; } else { newflag = 1; if (nfix == maxfix) { maxfix += DELTA; fix = (Fix **) memory->srealloc(fix,maxfix*sizeof(Fix *),"modify:fix"); memory->grow(fmask,maxfix,"modify:fmask"); } } // create the Fix // try first with suffix appended fix[ifix] = NULL; if (trysuffix && lmp->suffix_enable) { if (lmp->suffix) { char estyle[256]; sprintf(estyle,"%s/%s",arg[2],lmp->suffix); if (fix_map->find(estyle) != fix_map->end()) { FixCreator fix_creator = (*fix_map)[estyle]; fix[ifix] = fix_creator(lmp,narg,arg); } } if (fix[ifix] == NULL && lmp->suffix2) { char estyle[256]; sprintf(estyle,"%s/%s",arg[2],lmp->suffix2); if (fix_map->find(estyle) != fix_map->end()) { FixCreator fix_creator = (*fix_map)[estyle]; fix[ifix] = fix_creator(lmp,narg,arg); } } } if (fix[ifix] == NULL && fix_map->find(arg[2]) != fix_map->end()) { FixCreator fix_creator = (*fix_map)[arg[2]]; fix[ifix] = fix_creator(lmp,narg,arg); } if (fix[ifix] == NULL) error->all(FLERR,"Invalid fix style"); // check if Fix is in restart_global list // if yes, pass state info to the Fix so it can reset itself for (int i = 0; i < nfix_restart_global; i++) if (strcmp(id_restart_global[i],fix[ifix]->id) == 0 && strcmp(style_restart_global[i],fix[ifix]->style) == 0) { fix[ifix]->restart(state_restart_global[i]); if (comm->me == 0) { char *str = (char *) ("Resetting global state of Fix %s Style %s " "from restart file info\n"); if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style); if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style); } } // check if Fix is in restart_peratom list // if yes, loop over atoms so they can extract info from atom->extra array for (int i = 0; i < nfix_restart_peratom; i++) if (strcmp(id_restart_peratom[i],fix[ifix]->id) == 0 && strcmp(style_restart_peratom[i],fix[ifix]->style) == 0) { for (int j = 0; j < atom->nlocal; j++) fix[ifix]->unpack_restart(j,index_restart_peratom[i]); fix[ifix]->restart_reset = 1; if (comm->me == 0) { char *str = (char *) ("Resetting per-atom state of Fix %s Style %s " "from restart file info\n"); if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style); if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style); } } // increment nfix (if new) // set fix mask values // post_construct() allows new fix to create other fixes // nfix increment comes first so that recursive call to add_fix within // post_constructor() will see updated nfix if (newflag) nfix++; fmask[ifix] = fix[ifix]->setmask(); fix[ifix]->post_constructor(); } /* ---------------------------------------------------------------------- one instance per fix in style_fix.h ------------------------------------------------------------------------- */ template Fix *Modify::fix_creator(LAMMPS *lmp, int narg, char **arg) { return new T(lmp,narg,arg); } /* ---------------------------------------------------------------------- modify a Fix's parameters ------------------------------------------------------------------------- */ void Modify::modify_fix(int narg, char **arg) { if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); // lookup Fix ID int ifix; for (ifix = 0; ifix < nfix; ifix++) if (strcmp(arg[0],fix[ifix]->id) == 0) break; if (ifix == nfix) error->all(FLERR,"Could not find fix_modify ID"); fix[ifix]->modify_params(narg-1,&arg[1]); } /* ---------------------------------------------------------------------- delete a Fix from list of Fixes Atom class must update indices in its list of callbacks to fixes ------------------------------------------------------------------------- */ void Modify::delete_fix(const char *id) { int ifix = find_fix(id); if (ifix < 0) error->all(FLERR,"Could not find fix ID to delete"); delete fix[ifix]; atom->update_callback(ifix); // move other Fixes and fmask down in list one slot for (int i = ifix+1; i < nfix; i++) fix[i-1] = fix[i]; for (int i = ifix+1; i < nfix; i++) fmask[i-1] = fmask[i]; nfix--; } /* ---------------------------------------------------------------------- find a fix by ID return index of fix or -1 if not found ------------------------------------------------------------------------- */ int Modify::find_fix(const char *id) { int ifix; for (ifix = 0; ifix < nfix; ifix++) if (strcmp(id,fix[ifix]->id) == 0) break; if (ifix == nfix) return -1; return ifix; } /* ---------------------------------------------------------------------- check for fix associated with package name in compiled list return 1 if found else 0 used to determine whether LAMMPS was built with GPU, USER-INTEL, USER-OMP packages, which have their own fixes ------------------------------------------------------------------------- */ int Modify::check_package(const char *package_fix_name) { if (fix_map->find(package_fix_name) == fix_map->end()) return 0; return 1; } /* ---------------------------------------------------------------------- add a new compute ------------------------------------------------------------------------- */ void Modify::add_compute(int narg, char **arg, int trysuffix) { if (narg < 3) error->all(FLERR,"Illegal compute command"); // error check for (int icompute = 0; icompute < ncompute; icompute++) if (strcmp(arg[0],compute[icompute]->id) == 0) error->all(FLERR,"Reuse of compute ID"); // extend Compute list if necessary if (ncompute == maxcompute) { maxcompute += DELTA; compute = (Compute **) memory->srealloc(compute,maxcompute*sizeof(Compute *),"modify:compute"); } // create the Compute // try first with suffix appended compute[ncompute] = NULL; if (trysuffix && lmp->suffix_enable) { if (lmp->suffix) { char estyle[256]; sprintf(estyle,"%s/%s",arg[2],lmp->suffix); if (compute_map->find(estyle) != compute_map->end()) { ComputeCreator compute_creator = (*compute_map)[estyle]; compute[ncompute] = compute_creator(lmp,narg,arg); } } if (compute[ncompute] == NULL && lmp->suffix2) { char estyle[256]; sprintf(estyle,"%s/%s",arg[2],lmp->suffix2); if (compute_map->find(estyle) != compute_map->end()) { ComputeCreator compute_creator = (*compute_map)[estyle]; compute[ncompute] = compute_creator(lmp,narg,arg); } } } if (compute[ncompute] == NULL && compute_map->find(arg[2]) != compute_map->end()) { ComputeCreator compute_creator = (*compute_map)[arg[2]]; compute[ncompute] = compute_creator(lmp,narg,arg); } if (compute[ncompute] == NULL) error->all(FLERR,"Invalid compute style"); ncompute++; } /* ---------------------------------------------------------------------- one instance per compute in style_compute.h ------------------------------------------------------------------------- */ template Compute *Modify::compute_creator(LAMMPS *lmp, int narg, char **arg) { return new T(lmp,narg,arg); } /* ---------------------------------------------------------------------- modify a Compute's parameters ------------------------------------------------------------------------- */ void Modify::modify_compute(int narg, char **arg) { if (narg < 2) error->all(FLERR,"Illegal compute_modify command"); // lookup Compute ID int icompute; for (icompute = 0; icompute < ncompute; icompute++) if (strcmp(arg[0],compute[icompute]->id) == 0) break; if (icompute == ncompute) error->all(FLERR,"Could not find compute_modify ID"); compute[icompute]->modify_params(narg-1,&arg[1]); } /* ---------------------------------------------------------------------- delete a Compute from list of Computes ------------------------------------------------------------------------- */ void Modify::delete_compute(const char *id) { int icompute = find_compute(id); if (icompute < 0) error->all(FLERR,"Could not find compute ID to delete"); delete compute[icompute]; // move other Computes down in list one slot for (int i = icompute+1; i < ncompute; i++) compute[i-1] = compute[i]; ncompute--; } /* ---------------------------------------------------------------------- find a compute by ID return index of compute or -1 if not found ------------------------------------------------------------------------- */ int Modify::find_compute(const char *id) { int icompute; for (icompute = 0; icompute < ncompute; icompute++) if (strcmp(id,compute[icompute]->id) == 0) break; if (icompute == ncompute) return -1; return icompute; } /* ---------------------------------------------------------------------- clear invoked flag of all computes called everywhere that computes are used, before computes are invoked invoked flag used to avoid re-invoking same compute multiple times and to flag computes that store invocation times as having been invoked ------------------------------------------------------------------------- */ void Modify::clearstep_compute() { for (int icompute = 0; icompute < ncompute; icompute++) compute[icompute]->invoked_flag = 0; } /* ---------------------------------------------------------------------- loop over computes that store invocation times if its invoked flag set on this timestep, schedule next invocation called everywhere that computes are used, after computes are invoked ------------------------------------------------------------------------- */ void Modify::addstep_compute(bigint newstep) { for (int icompute = 0; icompute < n_timeflag; icompute++) if (compute[list_timeflag[icompute]]->invoked_flag) compute[list_timeflag[icompute]]->addstep(newstep); } /* ---------------------------------------------------------------------- loop over all computes schedule next invocation for those that store invocation times called when not sure what computes will be needed on newstep do not loop only over n_timeflag, since may not be set yet ------------------------------------------------------------------------- */ void Modify::addstep_compute_all(bigint newstep) { for (int icompute = 0; icompute < ncompute; icompute++) if (compute[icompute]->timeflag) compute[icompute]->addstep(newstep); } /* ---------------------------------------------------------------------- write to restart file for all Fixes with restart info (1) fixes that have global state (2) fixes that store per-atom quantities ------------------------------------------------------------------------- */ void Modify::write_restart(FILE *fp) { int me = comm->me; int count = 0; for (int i = 0; i < nfix; i++) if (fix[i]->restart_global) count++; if (me == 0) fwrite(&count,sizeof(int),1,fp); int n; for (int i = 0; i < nfix; i++) if (fix[i]->restart_global) { if (me == 0) { n = strlen(fix[i]->id) + 1; fwrite(&n,sizeof(int),1,fp); fwrite(fix[i]->id,sizeof(char),n,fp); n = strlen(fix[i]->style) + 1; fwrite(&n,sizeof(int),1,fp); fwrite(fix[i]->style,sizeof(char),n,fp); } fix[i]->write_restart(fp); } count = 0; for (int i = 0; i < nfix; i++) if (fix[i]->restart_peratom) count++; if (me == 0) fwrite(&count,sizeof(int),1,fp); for (int i = 0; i < nfix; i++) if (fix[i]->restart_peratom) { int maxsize_restart = fix[i]->maxsize_restart(); if (me == 0) { n = strlen(fix[i]->id) + 1; fwrite(&n,sizeof(int),1,fp); fwrite(fix[i]->id,sizeof(char),n,fp); n = strlen(fix[i]->style) + 1; fwrite(&n,sizeof(int),1,fp); fwrite(fix[i]->style,sizeof(char),n,fp); fwrite(&maxsize_restart,sizeof(int),1,fp); } } } /* ---------------------------------------------------------------------- read in restart file data on all previously defined Fixes with restart info (1) fixes that have global state (2) fixes that store per-atom quantities return maxsize of extra info that will be stored with any atom ------------------------------------------------------------------------- */ int Modify::read_restart(FILE *fp) { // nfix_restart_global = # of restart entries with global state info int me = comm->me; if (me == 0) fread(&nfix_restart_global,sizeof(int),1,fp); MPI_Bcast(&nfix_restart_global,1,MPI_INT,0,world); // allocate space for each entry if (nfix_restart_global) { id_restart_global = new char*[nfix_restart_global]; style_restart_global = new char*[nfix_restart_global]; state_restart_global = new char*[nfix_restart_global]; } // read each entry and Bcast to all procs // each entry has id string, style string, chunk of state data int n; for (int i = 0; i < nfix_restart_global; i++) { if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); id_restart_global[i] = new char[n]; if (me == 0) fread(id_restart_global[i],sizeof(char),n,fp); MPI_Bcast(id_restart_global[i],n,MPI_CHAR,0,world); if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); style_restart_global[i] = new char[n]; if (me == 0) fread(style_restart_global[i],sizeof(char),n,fp); MPI_Bcast(style_restart_global[i],n,MPI_CHAR,0,world); if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); state_restart_global[i] = new char[n]; if (me == 0) fread(state_restart_global[i],sizeof(char),n,fp); MPI_Bcast(state_restart_global[i],n,MPI_CHAR,0,world); } // nfix_restart_peratom = # of restart entries with peratom info int maxsize = 0; if (me == 0) fread(&nfix_restart_peratom,sizeof(int),1,fp); MPI_Bcast(&nfix_restart_peratom,1,MPI_INT,0,world); // allocate space for each entry if (nfix_restart_peratom) { id_restart_peratom = new char*[nfix_restart_peratom]; style_restart_peratom = new char*[nfix_restart_peratom]; index_restart_peratom = new int[nfix_restart_peratom]; } // read each entry and Bcast to all procs // each entry has id string, style string, maxsize of one atom's data // set index = which set of extra data this fix represents for (int i = 0; i < nfix_restart_peratom; i++) { if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); id_restart_peratom[i] = new char[n]; if (me == 0) fread(id_restart_peratom[i],sizeof(char),n,fp); MPI_Bcast(id_restart_peratom[i],n,MPI_CHAR,0,world); if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); style_restart_peratom[i] = new char[n]; if (me == 0) fread(style_restart_peratom[i],sizeof(char),n,fp); MPI_Bcast(style_restart_peratom[i],n,MPI_CHAR,0,world); if (me == 0) fread(&n,sizeof(int),1,fp); MPI_Bcast(&n,1,MPI_INT,0,world); maxsize += n; index_restart_peratom[i] = i; } return maxsize; } /* ---------------------------------------------------------------------- delete all lists of restart file Fix info ------------------------------------------------------------------------- */ void Modify::restart_deallocate() { if (nfix_restart_global) { for (int i = 0; i < nfix_restart_global; i++) { delete [] id_restart_global[i]; delete [] style_restart_global[i]; delete [] state_restart_global[i]; } delete [] id_restart_global; delete [] style_restart_global; delete [] state_restart_global; } if (nfix_restart_peratom) { for (int i = 0; i < nfix_restart_peratom; i++) { delete [] id_restart_peratom[i]; delete [] style_restart_peratom[i]; } delete [] id_restart_peratom; delete [] style_restart_peratom; delete [] index_restart_peratom; } nfix_restart_global = nfix_restart_peratom = 0; } /* ---------------------------------------------------------------------- create list of fix indices for fixes which match mask ------------------------------------------------------------------------- */ void Modify::list_init(int mask, int &n, int *&list) { delete [] list; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++; list = new int[n]; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask) list[n++] = i; } /* ---------------------------------------------------------------------- create list of fix indices for end_of_step fixes also create end_of_step_every[] ------------------------------------------------------------------------- */ void Modify::list_init_end_of_step(int mask, int &n, int *&list) { delete [] list; delete [] end_of_step_every; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++; list = new int[n]; end_of_step_every = new int[n]; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask) { list[n] = i; end_of_step_every[n++] = fix[i]->nevery; } } /* ---------------------------------------------------------------------- create list of fix indices for thermo energy fixes only added to list if fix has THERMO_ENERGY mask and its thermo_energy flag was set via fix_modify ------------------------------------------------------------------------- */ void Modify::list_init_thermo_energy(int mask, int &n, int *&list) { delete [] list; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask && fix[i]->thermo_energy) n++; list = new int[n]; n = 0; for (int i = 0; i < nfix; i++) if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i; } /* ---------------------------------------------------------------------- create list of compute indices for computes which store invocation times ------------------------------------------------------------------------- */ void Modify::list_init_compute() { delete [] list_timeflag; n_timeflag = 0; for (int i = 0; i < ncompute; i++) if (compute[i]->timeflag) n_timeflag++; list_timeflag = new int[n_timeflag]; n_timeflag = 0; for (int i = 0; i < ncompute; i++) if (compute[i]->timeflag) list_timeflag[n_timeflag++] = i; } /* ---------------------------------------------------------------------- return # of bytes of allocated memory from all fixes ------------------------------------------------------------------------- */ bigint Modify::memory_usage() { bigint bytes = 0; for (int i = 0; i < nfix; i++) bytes += static_cast (fix[i]->memory_usage()); for (int i = 0; i < ncompute; i++) bytes += static_cast (compute[i]->memory_usage()); return bytes; }