From 3db8850f0992bedcaccb8ea9607bbd6e9361beb8 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 15 Feb 2016 15:47:38 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14598 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/BODY/compute_temp_body.cpp | 394 +++++++++++++++++++++++++++++++++ src/BODY/compute_temp_body.h | 96 ++++++++ src/BODY/fix_nh_body.cpp | 154 +++++++++++++ src/BODY/fix_nh_body.h | 50 +++++ src/BODY/fix_nph_body.cpp | 72 ++++++ src/BODY/fix_nph_body.h | 48 ++++ src/BODY/fix_npt_body.cpp | 72 ++++++ src/BODY/fix_npt_body.h | 48 ++++ src/BODY/fix_nvt_body.cpp | 53 +++++ src/BODY/fix_nvt_body.h | 48 ++++ 10 files changed, 1035 insertions(+) create mode 100755 src/BODY/compute_temp_body.cpp create mode 100755 src/BODY/compute_temp_body.h create mode 100644 src/BODY/fix_nh_body.cpp create mode 100644 src/BODY/fix_nh_body.h create mode 100644 src/BODY/fix_nph_body.cpp create mode 100644 src/BODY/fix_nph_body.h create mode 100755 src/BODY/fix_npt_body.cpp create mode 100755 src/BODY/fix_npt_body.h create mode 100755 src/BODY/fix_nvt_body.cpp create mode 100755 src/BODY/fix_nvt_body.h diff --git a/src/BODY/compute_temp_body.cpp b/src/BODY/compute_temp_body.cpp new file mode 100755 index 0000000000..75ae512e5c --- /dev/null +++ b/src/BODY/compute_temp_body.cpp @@ -0,0 +1,394 @@ +/* ---------------------------------------------------------------------- + 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: Trung Dac Nguyen (ndactrung@gmail.com) + based on ComputeTempAsphere +------------------------------------------------------------------------- */ + +#include +#include +#include "compute_temp_body.h" +#include "math_extra.h" +#include "atom.h" +#include "atom_vec_body.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}; + +/* ---------------------------------------------------------------------- */ + +ComputeTempBody::ComputeTempBody(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg < 3) error->all(FLERR,"Illegal compute temp/body 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/body 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/body 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/body command"); + iarg += 2; + } else error->all(FLERR,"Illegal compute temp/body command"); + } + + vector = new double[6]; + +} + +/* ---------------------------------------------------------------------- */ + +ComputeTempBody::~ComputeTempBody() +{ + delete [] id_bias; + delete [] vector; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempBody::init() +{ + // error check + + avec = (AtomVecBody *) atom->style_match("body"); + if (!avec) + error->all(FLERR,"Compute temp/body requires atom style body"); + + // check that all particles are finite-size, no point particles allowed + + int *body = atom->body; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (body[i] < 0) + error->one(FLERR,"Compute temp/body requires bodies"); + + 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"); + if (strcmp(tbias->style,"temp/region") == 0) tempbias = 2; + else tempbias = 1; + + // init and setup bias compute because + // this compute's setup()->dof_compute() may be called first + + tbias->init(); + tbias->setup(); + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempBody::setup() +{ + dynamic = 0; + if (dynamic_user || group->dynamic[igroup]) dynamic = 1; + dof_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempBody::dof_compute() +{ + 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 + + natoms_temp = 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_temp; + + // additional adjustments to dof + + if (tempbias == 1) { + if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms_temp; + + } 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 ComputeTempBody::compute_scalar() +{ + invoked_scalar = update->ntimestep; + + if (tempbias) { + if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar(); + tbias->remove_bias_all(); + } + + AtomVecBody::Bonus *bonus = avec->bonus; + double **v = atom->v; + double **angmom = atom->angmom; + double *rmass = atom->rmass; + int *body = atom->body; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *inertia,*quat; + double wbody[3]; + double rot[3][3]; + + // sum translational and rotational energy for each particle + + 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 + + inertia = bonus[body[i]].inertia; + quat = bonus[body[i]].quat; + + // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat,rot); + MathExtra::transpose_matvec(rot,angmom[i],wbody); + if (inertia[0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[0]; + if (inertia[1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[1]; + if (inertia[2] == 0.0) wbody[2] = 0.0; + else 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 + + inertia = bonus[body[i]].inertia; + quat = bonus[body[i]].quat; + + // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat,rot); + MathExtra::transpose_matvec(rot,angmom[i],wbody); + if (inertia[0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[0]; + if (inertia[1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[1]; + if (inertia[2] == 0.0) wbody[2] = 0.0; + else 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(); + if (dof < 0.0 && natoms_temp > 0.0) + error->all(FLERR,"Temperature compute degrees of freedom < 0"); + scalar *= tfactor; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTempBody::compute_vector() +{ + int i; + + invoked_vector = update->ntimestep; + + if (tempbias) { + if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector(); + tbias->remove_bias_all(); + } + + AtomVecBody::Bonus *bonus = avec->bonus; + double **v = atom->v; + double **angmom = atom->angmom; + double *rmass = atom->rmass; + int *body = atom->body; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double *inertia,*quat; + double wbody[3],t[6]; + double rot[3][3]; + double massone; + + // sum translational and rotational energy for each particle + + 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 + + inertia = bonus[body[i]].inertia; + quat = bonus[body[i]].quat; + + // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat,rot); + MathExtra::transpose_matvec(rot,angmom[i],wbody); + if (inertia[0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[0]; + if (inertia[1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[1]; + if (inertia[2] == 0.0) wbody[2] = 0.0; + else 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 + + inertia = bonus[body[i]].inertia; + quat = bonus[body[i]].quat; + massone = rmass[i]; + + // wbody = angular velocity in body frame + + MathExtra::quat_to_mat(quat,rot); + MathExtra::transpose_matvec(rot,angmom[i],wbody); + if (inertia[0] == 0.0) wbody[0] = 0.0; + else wbody[0] /= inertia[0]; + if (inertia[1] == 0.0) wbody[1] = 0.0; + else wbody[1] /= inertia[1]; + if (inertia[2] == 0.0) wbody[2] = 0.0; + else 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 ComputeTempBody::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 ComputeTempBody::restore_bias(int i, double *v) +{ + if (tbias) tbias->restore_bias(i,v); +} diff --git a/src/BODY/compute_temp_body.h b/src/BODY/compute_temp_body.h new file mode 100755 index 0000000000..a623fe5898 --- /dev/null +++ b/src/BODY/compute_temp_body.h @@ -0,0 +1,96 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(temp/body,ComputeTempBody) + +#else + +#ifndef LMP_COMPUTE_TEMP_BODY_H +#define LMP_COMPUTE_TEMP_BODY_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeTempBody : public Compute { + public: + ComputeTempBody(class LAMMPS *, int, char **); + ~ComputeTempBody(); + void init(); + void setup(); + double compute_scalar(); + void compute_vector(); + + void remove_bias(int, double *); + void restore_bias(int, double *); + + private: + int mode; + double tfactor; + char *id_bias; + class Compute *tbias; // ptr to additional bias compute + class AtomVecBody *avec; + + void dof_compute(); +}; + +} + +#endif +#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 temp/body requires atom style body + +Self-explanatory. + +E: Compute temp/body requires bodies + +This compute can only be applied to body particles. + +E: Could not find compute ID for temperature bias + +Self-explanatory. + +E: Bias compute does not calculate temperature + +The specified compute must compute temperature. + +E: Bias compute does not calculate a velocity bias + +The specified compute must compute a bias for temperature. + +E: Bias compute group does not match compute group + +The specified compute must operate on the same group as the parent +compute. + +E: Temperature compute degrees of freedom < 0 + +This should not happen if you are calculating the temperature +on a valid set of atoms. + +U: Compute temp/asphere requires atom style ellipsoid + +Self-explanatory. + +*/ diff --git a/src/BODY/fix_nh_body.cpp b/src/BODY/fix_nh_body.cpp new file mode 100644 index 0000000000..a6255a486d --- /dev/null +++ b/src/BODY/fix_nh_body.cpp @@ -0,0 +1,154 @@ +/* ---------------------------------------------------------------------- + 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: Trung Dac Nguyen (ndactrung@gmail.com) + based on FixNHAsphere +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "math_extra.h" +#include "fix_nh_body.h" +#include "atom.h" +#include "atom_vec_body.h" +#include "group.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixNHBody::FixNHBody(LAMMPS *lmp, int narg, char **arg) : + FixNH(lmp, narg, arg) +{ +} + +/* ---------------------------------------------------------------------- */ + +void FixNHBody::init() +{ + avec = (AtomVecBody *) atom->style_match("body"); + if (!avec) + error->all(FLERR, + "Compute nvt/nph/npt body requires atom style body"); + + // check that all particles are finite-size + // no point particles allowed, spherical is OK + + int *body = atom->body; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (body[i] < 0) + error->one(FLERR,"Fix nvt/nph/npt body requires bodies"); + + FixNH::init(); +} + +/* ---------------------------------------------------------------------- + perform half-step update of angular momentum +-----------------------------------------------------------------------*/ + +void FixNHBody::nve_v() +{ + // standard nve_v velocity update + + FixNH::nve_v(); + + double **angmom = atom->angmom; + double **torque = atom->torque; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // update angular momentum by 1/2 step for all particles + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + angmom[i][0] += dtf*torque[i][0]; + angmom[i][1] += dtf*torque[i][1]; + angmom[i][2] += dtf*torque[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + perform full-step update of orientation +-----------------------------------------------------------------------*/ + +void FixNHBody::nve_x() +{ + double omega[3]; + double *quat,*inertia; + + // standard nve_x position update + + FixNH::nve_x(); + + AtomVecBody::Bonus *bonus = avec->bonus; + int *body = atom->body; + double **angmom = atom->angmom; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + + dtq = 0.5 * dtv; + + // update quaternion a full step via Richardson iteration + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + + // compute omega at 1/2 step from angmom at 1/2 step and current q + // update quaternion a full step via Richardson iteration + // returns new normalized quaternion + + inertia = bonus[body[i]].inertia; + quat = bonus[body[i]].quat; + MathExtra::mq_to_omega(angmom[i],quat,inertia,omega); + MathExtra::richardson(quat,angmom[i],omega,inertia,dtq); + } +} + +/* ---------------------------------------------------------------------- + perform half-step temperature scaling of angular momentum +-----------------------------------------------------------------------*/ + +void FixNHBody::nh_v_temp() +{ + // standard nh_v_temp scaling + + FixNH::nh_v_temp(); + + double **angmom = atom->angmom; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + angmom[i][0] *= factor_eta; + angmom[i][1] *= factor_eta; + angmom[i][2] *= factor_eta; + } + } +} diff --git a/src/BODY/fix_nh_body.h b/src/BODY/fix_nh_body.h new file mode 100644 index 0000000000..80fa9d53e0 --- /dev/null +++ b/src/BODY/fix_nh_body.h @@ -0,0 +1,50 @@ +/* -*- 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_FIX_NH_BODY_H +#define LMP_FIX_NH_BODY_H + +#include "fix_nh.h" + +namespace LAMMPS_NS { + +class FixNHBody : public FixNH { + public: + FixNHBody(class LAMMPS *, int, char **); + virtual ~FixNHBody() {} + void init(); + + protected: + double dtq; + class AtomVecBody *avec; + + void nve_v(); + void nve_x(); + void nh_v_temp(); +}; + +} + +#endif + +/* ERROR/WARNING messages: + +E: Compute nvt/nph/npt body requires atom style body + +Self-explanatory. + +E: Fix nvt/nph/npt body requires bodies + +Self-explanatory. + +*/ diff --git a/src/BODY/fix_nph_body.cpp b/src/BODY/fix_nph_body.cpp new file mode 100644 index 0000000000..1808b9ade8 --- /dev/null +++ b/src/BODY/fix_nph_body.cpp @@ -0,0 +1,72 @@ +/* ---------------------------------------------------------------------- + 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: Trung Dac Nguyen (ndactrung@gmail.com) +------------------------------------------------------------------------- */ + +#include +#include "fix_nph_body.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixNPHBody::FixNPHBody(LAMMPS *lmp, int narg, char **arg) : + FixNHBody(lmp, narg, arg) +{ + if (tstat_flag) + error->all(FLERR,"Temperature control can not be used with fix nph/body"); + if (!pstat_flag) + error->all(FLERR,"Pressure control must be used with fix nph/body"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/body"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pflag = 1; +} diff --git a/src/BODY/fix_nph_body.h b/src/BODY/fix_nph_body.h new file mode 100644 index 0000000000..a2e80caf38 --- /dev/null +++ b/src/BODY/fix_nph_body.h @@ -0,0 +1,48 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(nph/body,FixNPHBody) + +#else + +#ifndef LMP_FIX_NPH_BODY_H +#define LMP_FIX_NPH_BODY_H + +#include "fix_nh_body.h" + +namespace LAMMPS_NS { + +class FixNPHBody : public FixNHBody { + public: + FixNPHBody(class LAMMPS *, int, char **); + ~FixNPHBody() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Temperature control can not be used with fix nph/body + +Self-explanatory. + +E: Pressure control must be used with fix nph/body + +Self-explanatory. + +*/ diff --git a/src/BODY/fix_npt_body.cpp b/src/BODY/fix_npt_body.cpp new file mode 100755 index 0000000000..17014ba17c --- /dev/null +++ b/src/BODY/fix_npt_body.cpp @@ -0,0 +1,72 @@ +/* ---------------------------------------------------------------------- + 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: Trung Dac Nguyen (ndactrung@gmail.com) +------------------------------------------------------------------------- */ + +#include +#include "fix_npt_body.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixNPTBody::FixNPTBody(LAMMPS *lmp, int narg, char **arg) : + FixNHBody(lmp, narg, arg) +{ + if (!tstat_flag) + error->all(FLERR,"Temperature control must be used with fix npt/body"); + if (!pstat_flag) + error->all(FLERR,"Pressure control must be used with fix npt/body"); + + // create a new compute temp style + // id = fix-ID + temp + // compute group = all since pressure is always global (group all) + // and thus its KE/temperature contribution should use group all + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "temp/body"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; + + // create a new compute pressure style + // id = fix-ID + press, compute group = all + // pass id_temp as 4th arg to pressure constructor + + n = strlen(id) + 7; + id_press = new char[n]; + strcpy(id_press,id); + strcat(id_press,"_press"); + + newarg = new char*[4]; + newarg[0] = id_press; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "pressure"; + newarg[3] = id_temp; + modify->add_compute(4,newarg); + delete [] newarg; + pflag = 1; +} diff --git a/src/BODY/fix_npt_body.h b/src/BODY/fix_npt_body.h new file mode 100755 index 0000000000..24f6701dad --- /dev/null +++ b/src/BODY/fix_npt_body.h @@ -0,0 +1,48 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(npt/body,FixNPTBody) + +#else + +#ifndef LMP_FIX_NPT_BODY_H +#define LMP_FIX_NPT_BODY_H + +#include "fix_nh_body.h" + +namespace LAMMPS_NS { + +class FixNPTBody : public FixNHBody { + public: + FixNPTBody(class LAMMPS *, int, char **); + ~FixNPTBody() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Temperature control must be used with fix npt/body + +Self-explanatory. + +E: Pressure control must be used with fix npt/body + +Self-explanatory. + +*/ diff --git a/src/BODY/fix_nvt_body.cpp b/src/BODY/fix_nvt_body.cpp new file mode 100755 index 0000000000..d6601f4966 --- /dev/null +++ b/src/BODY/fix_nvt_body.cpp @@ -0,0 +1,53 @@ +/* ---------------------------------------------------------------------- + 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: Trung Dac Nguyen (ndactrung@gmail.com) +------------------------------------------------------------------------- */ + +#include +#include "fix_nvt_body.h" +#include "group.h" +#include "modify.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixNVTBody::FixNVTBody(LAMMPS *lmp, int narg, char **arg) : + FixNHBody(lmp, narg, arg) +{ + if (!tstat_flag) + error->all(FLERR,"Temperature control must be used with fix nvt/body"); + if (pstat_flag) + error->all(FLERR,"Pressure control can not be used with fix nvt/body"); + + // create a new compute temp style + // id = fix-ID + temp + + int n = strlen(id) + 6; + id_temp = new char[n]; + strcpy(id_temp,id); + strcat(id_temp,"_temp"); + + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "temp/body"; + + modify->add_compute(3,newarg); + delete [] newarg; + tflag = 1; +} diff --git a/src/BODY/fix_nvt_body.h b/src/BODY/fix_nvt_body.h new file mode 100755 index 0000000000..26cc422baa --- /dev/null +++ b/src/BODY/fix_nvt_body.h @@ -0,0 +1,48 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(nvt/body,FixNVTBody) + +#else + +#ifndef LMP_FIX_NVT_BODY_H +#define LMP_FIX_NVT_BODY_H + +#include "fix_nh_body.h" + +namespace LAMMPS_NS { + +class FixNVTBody : public FixNHBody { + public: + FixNVTBody(class LAMMPS *, int, char **); + ~FixNVTBody() {} +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Temperature control must be used with fix nvt/body + +Self-explanatory. + +E: Pressure control can not be used with fix nvt/body + +Self-explanatory. + +*/