git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@831 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/fix_nve_limit.cpp

@@ -0,0 +1,195 @@
+/* ----------------------------------------------------------------------
+   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 "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_nve_limit.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "respa.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixNVELimit::FixNVELimit(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 4) error->all("Illegal fix nve/limit command");
+
+  xlimit = atof(arg[3]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixNVELimit::setmask()
+{
+  int mask = 0;
+  mask |= INITIAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE;
+  mask |= INITIAL_INTEGRATE_RESPA;
+  mask |= FINAL_INTEGRATE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVELimit::init()
+{
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * force->ftm2v;
+  vlimitsq = (xlimit/dtv) * (xlimit/dtv);
+
+  if (strcmp(update->integrate_style,"respa") == 0)
+    step_respa = ((Respa *) update->integrate)->step;
+}
+
+/* ----------------------------------------------------------------------
+   allow for both per-type and per-atom mass
+------------------------------------------------------------------------- */
+
+void FixNVELimit::initial_integrate()
+{
+  double dtfm,vsq,scale;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (mass) {
+    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];
+
+	vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+	if (vsq > vlimitsq) {
+	  scale = sqrt(vlimitsq/vsq);
+	  v[i][0] *= scale;
+	  v[i][1] *= scale;
+	  v[i][2] *= scale;
+	}
+
+	x[i][0] += dtv * v[i][0];
+	x[i][1] += dtv * v[i][1];
+	x[i][2] += dtv * v[i][2];
+      }
+    }
+
+  } else {
+    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];
+
+	vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+	if (vsq > vlimitsq) {
+	  scale = sqrt(vlimitsq/vsq);
+	  v[i][0] *= scale;
+	  v[i][1] *= scale;
+	  v[i][2] *= scale;
+	}
+
+	x[i][0] += dtv * v[i][0];
+	x[i][1] += dtv * v[i][1];
+	x[i][2] += dtv * v[i][2];
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVELimit::final_integrate()
+{
+  double dtfm,vsq,scale;
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (mass) {
+    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];
+
+	vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+	if (vsq > vlimitsq) {
+	  scale = sqrt(vlimitsq/vsq);
+	  v[i][0] *= scale;
+	  v[i][1] *= scale;
+	  v[i][2] *= scale;
+	}
+      }
+    }
+
+  } else {
+    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];
+
+	vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
+	if (vsq > vlimitsq) {
+	  scale = sqrt(vlimitsq/vsq);
+	  v[i][0] *= scale;
+	  v[i][1] *= scale;
+	  v[i][2] *= scale;
+	}
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVELimit::initial_integrate_respa(int ilevel, int flag)
+{
+  if (flag) return;             // only used by NPT,NPH
+
+  dtv = step_respa[ilevel];
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+
+  if (ilevel == 0) initial_integrate();
+  else final_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVELimit::final_integrate_respa(int ilevel)
+{
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+  final_integrate();
+}

src/fix_nve_limit.h

@@ -0,0 +1,40 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_NVE_LIMIT_H
+#define FIX_NVE_LIMIT_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixNVELimit : public Fix {
+ public:
+  FixNVELimit(class LAMMPS *, int, char **);
+  int setmask();
+  void init();
+  void initial_integrate();
+  void final_integrate();
+  void initial_integrate_respa(int, int);
+  void final_integrate_respa(int);
+
+ private:
+  double dtv,dtf;
+  double *step_respa;
+  int mass_require;
+  double xlimit,vlimitsq;
+};
+
+}
+
+#endif

src/style.h

@@ -157,7 +157,9 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_momentum.h"
 #include "fix_nph.h"
 #include "fix_npt.h"
+  //#include "fix_npt_pr.h"
 #include "fix_nve.h"
+#include "fix_nve_limit.h"
 #include "fix_nve_noforce.h"
 #include "fix_nvt.h"
 #include "fix_nvt_sllod.h"
@@ -204,9 +206,11 @@ FixStyle(MINIMIZE,FixMinimize)
 FixStyle(momentum,FixMomentum)
 FixStyle(msd,FixMSD)
 FixStyle(nph,FixNPH)
+  //FixStyle(npt/pr,FixNPTpr)
 FixStyle(npt,FixNPT)
 FixStyle(nve,FixNVE)
-FixStyle(nve,FixNVENoforce)
+FixStyle(nve/limit,FixNVELimit)
+FixStyle(nve/noforce,FixNVENoforce)
 FixStyle(nvt,FixNVT)
 FixStyle(nvt/sllod,FixNVTSlodd)
 FixStyle(orient/fcc,FixOrientFCC)
@@ -329,6 +333,8 @@ RegionStyle(union,RegUnion)
 #include "style_poems.h"
 #include "style_xtc.h"
 
+//#include "style_extra.h"
+
 // user add-ons
 
 #include "style_user.h"