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

src/ASPHERE/fix_npt_asphere.cpp

@@ -0,0 +1,300 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Mike Brown (SNL)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "stdlib.h"
+#include "math.h"
+#include "fix_npt_asphere.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "force.h"
+#include "compute.h"
+#include "kspace.h"
+#include "update.h"
+#include "domain.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixNPTASphere::FixNPTASphere(LAMMPS *lmp, int narg, char **arg) :
+  FixNPT(lmp, narg, arg)
+{
+  if (atom->quat == NULL || atom->angmom == NULL || atom->torque == NULL)
+    error->all("Fix npt/asphere requires atom attributes "
+	       "quat, angmom, torque");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNPTASphere::init() 
+{
+  FixNPT::init();
+  dtq = 0.5 * update->dt;
+}
+
+/* ----------------------------------------------------------------------
+   1st half of Verlet update 
+------------------------------------------------------------------------- */
+
+void FixNPTASphere::initial_integrate()
+{
+  int i;
+
+  double delta = update->ntimestep - update->beginstep;
+  delta /= update->endstep - update->beginstep;
+
+  // update eta_dot
+
+  t_target = t_start + delta * (t_stop-t_start);
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+  eta += dtv*eta_dot;
+
+  // update omega_dot
+  // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
+
+  double f_omega;
+  double denskt = (atom->natoms*boltz*t_target) / 
+    (domain->xprd*domain->yprd*domain->zprd) * nktv2p;
+
+  for (i = 0; i < 3; i++) {
+    p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
+    f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
+    omega_dot[i] += f_omega*dthalf;
+    omega_dot[i] *= drag_factor;
+    omega[i] += dtv*omega_dot[i];
+    factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
+    dilation[i] = exp(dthalf*omega_dot[i]);
+  }
+  ang_factor = exp(-dthalf*eta_dot);
+
+  // v update only for atoms in NPT group
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **quat = atom->quat;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double dtfm;
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+      v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
+      v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
+      v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
+    }
+  }
+
+  // rescale simulation box and all owned atoms by 1/2 step
+
+  box_dilate(0);
+
+  // x update by full step only for atoms in NPT group
+
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      x[i][0] += dtv * v[i][0];
+      x[i][1] += dtv * v[i][1];
+      x[i][2] += dtv * v[i][2];
+    }
+  }
+
+  // update angular momentum by 1/2 step
+  // update quaternion a full step via Richardson iteration
+  // returns new normalized quaternion
+  for (i = 0; i < nlocal; i++) {    
+    if (mask[i] & groupbit) {
+      angmom[i][0] = angmom[i][0] * ang_factor + dtf * torque[i][0];
+      angmom[i][1] = angmom[i][1] * ang_factor + dtf * torque[i][1];
+      angmom[i][2] = angmom[i][2] * ang_factor + dtf * torque[i][2];
+		
+      double inertia[3];
+      calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
+      richardson(quat[i],angmom[i],inertia);
+    }
+  }
+
+  // rescale simulation box and all owned atoms by 1/2 step
+  // redo KSpace coeffs since volume has changed
+
+  box_dilate(0);
+  if (kspace_flag) force->kspace->setup();
+}
+
+/* ----------------------------------------------------------------------
+   2nd half of Verlet update 
+------------------------------------------------------------------------- */
+
+void FixNPTASphere::final_integrate()
+{
+  int i;
+
+  // v update only for atoms in NPT group
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double dtfm;
+  for (i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+      v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
+      v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
+      v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
+  
+      angmom[i][0] = (angmom[i][0] + dtf * torque[i][0]) * ang_factor;
+      angmom[i][1] = (angmom[i][1] + dtf * torque[i][1]) * ang_factor;
+      angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * ang_factor;
+    }
+  }
+
+  // compute new T,P
+
+  t_current = temperature->compute_scalar();
+  if (press_couple == 0) {
+    if (ptemperature) double ptmp = ptemperature->compute_scalar();
+    double tmp = pressure->compute_scalar();
+  } else {
+    temperature->compute_vector();
+    if (ptemperature) ptemperature->compute_vector();
+    pressure->compute_vector();
+  }
+  couple();
+
+  // update eta_dot
+
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+
+  // update omega_dot
+  // for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
+
+  double f_omega;
+  double denskt = (atom->natoms*boltz*t_target) / 
+    (domain->xprd*domain->yprd*domain->zprd) * nktv2p;
+
+  for (i = 0; i < 3; i++) {
+    f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
+    omega_dot[i] += f_omega*dthalf;
+    omega_dot[i] *= drag_factor;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   Richardson iteration to update quaternion accurately
+------------------------------------------------------------------------- */
+
+void FixNPTASphere::richardson(double *q, double *m, double *moments)
+{
+  // compute omega at 1/2 step from m at 1/2 step and q at 0
+
+  double w[3];
+  omega_from_mq(q,m,moments,w);
+
+  // full update from dq/dt = 1/2 w q
+
+  double wq[4];
+  MathExtra::multiply_vec_quat(w,q,wq);
+
+  double qfull[4];
+  qfull[0] = q[0] + dtq * wq[0];
+  qfull[1] = q[1] + dtq * wq[1];
+  qfull[2] = q[2] + dtq * wq[2];
+  qfull[3] = q[3] + dtq * wq[3];
+  MathExtra::normalize4(qfull);
+
+  // 1st half of update from dq/dt = 1/2 w q
+
+  double qhalf[4];
+  qhalf[0] = q[0] + 0.5*dtq * wq[0];
+  qhalf[1] = q[1] + 0.5*dtq * wq[1];
+  qhalf[2] = q[2] + 0.5*dtq * wq[2];
+  qhalf[3] = q[3] + 0.5*dtq * wq[3];
+  MathExtra::normalize4(qhalf);
+
+  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
+  // recompute wq
+
+  omega_from_mq(qhalf,m,moments,w);
+  MathExtra::multiply_vec_quat(w,qhalf,wq);
+
+  // 2nd half of update from dq/dt = 1/2 w q
+
+  qhalf[0] += 0.5*dtq * wq[0];
+  qhalf[1] += 0.5*dtq * wq[1];
+  qhalf[2] += 0.5*dtq * wq[2];
+  qhalf[3] += 0.5*dtq * wq[3];
+  MathExtra::normalize4(qhalf);
+
+  // corrected Richardson update
+
+  q[0] = 2.0*qhalf[0] - qfull[0];
+  q[1] = 2.0*qhalf[1] - qfull[1];
+  q[2] = 2.0*qhalf[2] - qfull[2];
+  q[3] = 2.0*qhalf[3] - qfull[3];
+  MathExtra::normalize4(q);
+}
+
+/* ----------------------------------------------------------------------
+   compute omega from angular momentum
+   w = omega = angular velocity in space frame
+   wbody = angular velocity in body frame
+   project space-frame angular momentum onto body axes
+     and divide by principal moments
+------------------------------------------------------------------------- */
+
+void FixNPTASphere::omega_from_mq(double *q, double *m, double *inertia,
+				  double *w)
+{
+  double rot[3][3];
+  MathExtra::quat_to_mat(q,rot);
+  
+  double wbody[3];
+  MathExtra::transpose_times_column3(rot,m,wbody);
+  wbody[0] /= inertia[0];
+  wbody[1] /= inertia[1];
+  wbody[2] /= inertia[2];
+  MathExtra::times_column3(rot,wbody,w);
+}
+
+/* ----------------------------------------------------------------------
+   calculate the moment of inertia for an ellipsoid, from mass and radii
+   shape = x,y,z radii in body frame
+------------------------------------------------------------------------- */
+
+void FixNPTASphere::calculate_inertia(double mass, double *shape,
+				      double *inertia)
+{
+  inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0;
+  inertia[1] = mass*(shape[0]*shape[0]+shape[2]*shape[2])/5.0;
+  inertia[2] = mass*(shape[0]*shape[0]+shape[1]*shape[1])/5.0;
+}

src/ASPHERE/fix_npt_asphere.h

@@ -11,29 +11,28 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#ifndef FIX_VOL_RESCALE_H
+#ifndef FIX_NPT_ASPHERE_H
-#define FIX_VOL_RESCALE_H
+#define FIX_NPT_ASPHERE_H
 
-#include "fix.h"
+#include "fix_npt.h"
 
 namespace LAMMPS_NS {
 
-class FixVolRescale : public Fix {
+class FixNPTASphere : public FixNPT {
  public:
-  FixVolRescale(class LAMMPS *, int, char **);
+  FixNPTASphere(class LAMMPS *, int, char **);
-  ~FixVolRescale();
+  ~FixNPTASphere() {}
-  int setmask();
   void init();
-  void end_of_step();
+  void initial_integrate();
+  void final_integrate();
 
  private:
-  int xflag,yflag,zflag;
+  double dtq;
-  double xlo_start,xlo_stop,xhi_start,xhi_stop;
+  double ang_factor;
-  double ylo_start,ylo_stop,yhi_start,yhi_stop;
+
-  double zlo_start,zlo_stop,zhi_start,zhi_stop;
+  void richardson(double *, double *, double *);
-  int kspace_flag;                 // 1 if KSpace invoked, 0 if not
+  void omega_from_mq(double *, double *, double *, double *);
-  int nrigid;                      // number of rigid fixes
+  void calculate_inertia(double mass, double *shape, double *inertia);
-  int *rfix;                       // indices of rigid fixes
 };
 
 }

src/ASPHERE/fix_nvt_asphere.cpp

@@ -0,0 +1,237 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Mike Brown (SNL)
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "stdlib.h"
+#include "math.h"
+#include "fix_nvt_asphere.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "force.h"
+#include "comm.h"
+#include "group.h"
+#include "update.h"
+#include "modify.h"
+#include "compute.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixNVTASphere::FixNVTASphere(LAMMPS *lmp, int narg, char **arg) :
+  FixNVT(lmp, narg, arg)
+{
+  if (atom->quat == NULL || atom->angmom == NULL || atom->torque == NULL)
+    error->all("Fix nvt/asphere requires atom attributes "
+	       "quat, angmom, torque");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTASphere::init()
+{
+  FixNVT::init();
+  dtq = 0.5 * update->dt;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTASphere::initial_integrate()
+{
+  double dtfm;
+
+  double delta = update->ntimestep - update->beginstep;
+  delta /= update->endstep - update->beginstep;
+  t_target = t_start + delta * (t_stop-t_start);
+
+  // update eta_dot
+
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+  eta += dtv*eta_dot;
+  factor = exp(-dthalf*eta_dot);
+
+  // update v and x of only atoms in NVT group
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **quat = atom->quat;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+      v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+      v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+      v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+      x[i][0] += dtv * v[i][0];
+      x[i][1] += dtv * v[i][1];
+      x[i][2] += dtv * v[i][2];
+
+      // update angular momentum by 1/2 step
+      // update quaternion a full step via Richardson iteration
+      // returns new normalized quaternion
+      
+      angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0];
+      angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1];
+      angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2];
+
+      double inertia[3];
+      calculate_inertia(atom->mass[type[i]],atom->shape[type[i]],inertia);
+      richardson(quat[i],angmom[i],inertia);
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTASphere::final_integrate()
+{
+  double dtfm;
+
+  // update v of only atoms in NVT group
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double **angmom = atom->angmom;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]] * factor;
+      v[i][0] = v[i][0]*factor + dtfm*f[i][0];
+      v[i][1] = v[i][1]*factor + dtfm*f[i][1];
+      v[i][2] = v[i][2]*factor + dtfm*f[i][2];
+
+      angmom[i][0] = angmom[i][0] * factor + dtf * torque[i][0];
+      angmom[i][1] = angmom[i][1] * factor + dtf * torque[i][1];
+      angmom[i][2] = angmom[i][2] * factor + dtf * torque[i][2];
+    }
+  }
+
+  // compute current T
+
+  t_current = temperature->compute_scalar();
+
+  // update eta_dot
+
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+}
+
+/* ----------------------------------------------------------------------
+   Richardson iteration to update quaternion accurately
+------------------------------------------------------------------------- */
+
+void FixNVTASphere::richardson(double *q, double *m, double *moments)
+{
+  // compute omega at 1/2 step from m at 1/2 step and q at 0
+
+  double w[3];
+  omega_from_mq(q,m,moments,w);
+
+  // full update from dq/dt = 1/2 w q
+
+  double wq[4];
+  MathExtra::multiply_vec_quat(w,q,wq);
+
+  double qfull[4];
+  qfull[0] = q[0] + dtq * wq[0];
+  qfull[1] = q[1] + dtq * wq[1];
+  qfull[2] = q[2] + dtq * wq[2];
+  qfull[3] = q[3] + dtq * wq[3];
+  MathExtra::normalize4(qfull);
+
+  // 1st half of update from dq/dt = 1/2 w q
+
+  double qhalf[4];
+  qhalf[0] = q[0] + 0.5*dtq * wq[0];
+  qhalf[1] = q[1] + 0.5*dtq * wq[1];
+  qhalf[2] = q[2] + 0.5*dtq * wq[2];
+  qhalf[3] = q[3] + 0.5*dtq * wq[3];
+  MathExtra::normalize4(qhalf);
+
+  // re-compute omega at 1/2 step from m at 1/2 step and q at 1/2 step
+  // recompute wq
+
+  omega_from_mq(qhalf,m,moments,w);
+  MathExtra::multiply_vec_quat(w,qhalf,wq);
+
+  // 2nd half of update from dq/dt = 1/2 w q
+
+  qhalf[0] += 0.5*dtq * wq[0];
+  qhalf[1] += 0.5*dtq * wq[1];
+  qhalf[2] += 0.5*dtq * wq[2];
+  qhalf[3] += 0.5*dtq * wq[3];
+  MathExtra::normalize4(qhalf);
+
+  // corrected Richardson update
+
+  q[0] = 2.0*qhalf[0] - qfull[0];
+  q[1] = 2.0*qhalf[1] - qfull[1];
+  q[2] = 2.0*qhalf[2] - qfull[2];
+  q[3] = 2.0*qhalf[3] - qfull[3];
+  MathExtra::normalize4(q);
+}
+
+/* ----------------------------------------------------------------------
+   compute omega from angular momentum
+   w = omega = angular velocity in space frame
+   wbody = angular velocity in body frame
+   project space-frame angular momentum onto body axes
+     and divide by principal moments
+------------------------------------------------------------------------- */
+
+void FixNVTASphere::omega_from_mq(double *q, double *m, double *inertia,
+				  double *w)
+{
+  double rot[3][3];
+  MathExtra::quat_to_mat(q,rot);
+  
+  double wbody[3];
+  MathExtra::transpose_times_column3(rot,m,wbody);
+  wbody[0] /= inertia[0];
+  wbody[1] /= inertia[1];
+  wbody[2] /= inertia[2];
+  MathExtra::times_column3(rot,wbody,w);
+}
+
+/* ----------------------------------------------------------------------
+   calculate the moment of inertia for an ELLIPSOID, from mass and radii
+   shape = x,y,z radii in body frame
+------------------------------------------------------------------------- */
+
+void FixNVTASphere::calculate_inertia(double mass, double *shape,
+				      double *inertia)
+{
+  inertia[0] = mass*(shape[1]*shape[1]+shape[2]*shape[2])/5.0;
+  inertia[1] = mass*(shape[0]*shape[0]+shape[2]*shape[2])/5.0;
+  inertia[2] = mass*(shape[0]*shape[0]+shape[1]*shape[1])/5.0;
+}

src/ASPHERE/fix_nvt_asphere.h

@@ -0,0 +1,39 @@
+/* ----------------------------------------------------------------------
+   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_NVT_ASPHERE_H
+#define FIX_NVT_ASPHERE_H
+
+#include "fix_nvt.h"
+
+namespace LAMMPS_NS {
+
+class FixNVTASphere : public FixNVT {
+ public:
+  FixNVTASphere(class LAMMPS *, int, char **);
+  ~FixNVTASphere() {}
+  void init();
+  void initial_integrate();
+  void final_integrate();
+
+ private:
+  double dtq;
+
+  void richardson(double *, double *, double *);
+  void omega_from_mq(double *, double *, double *, double *);
+  void calculate_inertia(double mass, double *shape, double *inertia);
+};
+
+}
+
+#endif

src/DIPOLE/Install.csh

@@ -0,0 +1,32 @@
+# Install/unInstall package classes in LAMMPS
+
+if ($1 == 1) then
+
+  cp style_dipole.h ..
+
+  cp atom_vec_dipole.cpp ..
+  cp compute_temp_dipole.cpp ..
+  cp fix_nve_dipole.cpp ..
+  cp pair_dipole_cut.cpp ..
+
+  cp atom_vec_dipole.h ..
+  cp compute_temp_dipole.h ..
+  cp fix_nve_dipole.h ..
+  cp pair_dipole_cut.h ..
+
+else if ($1 == 0) then
+
+  rm ../style_dipole.h
+  touch ../style_dipole.h
+
+  rm ../atom_vec_dipole.cpp
+  rm ../compute_temp_dipole.cpp
+  rm ../fix_nve_dipole.cpp
+  rm ../pair_dipole_cut.cpp
+
+  rm ../atom_vec_dipole.h
+  rm ../compute_temp_dipole.h
+  rm ../fix_nve_dipole.h
+  rm ../pair_dipole_cut.h
+
+endif

src/DIPOLE/atom_vec_dipole.cpp

@@ -0,0 +1,657 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "atom_vec_dipole.h"
+#include "atom.h"
+#include "domain.h"
+#include "modify.h"
+#include "fix.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define DELTA 10000
+
+/* ---------------------------------------------------------------------- */
+
+AtomVecDipole::AtomVecDipole(LAMMPS *lmp, int narg, char **arg) :
+  AtomVec(lmp, narg, arg)
+{
+  mass_type = 1;
+  shape_type = 1;
+  dipole_type = 1;
+  comm_x_only = comm_f_only = 0;
+  size_comm = 6;
+  size_reverse = 6;
+  size_border = 10;
+  size_data_atom = 9;
+  size_data_vel = 7;
+  xcol_data = 4;
+}
+
+/* ----------------------------------------------------------------------
+   grow atom arrays
+   n = 0 grows arrays by DELTA
+   n > 0 allocates arrays to size n 
+------------------------------------------------------------------------- */
+
+void AtomVecDipole::grow(int n)
+{
+  if (n == 0) nmax += DELTA;
+  else nmax = n;
+  atom->nmax = nmax;
+
+  tag = atom->tag = (int *) 
+    memory->srealloc(atom->tag,nmax*sizeof(int),"atom:tag");
+  type = atom->type = (int *)
+    memory->srealloc(atom->type,nmax*sizeof(int),"atom:type");
+  mask = atom->mask = (int *) 
+    memory->srealloc(atom->mask,nmax*sizeof(int),"atom:mask");
+  image = atom->image = (int *) 
+    memory->srealloc(atom->image,nmax*sizeof(int),"atom:image");
+  x = atom->x = memory->grow_2d_double_array(atom->x,nmax,3,"atom:x");
+  v = atom->v = memory->grow_2d_double_array(atom->v,nmax,3,"atom:v");
+  f = atom->f = memory->grow_2d_double_array(atom->f,nmax,3,"atom:f");
+
+  q = atom->q = (double *) 
+    memory->srealloc(atom->q,nmax*sizeof(double),"atom:q");
+  mu = atom->mu = 
+    memory->grow_2d_double_array(atom->mu,nmax,3,"atom:mu");
+  omega = atom->omega = 
+    memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega");
+  torque = atom->torque =
+    memory->grow_2d_double_array(atom->torque,nmax,3,"atom:torque");
+
+  if (atom->nextra_grow)
+    for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
+      modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDipole::copy(int i, int j)
+{
+  tag[j] = tag[i];
+  type[j] = type[i];
+  mask[j] = mask[i];
+  image[j] = image[i];
+  x[j][0] = x[i][0];
+  x[j][1] = x[i][1];
+  x[j][2] = x[i][2];
+  v[j][0] = v[i][0];
+  v[j][1] = v[i][1];
+  v[j][2] = v[i][2];
+
+  q[j] = q[i];
+  mu[j][0] = mu[i][0];
+  mu[j][1] = mu[i][1];
+  mu[j][2] = mu[i][2];
+  omega[j][0] = omega[i][0];
+  omega[j][1] = omega[i][1];
+  omega[j][2] = omega[i][2];
+
+  if (atom->nextra_grow)
+    for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
+      modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_comm(int n, int *list, double *buf,
+			     int pbc_flag, int *pbc)
+{
+  int i,j,m;
+  double dx,dy,dz;
+
+  m = 0;
+  if (pbc_flag == 0) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = x[j][0];
+      buf[m++] = x[j][1];
+      buf[m++] = x[j][2];
+      buf[m++] = mu[j][0];
+      buf[m++] = mu[j][1];
+      buf[m++] = mu[j][2];
+    }
+  } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
+      dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
+      dz = pbc[2]*domain->zprd;
+    }
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
+      buf[m++] = mu[j][0];
+      buf[m++] = mu[j][1];
+      buf[m++] = mu[j][2];
+    }
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_comm_one(int i, double *buf)
+{
+  buf[0] = mu[i][0];
+  buf[1] = mu[i][1];
+  buf[2] = mu[i][2];
+  return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDipole::unpack_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    x[i][0] = buf[m++];
+    x[i][1] = buf[m++];
+    x[i][2] = buf[m++];
+    mu[i][0] = buf[m++];
+    mu[i][1] = buf[m++];
+    mu[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::unpack_comm_one(int i, double *buf)
+{
+  mu[i][0] = buf[0];
+  mu[i][1] = buf[1];
+  mu[i][2] = buf[2];
+  return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_reverse(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    buf[m++] = f[i][0];
+    buf[m++] = f[i][1];
+    buf[m++] = f[i][2];
+    buf[m++] = torque[i][0];
+    buf[m++] = torque[i][1];
+    buf[m++] = torque[i][2];
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_reverse_one(int i, double *buf)
+{
+  buf[0] = torque[i][0];
+  buf[1] = torque[i][1];
+  buf[2] = torque[i][2];
+  return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDipole::unpack_reverse(int n, int *list, double *buf)
+{
+  int i,j,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    f[j][0] += buf[m++];
+    f[j][1] += buf[m++];
+    f[j][2] += buf[m++];
+    torque[j][0] += buf[m++];
+    torque[j][1] += buf[m++];
+    torque[j][2] += buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::unpack_reverse_one(int i, double *buf)
+{
+  torque[i][0] = buf[0];
+  torque[i][1] = buf[1];
+  torque[i][2] = buf[2];
+  return 3;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_border(int n, int *list, double *buf,
+			       int pbc_flag, int *pbc)
+{
+  int i,j,m;
+  double dx,dy,dz;
+
+  m = 0;
+  if (pbc_flag == 0) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = x[j][0];
+      buf[m++] = x[j][1];
+      buf[m++] = x[j][2];
+      buf[m++] = tag[j];
+      buf[m++] = type[j];
+      buf[m++] = mask[j];
+      buf[m++] = q[j];
+      buf[m++] = mu[j][0];
+      buf[m++] = mu[j][1];
+      buf[m++] = mu[j][2];
+    }
+  } else {
+    if (domain->triclinic == 0) {
+      dx = pbc[0]*domain->xprd;
+      dy = pbc[1]*domain->yprd;
+      dz = pbc[2]*domain->zprd;
+    } else {
+      dx = pbc[0];
+      dy = pbc[1];
+      dz = pbc[2];
+    }
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = x[j][0] + dx;
+      buf[m++] = x[j][1] + dy;
+      buf[m++] = x[j][2] + dz;
+      buf[m++] = tag[j];
+      buf[m++] = type[j];
+      buf[m++] = mask[j];
+      buf[m++] = q[j];
+      buf[m++] = mu[j][0];
+      buf[m++] = mu[j][1];
+      buf[m++] = mu[j][2];
+    }
+  }
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_border_one(int i, double *buf)
+{
+  buf[0] = q[i];
+  buf[1] = mu[i][0];
+  buf[2] = mu[i][1];
+  buf[3] = mu[i][2];
+  return 4;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void AtomVecDipole::unpack_border(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    if (i == nmax) grow(0);
+    x[i][0] = buf[m++];
+    x[i][1] = buf[m++];
+    x[i][2] = buf[m++];
+    tag[i] = static_cast<int> (buf[m++]);
+    type[i] = static_cast<int> (buf[m++]);
+    mask[i] = static_cast<int> (buf[m++]);
+    q[i] = buf[m++];
+    mu[i][0] = buf[m++];
+    mu[i][1] = buf[m++];
+    mu[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::unpack_border_one(int i, double *buf)
+{
+  q[i] = buf[0];
+  mu[i][0] = buf[1];
+  mu[i][1] = buf[2];
+  mu[i][2] = buf[3];
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   pack all atom quantities for shipping to another proc
+   xyz must be 1st 3 values, so that comm::exchange can test on them
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_exchange(int i, double *buf)
+{
+  int m = 1;
+  buf[m++] = x[i][0];
+  buf[m++] = x[i][1];
+  buf[m++] = x[i][2];
+  buf[m++] = v[i][0];
+  buf[m++] = v[i][1];
+  buf[m++] = v[i][2];
+  buf[m++] = tag[i];
+  buf[m++] = type[i];
+  buf[m++] = mask[i];
+  buf[m++] = image[i];
+
+  buf[m++] = q[i];
+  buf[m++] = mu[i][0];
+  buf[m++] = mu[i][1];
+  buf[m++] = mu[i][2];
+  buf[m++] = omega[i][0];
+  buf[m++] = omega[i][1];
+  buf[m++] = omega[i][2];
+
+  if (atom->nextra_grow)
+    for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
+      m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
+
+  buf[0] = m;
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int AtomVecDipole::unpack_exchange(double *buf)
+{
+  int nlocal = atom->nlocal;
+  if (nlocal == nmax) grow(0);
+
+  int m = 1;
+  x[nlocal][0] = buf[m++];
+  x[nlocal][1] = buf[m++];
+  x[nlocal][2] = buf[m++];
+  v[nlocal][0] = buf[m++];
+  v[nlocal][1] = buf[m++];
+  v[nlocal][2] = buf[m++];
+  tag[nlocal] = static_cast<int> (buf[m++]);
+  type[nlocal] = static_cast<int> (buf[m++]);
+  mask[nlocal] = static_cast<int> (buf[m++]);
+  image[nlocal] = static_cast<int> (buf[m++]);
+
+  q[nlocal] = buf[m++];
+  mu[nlocal][0] = buf[m++];
+  mu[nlocal][1] = buf[m++];
+  mu[nlocal][2] = buf[m++];
+  omega[nlocal][0] = buf[m++];
+  omega[nlocal][1] = buf[m++];
+  omega[nlocal][2] = buf[m++];
+
+  if (atom->nextra_grow)
+    for (int iextra = 0; iextra < atom->nextra_grow; iextra++) 
+      m += modify->fix[atom->extra_grow[iextra]]->
+	unpack_exchange(nlocal,&buf[m]);
+
+  atom->nlocal++;
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   size of restart data for all atoms owned by this proc
+   include extra data stored by fixes
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::size_restart()
+{
+  int i;
+
+  int nlocal = atom->nlocal;
+  int n = 18 * nlocal;
+
+  if (atom->nextra_restart)
+    for (int iextra = 0; iextra < atom->nextra_restart; iextra++) 
+      for (i = 0; i < nlocal; i++)
+	n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
+
+  return n;
+}
+
+/* ----------------------------------------------------------------------
+   pack atom I's data for restart file including extra quantities
+   xyz must be 1st 3 values, so that read_restart can test on them
+   molecular types may be negative, but write as positive   
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::pack_restart(int i, double *buf)
+{
+  int m = 1;
+  buf[m++] = x[i][0];
+  buf[m++] = x[i][1];
+  buf[m++] = x[i][2];
+  buf[m++] = tag[i];
+  buf[m++] = type[i];
+  buf[m++] = mask[i];
+  buf[m++] = image[i];
+  buf[m++] = v[i][0];
+  buf[m++] = v[i][1];
+  buf[m++] = v[i][2];
+
+  buf[m++] = q[i];
+  buf[m++] = mu[i][0];
+  buf[m++] = mu[i][1];
+  buf[m++] = mu[i][2];
+  buf[m++] = omega[i][0];
+  buf[m++] = omega[i][1];
+  buf[m++] = omega[i][2];
+
+  if (atom->nextra_restart)
+    for (int iextra = 0; iextra < atom->nextra_restart; iextra++) 
+      m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
+
+  buf[0] = m;
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   unpack data for one atom from restart file including extra quantities
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::unpack_restart(double *buf)
+{
+  int nlocal = atom->nlocal;
+  if (nlocal == nmax) {
+    grow(0);
+    if (atom->nextra_store)
+      atom->extra = memory->grow_2d_double_array(atom->extra,nmax,
+						 atom->nextra_store,
+						 "atom:extra");
+  }
+
+  int m = 1;
+  x[nlocal][0] = buf[m++];
+  x[nlocal][1] = buf[m++];
+  x[nlocal][2] = buf[m++];
+  tag[nlocal] = static_cast<int> (buf[m++]);
+  type[nlocal] = static_cast<int> (buf[m++]);
+  mask[nlocal] = static_cast<int> (buf[m++]);
+  image[nlocal] = static_cast<int> (buf[m++]);
+  v[nlocal][0] = buf[m++];
+  v[nlocal][1] = buf[m++];
+  v[nlocal][2] = buf[m++];
+
+  q[nlocal] = buf[m++];
+  mu[nlocal][0] = buf[m++];
+  mu[nlocal][1] = buf[m++];
+  mu[nlocal][2] = buf[m++];
+  omega[nlocal][0] = buf[m++];
+  omega[nlocal][1] = buf[m++];
+  omega[nlocal][2] = buf[m++];
+
+  double **extra = atom->extra;
+  if (atom->nextra_store) {
+    int size = static_cast<int> (buf[0]) - m;
+    for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
+  }
+
+  atom->nlocal++;
+  return m;
+}
+
+/* ----------------------------------------------------------------------
+   create one atom of itype at coord
+   set other values to defaults
+------------------------------------------------------------------------- */
+
+void AtomVecDipole::create_atom(int itype, double *coord)
+{
+  int nlocal = atom->nlocal;
+  if (nlocal == nmax) grow(0);
+
+  tag[nlocal] = 0;
+  type[nlocal] = itype;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
+  mask[nlocal] = 1;
+  image[nlocal] = (512 << 20) | (512 << 10) | 512;
+  v[nlocal][0] = 0.0;
+  v[nlocal][1] = 0.0;
+  v[nlocal][2] = 0.0;
+
+  q[nlocal] = 0.0;
+  mu[nlocal][0] = 0.0;
+  mu[nlocal][1] = 0.0;
+  mu[nlocal][2] = 0.0;
+  omega[nlocal][0] = 0.0;
+  omega[nlocal][1] = 0.0;
+  omega[nlocal][2] = 0.0;
+
+  atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+   unpack one line from Atoms section of data file
+   initialize other atom quantities
+------------------------------------------------------------------------- */
+
+void AtomVecDipole::data_atom(double *coord, int imagetmp, char **values)
+{
+  int nlocal = atom->nlocal;
+  if (nlocal == nmax) grow(0);
+
+  tag[nlocal] = atoi(values[0]);
+  if (tag[nlocal] <= 0)
+    error->one("Invalid atom ID in Atoms section of data file");
+
+  type[nlocal] = atoi(values[1]);
+  if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
+    error->one("Invalid atom type in Atoms section of data file");
+
+  q[nlocal] = atof(values[2]);
+
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
+
+  mu[nlocal][0] = atof(values[6]);
+  mu[nlocal][1] = atof(values[7]);
+  mu[nlocal][2] = atof(values[8]);
+
+  image[nlocal] = imagetmp;
+
+  mask[nlocal] = 1;
+  v[nlocal][0] = 0.0;
+  v[nlocal][1] = 0.0;
+  v[nlocal][2] = 0.0;
+  omega[nlocal][0] = 0.0;
+  omega[nlocal][1] = 0.0;
+  omega[nlocal][2] = 0.0;
+
+  atom->nlocal++;
+}
+
+/* ----------------------------------------------------------------------
+   unpack hybrid quantities from one line in Atoms section of data file
+   initialize other atom quantities for this sub-style
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::data_atom_hybrid(int nlocal, char **values)
+{
+  q[nlocal] = atof(values[0]);
+  mu[nlocal][0] = atof(values[1]);
+  mu[nlocal][1] = atof(values[2]);
+  mu[nlocal][2] = atof(values[3]);
+
+  v[nlocal][0] = 0.0;
+  v[nlocal][1] = 0.0;
+  v[nlocal][2] = 0.0;
+  omega[nlocal][0] = 0.0;
+  omega[nlocal][1] = 0.0;
+  omega[nlocal][2] = 0.0;
+
+  return 4;
+}
+
+/* ----------------------------------------------------------------------
+   unpack one line from Velocities section of data file
+------------------------------------------------------------------------- */
+
+void AtomVecDipole::data_vel(int m, char **values)
+{
+  v[m][0] = atof(values[0]);
+  v[m][1] = atof(values[1]);
+  v[m][2] = atof(values[2]);
+  omega[m][0] = atof(values[3]);
+  omega[m][1] = atof(values[4]);
+  omega[m][2] = atof(values[5]);
+}
+
+/* ----------------------------------------------------------------------
+   unpack hybrid quantities from one line in Velocities section of data file
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::data_vel_hybrid(int m, char **values)
+{
+  omega[m][0] = atof(values[0]);
+  omega[m][1] = atof(values[1]);
+  omega[m][2] = atof(values[2]);
+  return 3;
+}
+
+/* ----------------------------------------------------------------------
+   return # of bytes of allocated memory 
+------------------------------------------------------------------------- */
+
+int AtomVecDipole::memory_usage()
+{
+  int bytes = 0;
+
+  if (atom->memcheck("tag")) bytes += nmax * sizeof(int);
+  if (atom->memcheck("type")) bytes += nmax * sizeof(int);
+  if (atom->memcheck("mask")) bytes += nmax * sizeof(int);
+  if (atom->memcheck("image")) bytes += nmax * sizeof(int);
+  if (atom->memcheck("x")) bytes += nmax*3 * sizeof(double);
+  if (atom->memcheck("v")) bytes += nmax*3 * sizeof(double);
+  if (atom->memcheck("f")) bytes += nmax*3 * sizeof(double);
+
+  if (atom->memcheck("q")) bytes += nmax * sizeof(double);
+  if (atom->memcheck("mu")) bytes += nmax*3 * sizeof(double);
+  if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double);
+  if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double);
+
+  return bytes;
+}

src/DIPOLE/atom_vec_dipole.h

@@ -0,0 +1,58 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef ATOM_VEC_DIPOLE_H
+#define ATOM_VEC_DIPOLE_H
+
+#include "atom_vec.h"
+
+namespace LAMMPS_NS {
+
+class AtomVecDipole : public AtomVec {
+ public:
+  AtomVecDipole(class LAMMPS *, int, char **);
+  void grow(int);
+  void copy(int, int);
+  int pack_comm(int, int *, double *, int, int *);
+  int pack_comm_one(int, double *);
+  void unpack_comm(int, int, double *);
+  int unpack_comm_one(int, double *);
+  int pack_reverse(int, int, double *);
+  int pack_reverse_one(int, double *);
+  void unpack_reverse(int, int *, double *);
+  int unpack_reverse_one(int, double *);
+  int pack_border(int, int *, double *, int, int *);
+  int pack_border_one(int, double *);
+  void unpack_border(int, int, double *);
+  int unpack_border_one(int, double *);
+  int pack_exchange(int, double *);
+  int unpack_exchange(double *);
+  int size_restart();
+  int pack_restart(int, double *);
+  int unpack_restart(double *);
+  void create_atom(int, double *);
+  void data_atom(double *, int, char **);
+  int data_atom_hybrid(int, char **);
+  void data_vel(int, char **);
+  int data_vel_hybrid(int, char **);
+  int memory_usage();
+
+ private:
+  int *tag,*type,*mask,*image;
+  double **x,**v,**f;
+  double *q,**mu,**omega,**torque;
+};
+
+}
+
+#endif

src/DIPOLE/compute_temp_dipole.cpp

@@ -0,0 +1,144 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "compute_temp_dipole.h"
+#include "atom.h"
+#include "force.h"
+#include "group.h"
+#include "modify.h"
+#include "fix.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+// moment of inertia for a sphere
+
+#define INERTIA 0.4
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTempDipole::ComputeTempDipole(LAMMPS *lmp, int narg, char **arg) : 
+  Compute(lmp, narg, arg)
+{
+  if (narg != 3) error->all("Illegal compute temp/dipole command");
+
+  if (atom->omega == NULL || atom->shape == NULL)
+    error->all("Compute temp/dipole requires atom attributes omega, shape");
+
+  scalar_flag = vector_flag = 1;
+  size_vector = 6;
+  extensive = 0;
+  tempflag = 1;
+
+  vector = new double[6];
+  inertia = new double[atom->ntypes + 1];
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTempDipole::~ComputeTempDipole()
+{
+  delete [] vector;
+  delete [] inertia;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDipole::init()
+{
+  fix_dof = 0;
+  for (int i = 0; i < modify->nfix; i++)
+    fix_dof += modify->fix[i]->dof(igroup);
+  recount();
+
+  // moment of inertia for each particle type
+
+  double *mass = atom->mass;
+  double **shape = atom->shape;
+
+  for (int i = 1; i <= atom->ntypes; i++)
+    inertia[i] = INERTIA * mass[i] * 0.25*shape[i][0]*shape[i][0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDipole::recount()
+{
+  double natoms = group->count(igroup);
+  dof = 2.0 * force->dimension * natoms;
+  dof -= extra_dof + fix_dof;
+  if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz);
+  else tfactor = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double ComputeTempDipole::compute_scalar()
+{
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double **omega = atom->omega;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // rotational and translational kinetic energy
+
+  double t = 0.0;
+  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]] 
+	+ (omega[i][0] * omega[i][0] + omega[i][1] * omega[i][1] +
+	   omega[i][2] * omega[i][2]) * inertia[type[i]];
+
+  MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+  if (dynamic) recount();
+  scalar *= tfactor;
+  return scalar;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDipole::compute_vector()
+{
+  int i;
+
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double **omega = atom->omega;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double rmass,imass,t[6];
+  for (i = 0; i < 6; i++) t[i] = 0.0;
+
+  // rotational and translational kinetic energy
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      rmass = mass[type[i]];
+      imass = inertia[type[i]];
+      t[0] += rmass*v[i][0]*v[i][0] + imass*omega[i][0]*omega[i][0];
+      t[1] += rmass*v[i][1]*v[i][1] + imass*omega[i][1]*omega[i][1];
+      t[2] += rmass*v[i][2]*v[i][2] + imass*omega[i][2]*omega[i][2];
+      t[3] += rmass*v[i][0]*v[i][1] + imass*omega[i][0]*omega[i][1];
+      t[4] += rmass*v[i][0]*v[i][2] + imass*omega[i][0]*omega[i][2];
+      t[5] += rmass*v[i][1]*v[i][2] + imass*omega[i][1]*omega[i][2];
+    }
+
+  MPI_Allreduce(&t,&vector,6,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < 6; i++) vector[i] *= force->mvv2e;
+}

src/DIPOLE/compute_temp_dipole.h

@@ -0,0 +1,39 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_TEMP_DIPOLE_H
+#define COMPUTE_TEMP_DIPOLE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeTempDipole : public Compute {
+ public:
+  ComputeTempDipole(class LAMMPS *, int, char **);
+  ~ComputeTempDipole();
+  void init();
+  double compute_scalar();
+  void compute_vector();
+
+ private:
+  int fix_dof;
+  double tfactor;
+  double *inertia;
+
+  void recount();
+};
+
+}
+
+#endif

src/DIPOLE/fix_nve_dipole.cpp

@@ -0,0 +1,189 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "fix_nve_dipole.h"
+#include "atom.h"
+#include "force.h"
+#include "update.h"
+#include "respa.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+// moment of inertia for a sphere
+
+#define INERTIA 0.4
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEDipole::FixNVEDipole(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 3) error->all("Illegal fix nve/dipole command");
+  if (atom->mu == NULL || atom->omega == NULL || 
+      atom->torque == NULL || atom->shape == NULL)
+    error->all("Fix nve/dipole requires atom attributes "
+	       "mu, omega, torque, shape");
+  inertia = new double[atom->ntypes+1];
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixNVEDipole::~FixNVEDipole()
+{
+  delete [] inertia;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixNVEDipole::setmask()
+{
+  int mask = 0;
+  mask |= INITIAL_INTEGRATE;
+  mask |= FINAL_INTEGRATE;
+  mask |= INITIAL_INTEGRATE_RESPA;
+  mask |= FINAL_INTEGRATE_RESPA;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDipole::init()
+{
+  dtv = update->dt;
+  dtf = 0.5 * update->dt * force->ftm2v;
+
+  if (strcmp(update->integrate_style,"respa") == 0)
+    step_respa = ((Respa *) update->integrate)->step;
+
+  // moment of inertia for each particle type
+
+  double *mass = atom->mass;
+  double **shape = atom->shape;
+
+  for (int i = 1; i <= atom->ntypes; i++)
+    inertia[i] = INERTIA * mass[i] * 0.25*shape[i][0]*shape[i][0];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDipole::initial_integrate()
+{
+  double dtfm,msq,scale;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double **mu = atom->mu;
+  double **omega = atom->omega;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  double *dipole = atom->dipole;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double g[3],h[3];
+
+  // update v,x for all particles
+  // update omega,mu for all dipoles
+  // d_omega/dt = torque / inertia
+  // d_mu/dt = omega cross mu
+  // renormalize mu to dipole length
+
+  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];
+      x[i][0] += dtv * v[i][0];
+      x[i][1] += dtv * v[i][1];
+      x[i][2] += dtv * v[i][2];
+      if (dipole[type[i]] > 0.0) {
+	dtfm = dtf / inertia[type[i]];
+	omega[i][0] += dtfm * torque[i][0];
+	omega[i][1] += dtfm * torque[i][1];
+	omega[i][2] += dtfm * torque[i][2];
+
+	g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2] - omega[i][2]*mu[i][1]);
+	g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0] - omega[i][0]*mu[i][2]);
+	g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1] - omega[i][1]*mu[i][0]);
+
+	msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
+	scale = dipole[type[i]]/sqrt(msq);
+	mu[i][0] = g[0]*scale;
+	mu[i][1] = g[1]*scale;
+	mu[i][2] = g[2]*scale;
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDipole::final_integrate()
+{
+  double dtfm;
+
+  double **v = atom->v;
+  double **f = atom->f;
+  double **omega = atom->omega;
+  double **torque = atom->torque;
+  double *mass = atom->mass;
+  double *dipole = atom->dipole;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // update v for all particles
+  // update omega for all dipoles
+
+  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];
+      if (dipole[type[i]] > 0.0) {
+	dtfm = dtf / inertia[type[i]];
+	omega[i][0] += dtfm * torque[i][0];
+	omega[i][1] += dtfm * torque[i][1];
+	omega[i][2] += dtfm * torque[i][2];
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVEDipole::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 FixNVEDipole::final_integrate_respa(int ilevel)
+{
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+  final_integrate();
+}

src/DIPOLE/fix_nve_dipole.h

@@ -0,0 +1,40 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_NVE_DIPOLE_H
+#define FIX_NVE_DIPOLE_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixNVEDipole : public Fix {
+ public:
+  FixNVEDipole(class LAMMPS *, int, char **);
+  ~FixNVEDipole();
+  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;
+  double *inertia;
+};
+
+}
+
+#endif

src/DIPOLE/pair_dipole_cut.cpp

@@ -0,0 +1,625 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdlib.h"
+#include "pair_dipole_cut.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "comm.h"
+#include "force.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/* ---------------------------------------------------------------------- */
+
+PairDipoleCut::PairDipoleCut(LAMMPS *lmp) : Pair(lmp) {}
+
+/* ---------------------------------------------------------------------- */
+
+PairDipoleCut::~PairDipoleCut()
+{
+  if (allocated) {
+    memory->destroy_2d_int_array(setflag);
+    memory->destroy_2d_double_array(cutsq);
+
+    memory->destroy_2d_double_array(cut_lj);
+    memory->destroy_2d_double_array(cut_ljsq);
+    memory->destroy_2d_double_array(cut_coul);
+    memory->destroy_2d_double_array(cut_coulsq);
+    memory->destroy_2d_double_array(epsilon);
+    memory->destroy_2d_double_array(sigma);
+    memory->destroy_2d_double_array(lj1);
+    memory->destroy_2d_double_array(lj2);
+    memory->destroy_2d_double_array(lj3);
+    memory->destroy_2d_double_array(lj4);
+    memory->destroy_2d_double_array(offset);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairDipoleCut::compute(int eflag, int vflag)
+{
+  int i,j,k,numneigh,itype,jtype;
+  double qtmp,xtmp,ytmp,ztmp,delx,dely,delz;
+  double rsq,rinv,r2inv,r6inv,r3inv,r5inv,r7inv;
+  double forcecoulx,forcecouly,forcecoulz,fforce,crossx,crossy,crossz;
+  double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;
+  double fq,fx,fy,fz;
+  double pdotp,pidotr,pjdotr,pre1,pre2,pre3,pre4;
+  double forcelj,factor_coul,factor_lj;
+  double factor,phicoul,philj;
+  int *neighs;
+  double **f;
+
+  eng_vdwl = eng_coul = 0.0;
+  if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0;
+
+  if (vflag == 2) f = update->f_pair;
+  else f = atom->f;
+  double **x = atom->x;
+  double *q = atom->q;
+  double **mu = atom->mu;
+  double **torque = atom->torque;
+  int *type = atom->type;
+  double *dipole = atom->dipole;
+  int nlocal = atom->nlocal;
+  int nall = atom->nlocal + atom->nghost;
+  double *special_coul = force->special_coul;
+  double *special_lj = force->special_lj;
+  int newton_pair = force->newton_pair;
+  double qqrd2e = force->qqrd2e;
+
+  // loop over neighbors of my atoms
+
+  for (i = 0; i < nlocal; i++) {
+    qtmp = q[i];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    neighs = neighbor->firstneigh[i];
+    numneigh = neighbor->numneigh[i];
+
+    for (k = 0; k < numneigh; k++) {
+      j = neighs[k];
+
+      if (j < nall) factor_coul = factor_lj = 1.0;
+      else {
+	factor_coul = special_coul[j/nall];
+	factor_lj = special_lj[j/nall];
+	j = j % nall;
+      }
+
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+      jtype = type[j];
+
+      if (rsq < cutsq[itype][jtype]) {
+	r2inv = 1.0/rsq;
+	rinv = sqrt(r2inv);
+
+	// atom can have both a charge and dipole
+	// i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole
+	
+	if (rsq < cut_coulsq[itype][jtype]) {
+	  
+	  forcecoulx = forcecouly = forcecoulz = 0.0;
+	  tixcoul = tiycoul = tizcoul = 0.0;
+	  tjxcoul = tjycoul = tjzcoul = 0.0;
+
+	  if (qtmp != 0.0 && q[j] != 0.0) {
+            r3inv = r2inv*rinv;
+	    pre1 = qtmp*q[j]*r3inv;
+
+	    forcecoulx += pre1*delx;
+	    forcecouly += pre1*dely;
+	    forcecoulz += pre1*delz;
+	  }
+
+	  if (dipole[itype] > 0.0 && dipole[jtype] > 0.0) { 
+            r3inv = r2inv*rinv;
+            r5inv = r3inv*r2inv;
+	    r7inv = r5inv*r2inv;
+
+            pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
+            pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
+            pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
+
+	    pre1 = 3.0*r5inv*pdotp - 15.0*r7inv*pidotr*pjdotr;
+	    pre2 = 3.0*r5inv*pjdotr;
+	    pre3 = 3.0*r5inv*pidotr;
+	    pre4 = -1.0*r3inv;
+
+	    forcecoulx += pre1*delx + pre2*mu[i][0] + pre3*mu[j][0];
+	    forcecouly += pre1*dely + pre2*mu[i][1] + pre3*mu[j][1];
+	    forcecoulz += pre1*delz + pre2*mu[i][2] + pre3*mu[j][2];
+	    
+	    crossx = pre4 * (mu[i][1]*mu[j][2] - mu[i][2]*mu[j][1]);
+	    crossy = pre4 * (mu[i][2]*mu[j][0] - mu[i][0]*mu[j][2]);
+	    crossz = pre4 * (mu[i][0]*mu[j][1] - mu[i][1]*mu[j][0]);
+
+	    tixcoul += crossx + pre2 * (mu[i][1]*delz - mu[i][2]*dely);
+	    tiycoul += crossy + pre2 * (mu[i][2]*delx - mu[i][0]*delz);
+	    tizcoul += crossz + pre2 * (mu[i][0]*dely - mu[i][1]*delx);
+	    tjxcoul += -crossx + pre3 * (mu[j][1]*delz - mu[j][2]*dely);
+	    tjycoul += -crossy + pre3 * (mu[j][2]*delx - mu[j][0]*delz);
+	    tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx);
+	  }
+
+	  if (dipole[itype] > 0.0 && q[j] != 0.0) { 
+            r3inv = r2inv*rinv;
+            r5inv = r3inv*r2inv;
+            pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
+	    pre1 = 3.0*q[j]*r5inv * pidotr;
+	    pre2 = q[j]*r3inv;
+
+	    forcecoulx += pre2*mu[i][0] - pre1*delx;
+            forcecouly += pre2*mu[i][1] - pre1*dely;
+            forcecoulz += pre2*mu[i][2] - pre1*delz;
+	    tixcoul += pre2 * (mu[i][1]*delz - mu[i][2]*dely);
+	    tiycoul += pre2 * (mu[i][2]*delx - mu[i][0]*delz);
+	    tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx);
+	  }
+
+	  if (dipole[jtype] > 0.0 && qtmp != 0.0) { 
+            r3inv = r2inv*rinv;
+            r5inv = r3inv*r2inv;
+            pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
+	    pre1 = 3.0*qtmp*r5inv * pjdotr;
+	    pre2 = qtmp*r3inv;
+
+	    forcecoulx += pre1*delx - pre2*mu[j][0];
+            forcecouly += pre1*dely - pre2*mu[j][1];
+            forcecoulz += pre1*delz - pre2*mu[j][2];
+	    tjxcoul += -pre2 * (mu[j][1]*delz - mu[j][2]*dely);
+	    tjycoul += -pre2 * (mu[j][2]*delx - mu[j][0]*delz);
+	    tjzcoul += -pre2 * (mu[j][0]*dely - mu[j][1]*delx);
+	  }
+	}
+
+	// LJ interaction
+
+	if (rsq < cut_ljsq[itype][jtype]) {
+	  r6inv = r2inv*r2inv*r2inv;
+	  forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+	  fforce = factor_lj * forcelj*r2inv;
+	} else fforce = 0.0;
+	  
+	// total force
+
+	fq = factor_coul*qqrd2e;
+	fx = fq*forcecoulx + delx*fforce;
+	fy = fq*forcecouly + dely*fforce;
+	fz = fq*forcecoulz + delz*fforce;
+	
+	// force & torque accumulation
+
+	f[i][0] += fx;
+	f[i][1] += fy;
+	f[i][2] += fz;
+	torque[i][0] += fq*tixcoul;
+	torque[i][1] += fq*tiycoul;
+	torque[i][2] += fq*tizcoul;
+
+	if (newton_pair || j < nlocal) {
+	  f[j][0] -= fx;
+	  f[j][1] -= fy;
+	  f[j][2] -= fz;
+	  torque[j][0] += fq*tjxcoul;
+	  torque[j][1] += fq*tjycoul;
+	  torque[j][2] += fq*tjzcoul;
+	}
+
+	if (eflag) {
+	  if (newton_pair || j < nlocal) factor = 1.0;
+	  else factor = 0.5;
+
+	  if (rsq < cut_coulsq[itype][jtype]) {
+	    phicoul = qtmp*q[j]*rinv;
+	    if (dipole[itype] > 0.0 && dipole[jtype] > 0.0)
+	      phicoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr;
+	    if (dipole[itype] > 0.0 && q[j] != 0.0) 
+	      phicoul += -q[j]*r3inv*pidotr;
+	    if (dipole[jtype] > 0.0 && qtmp != 0.0)
+	      phicoul += qtmp*r3inv*pjdotr;
+	    eng_coul += factor*factor_coul*qqrd2e*phicoul;
+	  }
+
+	  if (rsq < cut_ljsq[itype][jtype]) {
+	    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+	      offset[itype][jtype];
+	    eng_vdwl += factor*factor_lj*philj;
+	  }
+	}
+
+	if (vflag == 1) {
+	  if (newton_pair == 0 && j >= nlocal) {
+	    fx *= 0.5; fy *= 0.5; fz *= 0.5;
+	  }
+	  virial[0] += delx*fx;
+	  virial[1] += dely*fy;
+	  virial[2] += delz*fz;
+	  virial[3] += delx*fy;
+	  virial[4] += delx*fz;
+	  virial[5] += dely*fz;
+	}
+      }
+    }
+  }
+  if (vflag == 2) virial_compute();
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays 
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+
+  cut_lj = memory->create_2d_double_array(n+1,n+1,"pair:cut_lj");
+  cut_ljsq = memory->create_2d_double_array(n+1,n+1,"pair:cut_ljsq");
+  cut_coul = memory->create_2d_double_array(n+1,n+1,"pair:cut_coul");
+  cut_coulsq = memory->create_2d_double_array(n+1,n+1,"pair:cut_coulsq");
+  epsilon = memory->create_2d_double_array(n+1,n+1,"pair:epsilon");
+  sigma = memory->create_2d_double_array(n+1,n+1,"pair:sigma");
+  lj1 = memory->create_2d_double_array(n+1,n+1,"pair:lj1");
+  lj2 = memory->create_2d_double_array(n+1,n+1,"pair:lj2");
+  lj3 = memory->create_2d_double_array(n+1,n+1,"pair:lj3");
+  lj4 = memory->create_2d_double_array(n+1,n+1,"pair:lj4");
+  offset = memory->create_2d_double_array(n+1,n+1,"pair:offset");
+}
+
+/* ----------------------------------------------------------------------
+   global settings 
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::settings(int narg, char **arg)
+{
+  if (narg < 1 || narg > 2)
+    error->all("Incorrect args in pair_style command");
+
+  cut_lj_global = atof(arg[0]);
+  if (narg == 1) cut_coul_global = cut_lj_global;
+  else cut_coul_global = atof(arg[1]);
+
+  // reset cutoffs that have been explicitly set
+
+  if (allocated) {
+    int i,j;
+    for (i = 1; i <= atom->ntypes; i++)
+      for (j = i+1; j <= atom->ntypes; j++)
+	if (setflag[i][j]) {
+	  cut_lj[i][j] = cut_lj_global;
+	  cut_coul[i][j] = cut_coul_global;
+	}
+  }
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::coeff(int narg, char **arg)
+{
+  if (narg < 4 || narg > 6) 
+    error->all("Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi,jlo,jhi;
+  force->bounds(arg[0],atom->ntypes,ilo,ihi);
+  force->bounds(arg[1],atom->ntypes,jlo,jhi);
+
+  double epsilon_one = atof(arg[2]);
+  double sigma_one = atof(arg[3]);
+
+  double cut_lj_one = cut_lj_global;
+  double cut_coul_one = cut_coul_global;
+  if (narg >= 5) cut_coul_one = cut_lj_one = atof(arg[4]);
+  if (narg == 6) cut_coul_one = atof(arg[5]);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      epsilon[i][j] = epsilon_one;
+      sigma[i][j] = sigma_one;
+      cut_lj[i][j] = cut_lj_one;
+      cut_coul[i][j] = cut_coul_one;
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairDipoleCut::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) {
+    epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
+			       sigma[i][i],sigma[j][j]);
+    sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
+    cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
+    cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]);
+  }
+
+  double cut = MAX(cut_lj[i][j],cut_coul[i][j]);
+  cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
+  cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j];
+
+  lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+  lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+  lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
+  lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
+     
+  if (offset_flag) {
+    double ratio = sigma[i][j] / cut_lj[i][j];
+    offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
+  } else offset[i][j] = 0.0;
+  
+  cut_ljsq[j][i] = cut_ljsq[i][j];
+  cut_coulsq[j][i] = cut_coulsq[i][j];
+  lj1[j][i] = lj1[i][j];
+  lj2[j][i] = lj2[i][j];
+  lj3[j][i] = lj3[i][j];
+  lj4[j][i] = lj4[i][j];
+  offset[j][i] = offset[i][j];
+
+  return cut;
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::init_style()
+{
+  if (atom->q == NULL || atom->mu == NULL || 
+      atom->torque == NULL || atom->dipole == NULL)
+    error->all("Pair dipole/cut requires atom attributes "
+	       "q, mu, torque, dipole");
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+      if (setflag[i][j]) {
+	fwrite(&epsilon[i][j],sizeof(double),1,fp);
+	fwrite(&sigma[i][j],sizeof(double),1,fp);
+	fwrite(&cut_lj[i][j],sizeof(double),1,fp);
+	fwrite(&cut_coul[i][j],sizeof(double),1,fp);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+      if (setflag[i][j]) {
+	if (me == 0) {
+	  fread(&epsilon[i][j],sizeof(double),1,fp);
+	  fread(&sigma[i][j],sizeof(double),1,fp);
+	  fread(&cut_lj[i][j],sizeof(double),1,fp);
+	  fread(&cut_coul[i][j],sizeof(double),1,fp);
+	}
+	MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
+	MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
+	MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
+	MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
+      }
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_lj_global,sizeof(double),1,fp);
+  fwrite(&cut_coul_global,sizeof(double),1,fp);
+  fwrite(&offset_flag,sizeof(int),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairDipoleCut::read_restart_settings(FILE *fp)
+{
+  if (comm->me == 0) {
+    fread(&cut_lj_global,sizeof(double),1,fp);
+    fread(&cut_coul_global,sizeof(double),1,fp);
+    fread(&offset_flag,sizeof(int),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairDipoleCut::single(int i, int j, int itype, int jtype, double rsq,
+			   double factor_coul, double factor_lj, int eflag,
+			   One &one)
+{
+  double rinv,r2inv,r6inv,r3inv,r5inv,r7inv;
+  double forcecoulx,forcecouly,forcecoulz,fforce,crossx,crossy,crossz;
+  double tixcoul,tiycoul,tizcoul,tjxcoul,tjycoul,tjzcoul;
+  double pdotp,pdotr,pidotr,pjdotr,pre1,pre2,pre3,pre4;
+  double fq,forcelj,phicoul,philj;
+
+  double delx = atom->x[i][0] - atom->x[j][0];
+  double dely = atom->x[i][1] - atom->x[j][1];
+  double delz = atom->x[i][2] - atom->x[j][2];
+
+  r2inv = 1.0/rsq;
+  rinv = sqrt(r2inv);
+
+  forcecoulx = forcecouly = forcecoulz = 0.0;
+  tixcoul = tiycoul = tizcoul = 0.0;
+  tjxcoul = tjycoul = tjzcoul = 0.0;
+  
+  if (rsq < cut_coulsq[itype][jtype]) {
+    double **mu = atom->mu;
+
+    if (atom->q[i] != 0.0 && atom->q[j] != 0.0) {
+      r3inv = r2inv*rinv;
+      pre1 = atom->q[i]*atom->q[j]*r3inv;
+      
+      forcecoulx += pre1*delx;
+      forcecouly += pre1*dely;
+      forcecoulz += pre1*delz;
+    }
+
+    if (atom->dipole[itype] > 0.0 && atom->dipole[jtype] > 0.0) { 
+      r3inv = r2inv*rinv;
+      r5inv = r3inv*r2inv;
+      r7inv = r5inv*r2inv;
+      
+      pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
+      pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
+      pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
+      
+      pre1 = 3.0*r5inv*pdotp - 15.0*r7inv*pidotr*pjdotr;
+      pre2 = 3.0*r5inv*pjdotr;
+      pre3 = 3.0*r5inv*pidotr;
+      pre4 = -1.0*r3inv;
+      
+      forcecoulx += pre1*delx + pre2*mu[i][0] + pre3*mu[j][0];
+      forcecouly += pre1*dely + pre2*mu[i][1] + pre3*mu[j][1];
+      forcecoulz += pre1*delz + pre2*mu[i][2] + pre3*mu[j][2];
+      
+      crossx = pre4 * (mu[i][1]*mu[j][2] - mu[i][2]*mu[j][1]);
+      crossy = pre4 * (mu[i][2]*mu[j][0] - mu[i][0]*mu[j][2]);
+      crossz = pre4 * (mu[i][0]*mu[j][1] - mu[i][1]*mu[j][0]);
+      
+      tixcoul += crossx + pre2 * (mu[i][1]*delz - mu[i][2]*dely);
+      tiycoul += crossy + pre2 * (mu[i][2]*delx - mu[i][0]*delz);
+      tizcoul += crossz + pre2 * (mu[i][0]*dely - mu[i][1]*delx);
+      tjxcoul += -crossx + pre3 * (mu[j][1]*delz - mu[j][2]*dely);
+      tjycoul += -crossy + pre3 * (mu[j][2]*delx - mu[j][0]*delz);
+      tjzcoul += -crossz + pre3 * (mu[j][0]*dely - mu[j][1]*delx);
+      
+    } else if (atom->dipole[itype] > 0.0 && atom->q[j] != 0.0) { 
+      r3inv = r2inv*rinv;
+      r5inv = r3inv*r2inv;
+      pdotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
+      pre1 = 3.0*atom->q[j]*r5inv * pdotr;
+      pre2 = atom->q[j]*r3inv;
+      
+      forcecoulx += pre2*mu[i][0] - pre1*delx;
+      forcecouly += pre2*mu[i][1] - pre1*dely;
+      forcecoulz += pre2*mu[i][2] - pre1*delz;
+      tixcoul += pre2 * (mu[i][1]*delz - mu[i][2]*dely);
+      tiycoul += pre2 * (mu[i][2]*delx - mu[i][0]*delz);
+      tizcoul += pre2 * (mu[i][0]*dely - mu[i][1]*delx);
+      
+    } else if (atom->dipole[jtype] > 0.0 && atom->q[i] != 0.0) { 
+      r3inv = r2inv*rinv;
+      r5inv = r3inv*r2inv;
+      pdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
+      pre1 = 3.0*atom->q[i]*r5inv * pdotr;
+      pre2 = atom->q[i]*r3inv;
+      
+      forcecoulx += pre1*delx - pre2*mu[j][0];
+      forcecouly += pre1*dely - pre2*mu[j][1];
+      forcecoulz += pre1*delz - pre2*mu[j][2];
+      tjxcoul += -pre2 * (mu[j][1]*delz - mu[j][2]*dely);
+      tjycoul += -pre2 * (mu[j][2]*delx - mu[j][0]*delz);
+	    tjzcoul += -pre2 * (mu[j][0]*dely - mu[j][1]*delx);
+    }
+  }
+
+  if (rsq < cut_ljsq[itype][jtype]) {
+    r6inv = r2inv*r2inv*r2inv;
+    forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+    fforce = factor_lj * forcelj*r2inv;
+  } else fforce = 0.0;
+	  
+  fq = factor_coul*force->qqrd2e;
+  one.fx = fq*forcecoulx + delx*fforce;
+  one.fy = fq*forcecouly + dely*fforce;
+  one.fz = fq*forcecoulz + delz*fforce;
+  one.tix = fq*tixcoul;
+  one.tiy = fq*tiycoul;
+  one.tiz = fq*tizcoul;
+  one.tjx = fq*tjxcoul;
+  one.tjy = fq*tjycoul;
+  one.tjz = fq*tjzcoul;
+	
+  if (eflag) {
+    if (rsq < cut_coulsq[itype][jtype]) {
+      phicoul = atom->q[i]*atom->q[j]*rinv;
+      if (atom->dipole[itype] > 0.0 && atom->dipole[jtype] > 0.0)
+	phicoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr;
+      else if (atom->dipole[itype] > 0.0 && atom->q[j] != 0.0)
+	phicoul += -pre2*pdotr;
+      else if (atom->dipole[jtype] > 0.0 && atom->q[i] != 0.0)
+	phicoul += pre2*pdotr;
+      one.eng_coul = factor_coul*force->qqrd2e*phicoul;
+    } else one.eng_coul = 0.0;
+    if (rsq < cut_ljsq[itype][jtype]) {
+      philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
+	offset[itype][jtype];
+      one.eng_vdwl = factor_lj*philj;
+    } else one.eng_vdwl = 0.0;
+  }
+}

src/DIPOLE/pair_dipole_cut.h

@@ -0,0 +1,48 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef PAIR_DIPOLE_CUT_H
+#define PAIR_DIPOLE_CUT_H
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairDipoleCut : public Pair {
+ public:
+    PairDipoleCut(class LAMMPS *);
+  ~PairDipoleCut();
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  void single(int, int, int, int, double, double, double, int, One &);
+
+ private:
+  double cut_lj_global,cut_coul_global;
+  double **cut_lj,**cut_ljsq;
+  double **cut_coul,**cut_coulsq;
+  double **epsilon,**sigma;
+  double **lj1,**lj2,**lj3,**lj4,**offset;
+
+  void allocate();
+};
+
+}
+
+#endif

src/DIPOLE/style_dipole.h

@@ -0,0 +1,44 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef AtomInclude
+#include "atom_vec_dipole.h"
+#endif
+
+#ifdef AtomClass
+AtomStyle(dipole,AtomVecDipole)
+#endif
+
+#ifdef ComputeInclude
+#include "compute_temp_dipole.h"
+#endif
+
+#ifdef ComputeClass
+ComputeStyle(temp/dipole,ComputeTempDipole)
+#endif
+
+#ifdef FixInclude
+#include "fix_nve_dipole.h"
+#endif
+
+#ifdef FixClass
+FixStyle(nve/dipole,FixNVEDipole)
+#endif
+
+#ifdef PairInclude
+#include "pair_dipole_cut.h"
+#endif
+
+#ifdef PairClass
+PairStyle(dipole/cut,PairDipoleCut)
+#endif

src/MANYBODY/Install.csh

@@ -4,14 +4,12 @@ if ($1 == 1) then
 
   cp style_manybody.h ..
 
-  cp pair_airebo.cpp ..
   cp pair_eam.cpp ..
   cp pair_eam_alloy.cpp ..
   cp pair_eam_fs.cpp ..
   cp pair_sw.cpp ..
   cp pair_tersoff.cpp ..
 
-  cp pair_airebo.h ..
   cp pair_eam.h ..
   cp pair_eam_alloy.h ..
   cp pair_eam_fs.h ..
@@ -23,14 +21,12 @@ else if ($1 == 0) then
   rm ../style_manybody.h
   touch ../style_manybody.h
 
-  rm ../pair_airebo.cpp
   rm ../pair_eam.cpp
   rm ../pair_eam_alloy.cpp
   rm ../pair_eam_fs.cpp
   rm ../pair_sw.cpp
   rm ../pair_tersoff.cpp
 
-  rm ../pair_airebo.h
   rm ../pair_eam.h
   rm ../pair_eam_alloy.h
   rm ../pair_eam_fs.h

src/MANYBODY/style_manybody.h

@@ -12,7 +12,6 @@
 ------------------------------------------------------------------------- */
 
 #ifdef PairInclude
-#include "pair_airebo.h"
 #include "pair_eam.h"
 #include "pair_eam_alloy.h"
 #include "pair_eam_fs.h"
@@ -21,7 +20,6 @@
 #endif
 
 #ifdef PairClass
-PairStyle(airebo,PairAIREBO)
 PairStyle(eam,PairEAM)
 PairStyle(eam/alloy,PairEAMAlloy)
 PairStyle(eam/fs,PairEAMFS)

src/compute_temp_deform.cpp

@@ -0,0 +1,193 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Pieter in 't Veld (SNL)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "string.h"
+#include "compute_temp_deform.h"
+#include "domain.h"
+#include "atom.h"
+#include "force.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_deform.h"
+#include "group.h"
+#include "comm.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("Illegal compute temp/deform command");
+
+  scalar_flag = vector_flag = 1;
+  size_vector = 6;
+  extensive = 0;
+  tempflag = 1;
+
+  vector = new double[6];
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTempDeform::~ComputeTempDeform()
+{
+  delete [] vector;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDeform::init()
+{
+  int i;
+
+  fix_dof = 0;
+  for (i = 0; i < modify->nfix; i++)
+    fix_dof += modify->fix[i]->dof(igroup);
+  recount();
+
+  // 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 != V_REMAP && 
+	  comm->me == 0)
+	error->warning("Using fix nvt/sllod with inconsistent fix deform remapping");
+      break;
+    }
+  if (i == modify->nfix && comm->me == 0)
+    error->warning("Using fix nvt/sllod with no fix deform defined");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDeform::recount()
+{
+  double natoms = group->count(igroup);
+  dof = force->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];
+
+  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;
+
+  if (mass) {
+    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];
+	t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + 
+	      vthermal[2]*vthermal[2]) * mass[type[i]];
+      }
+  }
+  else
+    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];
+	t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] + 
+	      vthermal[2]*vthermal[2]) * rmass[i];
+      }
+
+  MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
+  scalar *= tfactor;
+  return scalar;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeTempDeform::compute_vector()
+{
+  double lamda[3],vstream[3],vthermal[3];
+
+  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 (mass) massone = mass[type[i]];
+      else massone = rmass[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;
+}

src/compute_temp_deform.h

@@ -0,0 +1,38 @@
+/* ----------------------------------------------------------------------
+   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 COMPUTE_TEMP_DEFORM_H
+#define COMPUTE_TEMP_DEFORM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeTempDeform : public Compute {
+ public:
+  ComputeTempDeform(class LAMMPS *, int, char **);
+  ~ComputeTempDeform();
+  void init();
+  double compute_scalar();
+  void compute_vector();
+
+ private:
+  int fix_dof;
+  double tfactor;
+
+  void recount();
+};
+
+}
+
+#endif

src/fix_deform.cpp

@@ -0,0 +1,674 @@
+/* ----------------------------------------------------------------------
+   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 "string.h"
+#include "stdlib.h"
+#include "math.h"
+#include "fix_deform.h"
+#include "atom.h"
+#include "update.h"
+#include "comm.h"
+#include "domain.h"
+#include "lattice.h"
+#include "force.h"
+#include "modify.h"
+#include "kspace.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{NONE,FINAL,DELTA,SCALE,VEL,RATE,VOLUME};
+enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE};
+
+// same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp
+
+enum{NO_REMAP,X_REMAP,V_REMAP};
+
+/* ---------------------------------------------------------------------- */
+
+FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+  if (narg < 4) error->all("Illegal fix deform command");
+
+  box_change = 1;
+
+  nevery = atoi(arg[3]);
+  if (nevery <= 0) error->all("Illegal fix deform command");
+
+  // set defaults
+
+  set = new Set[6];
+  set[0].style = set[1].style = set[2].style = 
+    set[3].style = set[4].style = set[5].style = NONE;
+
+  // parse arguments
+
+  triclinic = domain->triclinic;
+
+  int index;
+  int iarg = 4;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"x") == 0 || strcmp(arg[iarg],"y") == 0 ||
+	strcmp(arg[iarg],"z") == 0) {
+      if (strcmp(arg[iarg],"x") == 0) index = 0;
+      else if (strcmp(arg[iarg],"y") == 0) index = 1;
+      else if (strcmp(arg[iarg],"z") == 0) index = 2;
+
+      if (iarg+2 > narg) error->all("Illegal fix deform command");
+      if (strcmp(arg[iarg+1],"final") == 0) {
+	if (iarg+4 > narg) error->all("Illegal fix deform command");
+	set[index].style = FINAL;
+	set[index].flo = atof(arg[iarg+2]);
+	set[index].fhi = atof(arg[iarg+3]);
+	iarg += 4;
+      } else if (strcmp(arg[iarg+1],"delta") == 0) {
+	if (iarg+4 > narg) error->all("Illegal fix deform command");
+	set[index].style = DELTA;
+	set[index].dlo = atof(arg[iarg+2]);
+	set[index].dhi = atof(arg[iarg+3]);
+	iarg += 4;
+      } else if (strcmp(arg[iarg+1],"scale") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = SCALE;
+	set[index].scale = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"vel") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = VEL;
+	set[index].vel = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"rate") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = RATE;
+	set[index].rate = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"volume") == 0) {
+	set[index].style = VOLUME;
+	iarg += 2;
+      } else error->all("Illegal fix deform command");
+
+    } else if (strcmp(arg[iarg],"xy") == 0 || strcmp(arg[iarg],"xz") == 0 ||
+	strcmp(arg[iarg],"yz") == 0) {
+      if (triclinic == 0)
+	error->all("Fix deform tilt factors require triclinic box");
+      if (strcmp(arg[iarg],"xy") == 0) index = 5;
+      else if (strcmp(arg[iarg],"xz") == 0) index = 4;
+      else if (strcmp(arg[iarg],"yz") == 0) index = 3;
+      if (iarg+2 > narg) error->all("Illegal fix deform command");
+      if (strcmp(arg[iarg+1],"final") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = FINAL;
+	set[index].ftilt = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"delta") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = DELTA;
+	set[index].dtilt = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"vel") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = VEL;
+	set[index].vel = atof(arg[iarg+2]);
+	iarg += 3;
+      } else if (strcmp(arg[iarg+1],"rate") == 0) {
+	if (iarg+3 > narg) error->all("Illegal fix deform command");
+	set[index].style = RATE;
+	set[index].rate = atof(arg[iarg+2]);
+	iarg += 3;
+      } else error->all("Illegal fix deform command");
+    } else break;
+  }
+
+  // read options from end of input line
+
+  options(narg-iarg,&arg[iarg]);
+
+  // check periodicity
+
+  if ((set[0].style && domain->xperiodic == 0) ||
+      (set[1].style && domain->yperiodic == 0) ||
+      (set[2].style && domain->zperiodic == 0))
+    error->all("Cannot fix deform on a non-periodic boundary");
+
+  if (set[3].style && (domain->yperiodic == 0 || domain->zperiodic == 0))
+    error->all("Cannot fix deform on a non-periodic boundary");
+  if (set[4].style && (domain->xperiodic == 0 || domain->zperiodic == 0))
+    error->all("Cannot fix deform on a non-periodic boundary");
+  if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0))
+    error->all("Cannot fix deform on a non-periodic boundary");
+
+  // setup scaling
+
+  if (scaleflag && domain->lattice == NULL)
+    error->all("Use of fix deform with undefined lattice");
+
+  if (scaleflag) {
+    xscale = domain->lattice->xlattice;
+    yscale = domain->lattice->ylattice;
+    zscale = domain->lattice->zlattice;
+  }
+  else xscale = yscale = zscale = 1.0;
+
+  // apply scaling to input parameters with dist/vel units
+  // for 3,4,5 scale is in 1st dimension, e.g. x for xz
+
+  double map[6];
+  map[0] = xscale; map[1] = yscale; map[2] = zscale;
+  map[3] = yscale; map[4] = xscale; map[5] = xscale;
+
+  for (int i = 0; i < 3; i++) {
+    if (set[i].style == FINAL) {
+      set[i].flo *= map[i];
+      set[i].fhi *= map[i];
+    } else if (set[i].style == DELTA) {
+      set[i].dlo *= map[i];
+      set[i].dhi *= map[i];
+    } else if (set[i].style == VEL) {
+      set[i].vel *= map[i];
+    }
+  }
+
+  for (int i = 3; i < 6; i++) {
+    if (set[i].style == FINAL) set[i].ftilt *= map[i];
+    else if (set[i].style == DELTA) set[i].dtilt *= map[i];
+    else if (set[i].style == VEL) set[i].vel *= map[i];
+  }
+
+  // set initial/final values for box size and shape for FINAL,DELTA,SCALE
+  // final not possible for VEL,RATE since don't know length of run yet
+  // final = initial if no setting
+
+  for (int i = 0; i < 3; i++) {
+    set[i].lo_target = set[i].lo_stop = set[i].lo_start = domain->boxlo[i];
+    set[i].hi_target = set[i].hi_stop = set[i].hi_start = domain->boxhi[i];
+
+    if (set[i].style == FINAL) {
+      set[i].lo_stop = set[i].flo;
+      set[i].hi_stop = set[i].fhi;
+    } else if (set[i].style == DELTA) {
+      set[i].lo_stop = set[i].lo_start + set[i].dlo;
+      set[i].hi_stop = set[i].hi_start + set[i].dhi;
+    } else if (set[i].style == SCALE) {
+      set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - 
+	0.5*set[i].scale*(set[i].hi_start-set[i].lo_start);
+      set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + 
+	0.5*set[i].scale*(set[i].hi_start-set[i].lo_start);
+    }
+  }
+
+  for (int i = 3; i < 6; i++) {
+    if (i == 5) set[i].tilt_start = domain->xy;
+    else if (i == 4) set[i].tilt_start = domain->xz;
+    else if (i == 3) set[i].tilt_start = domain->yz;
+    set[i].tilt_flip = set[i].tilt_target = 
+      set[i].tilt_stop = set[i].tilt_start;
+
+    if (set[i].style == FINAL) {
+      set[i].tilt_stop = set[i].ftilt;
+    } else if (set[i].style == DELTA) {
+      set[i].tilt_stop = set[i].tilt_start + set[i].dtilt;
+    }
+  }
+
+  // for VOLUME, setup links to other dims
+  // fixed, dynamic1,2, vol_start
+
+  for (int i = 0; i < 3; i++) {
+    set[i].vol_start = domain->xprd * domain->yprd * domain->zprd;
+
+    if (set[i].style != VOLUME) continue;
+    int other1 = (i+1) % 3;
+    int other2 = (i+2) % 3;
+
+    if (set[other1].style == NONE) {
+      if (set[other2].style == NONE || set[other2].style == VOLUME)
+	error->all("Fix deform volume setting is invalid");
+      set[i].substyle = ONE_FROM_ONE;
+      set[i].fixed = other1;
+      set[i].dynamic1 = other2;
+    } else if (set[other2].style == NONE) {
+      if (set[other1].style == NONE || set[other1].style == VOLUME)
+	error->all("Fix deform volume setting is invalid");
+      set[i].substyle = ONE_FROM_ONE;
+      set[i].fixed = other2;
+      set[i].dynamic1 = other1;
+    } else if (set[other1].style == VOLUME) {
+      if (set[other2].style == NONE || set[other2].style == VOLUME)
+	error->all("Fix deform volume setting is invalid");
+      set[i].substyle = TWO_FROM_ONE;
+      set[i].fixed = other1;
+      set[i].dynamic1 = other2;
+    } else if (set[other2].style == VOLUME) {
+      if (set[other1].style == NONE || set[other1].style == VOLUME)
+	error->all("Fix deform volume setting is invalid");
+      set[i].substyle = TWO_FROM_ONE;
+      set[i].fixed = other2;
+      set[i].dynamic1 = other1;
+    } else {
+      set[i].substyle = ONE_FROM_TWO;
+      set[i].dynamic2 = other1;
+      set[i].dynamic2 = other2;
+    }
+  }
+
+  // initial settings
+  // reneighboring only forced if flips will occur due to shape changes
+
+  if (set[3].style || set[4].style || set[5].style) force_reneighbor = 1;
+  next_reneighbor = -1;
+
+  nrigid = 0;
+  rfix = NULL;
+  flip = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixDeform::~FixDeform()
+{
+  delete [] set;
+  delete [] rfix;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixDeform::setmask()
+{
+  int mask = 0;
+  mask |= PRE_EXCHANGE;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixDeform::init()
+{
+  // error if more than one fix deform
+  // domain, fix nvt/sllod, compute temp/deform only work on single h_rate
+
+  int count = 0;
+  for (int i = 0; i < modify->nfix; i++)
+    if (strcmp(modify->fix[i]->style,"deform") == 0) count++;
+  if (count > 1) error->warning("More than one fix deform");
+
+  // Kspace setting
+
+  if (force->kspace) kspace_flag = 1;
+  else kspace_flag = 0;
+
+  // elapsed time for entire simulation, including multiple runs if defined
+  // cannot be 0.0 since can't deform a finite distance in 0 time
+
+  double delt = (update->endstep - update->beginstep) * update->dt;
+  if (delt == 0.0) error->all("Cannot use fix deform for 0 timestep run");
+
+  // set final values for box size and shape for VEL,RATE
+  // now possible since length of run is known
+
+  for (int i = 0; i < 3; i++) {
+    if (set[i].style == VEL) {
+      set[i].lo_stop = set[i].lo_start - 0.5*delt*set[i].vel;
+      set[i].hi_stop = set[i].hi_start + 0.5*delt*set[i].vel;
+    } else if (set[i].style == RATE) {
+      set[i].lo_stop = 0.5*(set[i].lo_start+set[i].hi_start) - 
+	0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt));
+      set[i].hi_stop = 0.5*(set[i].lo_start+set[i].hi_start) + 
+	0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt));
+    }
+  }
+  for (int i = 3; i < 6; i++) {
+    if (set[i].style == VEL) {
+      set[i].tilt_stop = set[i].tilt_start + delt*set[i].vel;
+    } else if (set[i].style == RATE) {
+      set[i].tilt_stop = set[i].tilt_start * pow(1.0+set[i].rate,delt);
+    }
+  }
+
+  // if using tilt RATE, then initial tilt must be non-zero
+
+  for (int i = 3; i < 6; i++)
+    if (set[i].style == RATE)
+      if ((i == 5 && domain->xy == 0.0) || 
+	  (i == 4 && domain->xz == 0.0) || 
+	  (i == 3 && domain->yz == 0.0)) 
+	error->all("Cannot use fix deform rate to tilt a box with zero tilt");
+
+  // if yz changes and will cause box flip, then xy cannot be changing
+  // this is b/c the flips would induce continuous changes in xz
+  //   in order to keep the edge vectors of the flipped shape matrix
+  //   a linear combination of the edge vectors of the unflipped shape matrix
+
+  if (set[3].style && set[5].style)
+    if (set[3].tilt_stop < -0.5*(set[1].hi_start-set[1].lo_start) ||
+	set[3].tilt_stop > 0.5*(set[1].hi_start-set[1].lo_start) ||
+	set[3].tilt_stop < -0.5*(set[1].hi_stop-set[1].lo_stop) ||
+	set[3].tilt_stop > 0.5*(set[1].hi_stop-set[1].lo_stop))
+      error->all("Fix deform is changing yz by too much with changing xy");
+
+  // set domain->h_rate values for use by domain and other fixes/computes
+  // initialize all rates to 0.0
+  // cannot set rate now for RATE,VOLUME styles since not constant
+
+  h_rate = domain->h_rate;
+  h_ratelo = domain->h_ratelo;
+
+  for (int i = 0; i < 3; i++) {
+    h_rate[i] = h_ratelo[i] = 0.0;
+    if (set[i].style == FINAL || set[i].style == DELTA ||
+	set[i].style == SCALE || set[i].style == VEL) {
+      double dlo_dt = (set[i].lo_stop - set[i].lo_start) / delt;
+      double dhi_dt = (set[i].hi_stop - set[i].hi_start) / delt;
+      h_rate[i] = dhi_dt - dlo_dt;
+      h_ratelo[i] = dlo_dt;
+    }
+  }    
+  for (int i = 3; i < 6; i++) {
+    h_rate[i] = 0.0;
+    if (set[i].style == FINAL || set[i].style == DELTA || 
+	set[i].style == VEL)
+      h_rate[i] = (set[i].tilt_stop - set[i].tilt_start) / delt;
+  }    
+
+  // detect if any rigid fixes exist so rigid bodies can be rescaled
+  // 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;
+  }
+}
+
+/* ----------------------------------------------------------------------
+  box flipped on previous step
+  perform irregular comm to migrate atoms to new procs
+  reset box tilts for flipped config and create new box in domain
+  remap to put far-away atoms back into new box
+  perform irregular on atoms in lamda coords to get atoms to new procs
+  force reneighboring on next timestep
+------------------------------------------------------------------------- */
+
+void FixDeform::pre_exchange()
+{
+  if (flip == 0) return;
+  
+  domain->yz = set[3].tilt_target = set[3].tilt_flip;
+  domain->xz = set[4].tilt_target = set[4].tilt_flip;
+  domain->xy = set[5].tilt_target = set[5].tilt_flip;
+  domain->set_global_box();
+  domain->set_local_box();
+  
+  double **x = atom->x;
+  int *image = atom->image;
+  int nlocal = atom->nlocal;
+  for (int i = 0; i < nlocal; i++)
+    domain->remap(x[i],image[i]);
+  
+  domain->x2lamda(atom->nlocal);
+  comm->irregular();
+  domain->lamda2x(atom->nlocal);
+
+  flip = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixDeform::end_of_step()
+{
+  int i;
+
+  double delta = update->ntimestep - update->beginstep;
+  delta /= update->endstep - update->beginstep;
+
+  // set new box size
+  // for RATE, set target directly based on current time and set h_rate
+  // for others except VOLUME, target is linear value between start and stop
+
+  for (i = 0; i < 3; i++) {
+    if (set[i].style == NONE) continue;
+
+    if (set[i].style == RATE) {
+      double delt = (update->ntimestep - update->beginstep) * update->dt;
+      set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) - 
+	0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt));
+      set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) + 
+	0.5*((set[i].hi_start-set[i].lo_start)*pow(1.0+set[i].rate,delt));
+      h_rate[i] = set[i].rate * domain->h[i];
+      h_ratelo[i] = -0.5*h_rate[i];
+
+    } else if (set[i].style != VOLUME) {
+      set[i].lo_target = set[i].lo_start + 
+	delta*(set[i].lo_stop - set[i].lo_start);
+      set[i].hi_target = set[i].hi_start + 
+	delta*(set[i].hi_stop - set[i].hi_start);
+    }
+  }
+
+  // set new box size for VOLUME dims that are linked to other dims
+  // also need to set h_rate
+
+  for (int i = 0; i < 3; i++) {
+    if (set[i].style != VOLUME) continue;
+
+    if (set[i].substyle == ONE_FROM_ONE) {
+      set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) -
+	0.5*(set[i].vol_start /
+	     (set[set[i].dynamic1].hi_target -
+	      set[set[i].dynamic1].lo_target) /
+	     (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start));
+      set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) +
+	0.5*(set[i].vol_start /
+	     (set[set[i].dynamic1].hi_target -
+	      set[set[i].dynamic1].lo_target) /
+	     (set[set[i].fixed].hi_start-set[set[i].fixed].lo_start));
+
+    } else if (set[i].substyle == ONE_FROM_TWO) {
+      set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) -
+	0.5*(set[i].vol_start /
+	     (set[set[i].dynamic1].hi_target - 
+	      set[set[i].dynamic1].lo_target) /
+	     (set[set[i].dynamic2].hi_target - 
+	      set[set[i].dynamic2].lo_target));
+      set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) +
+	0.5*(set[i].vol_start /
+	     (set[set[i].dynamic1].hi_target -
+	      set[set[i].dynamic1].lo_target) /
+	     (set[set[i].dynamic2].hi_target - 
+	      set[set[i].dynamic2].lo_target));
+      
+    } else if (set[i].substyle == TWO_FROM_ONE) {
+      set[i].lo_target = 0.5*(set[i].lo_start+set[i].hi_start) -
+	0.5*sqrt(set[i].vol_start /
+		 (set[set[i].dynamic1].hi_target - 
+		  set[set[i].dynamic1].lo_target) /
+		 (set[set[i].fixed].hi_start - 
+		  set[set[i].fixed].lo_start) *
+		 (set[i].hi_start - set[i].lo_start));
+      set[i].hi_target = 0.5*(set[i].lo_start+set[i].hi_start) +
+	0.5*sqrt(set[i].vol_start /
+		 (set[set[i].dynamic1].hi_target - 
+		  set[set[i].dynamic1].lo_target) /
+		 (set[set[i].fixed].hi_start - 
+		  set[set[i].fixed].lo_start) *
+		 (set[i].hi_start - set[i].lo_start));
+    }
+  }
+
+  // for triclinic, set new box shape
+  // for RATE, set target directly based on current time and set h_rate
+  // for all others, target is linear value between start and stop values
+
+  if (triclinic) {
+    double *h = domain->h;
+
+    for (i = 3; i < 6; i++) {
+      if (set[i].style == NONE) continue;
+
+      if (set[i].style == RATE) {
+	double delt = (update->ntimestep - update->beginstep) * update->dt;
+	set[i].tilt_target = set[i].tilt_start * pow(1.0+set[i].rate,delt);
+	h_rate[i] = set[i].rate * domain->h[i];
+      } else if (set[i].style) {
+	set[i].tilt_target = set[i].tilt_start + 
+	  delta*(set[i].tilt_stop - set[i].tilt_start);
+      }
+
+      // tilt_target can be large positive or large negative value
+      // add/subtract box lengths until tilt_target is closest to current value
+
+      int idenom;
+      if (i == 5) idenom = 0;
+      else if (i == 4) idenom = 0;
+      else if (i == 3) idenom = 1;
+      double denom = set[idenom].hi_target - set[idenom].lo_target;
+
+      double current = h[i]/h[idenom];
+      while (set[i].tilt_target/denom - current > 0.0)
+	set[i].tilt_target -= denom;
+      while (set[i].tilt_target/denom - current < 0.0)
+	set[i].tilt_target += denom;
+      if (fabs(set[i].tilt_target/denom - 1.0 - current) < 
+	  fabs(set[i].tilt_target/denom - current))
+	set[i].tilt_target -= denom;
+    }
+  }
+
+  // if any tilt targets exceed bounds, set flip flag and new tilt_flip values
+  // flip will be performed on next timestep before reneighboring
+  // when yz flips and xy is non-zero, xz must also change
+  // this is to keep the edge vectors of the flipped shape matrix
+  //   a linear combination of the edge vectors of the unflipped shape matrix
+
+  if (triclinic) {
+    double xprd = set[0].hi_target - set[0].lo_target;
+    double yprd = set[1].hi_target - set[1].lo_target;
+    if (set[3].tilt_target < -0.5*yprd || set[3].tilt_target > 0.5*yprd ||
+	set[4].tilt_target < -0.5*xprd || set[4].tilt_target > 0.5*xprd ||
+	set[5].tilt_target < -0.5*xprd || set[5].tilt_target > 0.5*xprd) {
+
+      flip = 1;
+      next_reneighbor = update->ntimestep + 1;
+      set[3].tilt_flip = set[3].tilt_target;
+      set[4].tilt_flip = set[4].tilt_target;
+      set[5].tilt_flip = set[5].tilt_target;
+
+      if (set[3].tilt_flip < -0.5*yprd) {
+	set[3].tilt_flip += yprd;
+	set[4].tilt_flip += set[5].tilt_flip;
+      } else if (set[3].tilt_flip >= 0.5*yprd) {
+	set[3].tilt_flip -= yprd;
+	set[4].tilt_flip -= set[5].tilt_flip;
+      }
+      if (set[4].tilt_flip < -0.5*xprd) set[4].tilt_flip += xprd;
+      if (set[4].tilt_flip > 0.5*xprd) set[4].tilt_flip -= xprd;
+      if (set[5].tilt_flip < -0.5*xprd) set[5].tilt_flip += xprd;
+      if (set[5].tilt_flip > 0.5*xprd) set[5].tilt_flip -= xprd;
+    }
+  }
+
+  // convert atoms to lamda coords
+
+  if (remapflag == X_REMAP) {
+    double **x = atom->x;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	domain->x2lamda(x[i],x[i]);
+ }
+
+  // reset global and local box to new size/shape
+
+  domain->boxlo[0] = set[0].lo_target;
+  domain->boxlo[1] = set[1].lo_target;
+  domain->boxlo[2] = set[2].lo_target;
+  domain->boxhi[0] = set[0].hi_target;
+  domain->boxhi[1] = set[1].hi_target;
+  domain->boxhi[2] = set[2].hi_target;
+
+  if (triclinic) {
+    domain->yz = set[3].tilt_target;
+    domain->xz = set[4].tilt_target;
+    domain->xy = set[5].tilt_target;
+  }
+
+  domain->set_global_box();
+  domain->set_local_box();
+
+  // convert affine atoms back to box coords
+
+  if (remapflag == X_REMAP) {
+    double **x = atom->x;
+    int *mask = atom->mask;
+    int nlocal = atom->nlocal;
+
+    for (i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit)
+	domain->lamda2x(x[i],x[i]);
+  }
+
+  // remap centers-of-mass of rigid bodies
+  // needs to be new call, not dilate
+  // check where else fix dilate is called from
+
+  /*
+  if (nrigid)
+    for (i = 0; i < nrigid; i++)
+      modify->fix[rfix[i]]->dilate(2,oldlo,oldhi,newlo,newhi);
+  */
+
+  // redo KSpace coeffs since box has changed
+
+  if (kspace_flag) force->kspace->setup();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixDeform::options(int narg, char **arg)
+{
+  if (narg < 0) error->all("Illegal fix deform command");
+
+  remapflag = X_REMAP;
+  scaleflag = 1;
+
+  int iarg = 0;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"remap") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix deform command");
+      if (strcmp(arg[iarg+1],"x") == 0) remapflag = X_REMAP;
+      else if (strcmp(arg[iarg+1],"v") == 0) remapflag = V_REMAP;
+      else if (strcmp(arg[iarg+1],"none") == 0) remapflag = NO_REMAP;
+      else error->all("Illegal fix deform command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"units") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix deform command");
+      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+      else error->all("Illegal fix deform command");
+      iarg += 2;
+    } else error->all("Illegal fix deform command");
+  }
+}

src/fix_deform.h

@@ -11,41 +11,54 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#ifndef FIX_UNIAXIAL_H
+#ifndef FIX_DEFORM_H
-#define FIX_UNIAXIAL_H
+#define FIX_DEFORM_H
 
 #include "fix.h"
 
 namespace LAMMPS_NS {
 
-class FixUniaxial : public Fix {
+class FixDeform : public Fix {
  public:
-  FixUniaxial(class LAMMPS *, int, char **);
+  int remapflag;
-  ~FixUniaxial();
+
+  FixDeform(class LAMMPS *, int, char **);
+  ~FixDeform();
   int setmask();
   void init();
+  void pre_exchange();
   void end_of_step();
 
  private:
-  int dir;
+  double xlo_start,xhi_start,ylo_start,yhi_start,zlo_start,zhi_start;
-  double lambda_final;
+  double xlo_stop,xhi_stop,ylo_stop,yhi_stop,zlo_stop,zhi_stop;
-  double alpha0;                   // asymmetry parameter
+  double xy_start,xz_start,yz_start;
+  double xy_stop,xz_stop,yz_stop;
+  double xscale,yscale,zscale;
+  int triclinic,scaleflag,flip;
+  double *h_rate,*h_ratelo;
 
-  double lambdai[3];               // initial
+  struct Set {
-  double lambdaf[3];               // final
+    int style,substyle;
-  double domainloi[3];             // initial box lo/hi
+    double flo,fhi,ftilt;
-  double domainhii[3];
+    double dlo,dhi,dtilt;
-  double *domainlo[3];             // pointers to make loop over dims possible
+    double scale,vel,rate;
-  double *domainhi[3];
+    double lo_start,hi_start;
-  double *domainprd[3];
+    double lo_stop,hi_stop;
+    double lo_target,hi_target;
+    double tilt_start,tilt_stop,tilt_target;
+    double vol_start,tilt_flip;
+    int fixed,dynamic1,dynamic2;
+  };
+  Set *set;
 
   int kspace_flag;                 // 1 if KSpace invoked, 0 if not
   int nrigid;                      // number of rigid fixes
   int *rfix;                       // indices of rigid fixes
-};
 
-#endif
+  void options(int, char **);
+};
 
 }
 
-
+#endif

src/fix_nvt_sllod.cpp

@@ -0,0 +1,278 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Pieter in 't Veld (SNL)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_nvt_sllod.h"
+#include "math_extra.h"
+#include "atom.h"
+#include "force.h"
+#include "domain.h"
+#include "comm.h"
+#include "group.h"
+#include "update.h"
+#include "respa.h"
+#include "modify.h"
+#include "fix.h"
+#include "fix_deform.h"
+#include "compute.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{NO_REMAP,X_REMAP,V_REMAP};                   // same as fix_deform.cpp
+
+/* ---------------------------------------------------------------------- */
+
+FixNVTSlodd::FixNVTSlodd(LAMMPS *lmp, int narg, char **arg) :
+  FixNVT(lmp, narg, arg) {}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTSlodd::init()
+{
+  FixNVT::init();
+
+  // 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("Using compute temp/deform with fix deform remapping coords");
+      break;
+    }
+  if (i == modify->nfix && comm->me == 0)
+    error->warning("Using compute temp/deform with no fix deform defined");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTSlodd::initial_integrate()
+{
+  double dtfm;
+
+  double delta = update->ntimestep - update->firststep;
+  delta /= update->nsteps;
+  t_target = t_start + delta * (t_stop-t_start);
+
+  // update eta_dot
+
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+  eta += dtv*eta_dot;
+  factor = exp(-dthalf*eta_dot);
+
+  // update vthermal and x of only atoms in NVT group
+  // lamda = 0-1 triclinic lamda coords
+  // vstream = streaming velocity = Hrate*lamda + Hratelo
+  // vthermal = thermal velocity = v - vstream
+  // vdelu = Hrate*Hinv*vthermal
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double *h_rate = domain->h_rate;
+  double *h_ratelo = domain->h_ratelo;
+  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
+  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+
+      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];
+      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
+	h_two[4]*vthermal[2];
+      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
+      vdelu[2] = h_two[2]*vthermal[2];
+
+      v[i][0] = vstream[0] + 
+	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
+      v[i][1] = vstream[1] + 
+	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
+      v[i][2] = vstream[2] + 
+	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
+      x[i][0] += dtv * v[i][0];
+      x[i][1] += dtv * v[i][1];
+      x[i][2] += dtv * v[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTSlodd::final_integrate()
+{
+  double dtfm;
+
+  // update vthermal of only atoms in NVT group
+  // lamda = 0-1 triclinic lamda coords
+  // vstream = streaming velocity = Hrate*lamda + Hratelo
+  // vthermal = thermal velocity = v - vstream
+  // vdelu = Hrate*Hinv*vthermal
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  double *h_rate = domain->h_rate;
+  double *h_ratelo = domain->h_ratelo;
+  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
+  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+
+      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];
+      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
+	h_two[4]*vthermal[2];
+      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
+      vdelu[2] = h_two[2]*vthermal[2];
+
+      v[i][0] = vstream[0] + 
+	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
+      v[i][1] = vstream[1] + 
+	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
+      v[i][2] = vstream[2] + 
+	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
+    }
+  }
+
+  // compute current T
+
+  t_current = temperature->compute_scalar();
+
+  // update eta_dot
+
+  f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+  eta_dot += f_eta*dthalf;
+  eta_dot *= drag_factor;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixNVTSlodd::initial_integrate_respa(int ilevel, int flag)
+{
+  if (flag) return;             // only used by NPT,NPH
+
+  // set timesteps by level
+
+  double dtfm;
+  dtv = step_respa[ilevel];
+  dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
+  dthalf = 0.5 * step_respa[ilevel];
+
+  // atom quantities
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double **f = atom->f;
+  double *mass = atom->mass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  // outermost level - update eta_dot and apply to v
+  // all other levels - NVE update of v
+
+  if (ilevel == nlevels_respa-1) {
+    double delta = update->ntimestep - update->firststep;
+    delta /= update->nsteps;
+    t_target = t_start + delta * (t_stop-t_start);
+    
+    // update eta_dot
+
+    f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
+    eta_dot += f_eta*dthalf;
+    eta_dot *= drag_factor;
+    eta += dtv*eta_dot;
+    factor = exp(-dthalf*eta_dot);
+  } else factor = 1.0;
+
+  // update v of only atoms in NVT group
+
+  double *h_rate = domain->h_rate;
+  double *h_ratelo = domain->h_ratelo;
+  double h_two[6],lamda[3],vstream[3],vthermal[3],vdelu[3];
+  MathExtra::multiply_shape_shape(h_rate,domain->h_inv,h_two);
+
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      dtfm = dtf / mass[type[i]];
+
+      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];
+      vdelu[0] = h_two[0]*vthermal[0] + h_two[5]*vthermal[1] + 
+	h_two[4]*vthermal[2];
+      vdelu[1] = h_two[1]*vthermal[1] + h_two[3]*vthermal[2];
+      vdelu[2] = h_two[2]*vthermal[2];
+
+      v[i][0] = vstream[0] + 
+	vthermal[0]*factor + dtfm*f[i][0] - dthalf*vdelu[0];
+      v[i][1] = vstream[1] + 
+	vthermal[1]*factor + dtfm*f[i][1] - dthalf*vdelu[1];
+      v[i][2] = vstream[2] + 
+	vthermal[2]*factor + dtfm*f[i][2] - dthalf*vdelu[2];
+    }
+  }
+
+  // innermost level - also update x of only atoms in NVT group
+
+  if (ilevel == 0) {
+    for (int i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	x[i][0] += dtv * v[i][0];
+	x[i][1] += dtv * v[i][1];
+	x[i][2] += dtv * v[i][2];
+      }
+    }
+  }
+}

src/fix_nvt_sllod.h

@@ -0,0 +1,33 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_NVT_SLODD_H
+#define FIX_NVT_SLODD_H
+
+#include "fix_nvt.h"
+
+namespace LAMMPS_NS {
+
+class FixNVTSlodd : public FixNVT {
+ public:
+  FixNVTSlodd(class LAMMPS *, int, char **);
+  ~FixNVTSlodd() {}
+  void init();
+  void initial_integrate();
+  void final_integrate();
+  void initial_integrate_respa(int,int);
+};
+
+}
+
+#endif

src/fix_uniaxial.cpp

@@ -1,212 +0,0 @@
-/* ----------------------------------------------------------------------
-   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: Carsten Svaneborg
-   (Max Planck Institute for Complex Systems, Dresden, Germany)
-------------------------------------------------------------------------- */
-
-#include "string.h"
-#include "stdlib.h"
-#include "math.h"
-#include "fix_uniaxial.h"
-#include "atom.h"
-#include "update.h"
-#include "force.h"
-#include "domain.h"
-#include "modify.h"
-#include "comm.h"
-#include "kspace.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-FixUniaxial::FixUniaxial(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
-{
-  if (narg != 6) error->all("Illegal fix uniaxial command");
-
-  box_change = 1;
-
-  nevery = atoi(arg[3]);
-  if (nevery <= 0) error->all("Illegal fix uniaxial command");
-
-  if (strcmp(arg[4],"x") == 0) dir = 0;
-  else if (strcmp(arg[4],"y") == 0) dir = 1;
-  else if (strcmp(arg[4],"z") == 0) dir = 2;
-  else error->all("Illegal fix uniaxial command");
-  lambda_final = atof(arg[5]);
-
-  if (lambda_final <= 0) error->all("Illegal fix uniaxial command");
-
-  if (domain->nonperiodic) 
-    error->all("Cannot use fix uniaxial on non-periodic system");
-  if (domain->triclinic)
-    error->all("Cannot use fix uniaxial with triclinic box");
-
-  nrigid = 0;
-  rfix = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixUniaxial::~FixUniaxial()
-{
-  delete [] rfix;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixUniaxial::setmask()
-{
-  int mask = 0;
-  mask |= END_OF_STEP;
-  return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixUniaxial::init()
-{
-  // store pointers to domain variable so can loop over dimensions
-
-  domainlo[0] = &domain->boxlo[0];
-  domainlo[1] = &domain->boxlo[1];
-  domainlo[2] = &domain->boxlo[2];
-  
-  domainhi[0] = &domain->boxhi[0];
-  domainhi[1] = &domain->boxhi[1];
-  domainhi[2] = &domain->boxhi[2];
-
-  domainprd[0] = &domain->xprd;
-  domainprd[1] = &domain->yprd;
-  domainprd[2] = &domain->zprd;
-
-  double L = pow((domain->boxhi[0]-domain->boxlo[0])*
-                 (domain->boxhi[1]-domain->boxlo[1])*
-                 (domain->boxhi[2]-domain->boxlo[2]) ,1.0/3.0);
-   
-  // save box sizes for coordinate rescaling
-  // calculate strains and asymmetry parameter
-  // alpha=lampdai[first]/lampbdai[second] for the two perp directions
-
-  alpha0 = 1.0;
-  for (int m = 0; m < 3; m++) {
-    domainloi[m] = *domainlo[m];
-    domainhii[m] = *domainhi[m];
-    lambdai[m] = (*domainhi[m] - *domainlo[m])/L;    
-    lambdaf[m] = ( m==dir ? lambda_final : 1.0/sqrt(lambda_final) ) ;    
-    if (m != dir) {
-      alpha0*= lambdai[m];
-      alpha0=1.0/alpha0;
-    }
-  } 
-  
-  if (comm->me == 0) {
-    if (screen) {
-      fprintf(screen,"Initial strain = %g %g %g\n",
-	      lambdai[0],lambdai[1],lambdai[2]);
-      fprintf(screen,"Target strain = %g %g %g\n",
-	      lambdaf[0],lambdaf[1],lambdaf[2]);
-    }
-    if (logfile) {
-      fprintf(logfile,"Initial strain = %g %g %g\n",
-	      lambdai[0],lambdai[1],lambdai[2]);
-      fprintf(logfile,"Target strain = %g %g %g\n",
-	      lambdaf[0],lambdaf[1],lambdaf[2]);
-    }
-  }
-
-  if (force->kspace) kspace_flag = 1;
-  else kspace_flag = 0;
-
-  // detect if any fix rigid exist so rigid bodies can be re-scaled
-  // rfix[] = indices to each fix rigid
-
-  delete [] rfix;
-  nrigid = 0;
-  rfix = NULL;
-
-  for (int i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"rigid") == 0 ||
-	strcmp(modify->fix[i]->style,"poems") == 0) nrigid++;
-  if (nrigid) {
-    rfix = new int[nrigid];
-    nrigid = 0;
-    for (int i = 0; i < modify->nfix; i++)
-      if (strcmp(modify->fix[i]->style,"rigid") == 0 ||
-	  strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixUniaxial::end_of_step()
-{
-  int i,m;
-  double oldlo,oldhi,newlo,newhi,ratio;
-
-  double delta = update->ntimestep - update->beginstep;
-  delta /= update->endstep - update->beginstep;
- 
-  double lvalue[3];
-
-  // linear interpolation of strain in specified direction
-  
-  lvalue[dir] = lambdai[dir]*(1.0-delta) + lambdaf[dir]*delta;
-
-  // linear interpolation of asymmetry parameter in the perp direction
-
-  double alpha = alpha0*(1-delta) + delta;
-
-  // calculate strains perpendicular to dir
-
-  for (m = 0; m < 3; m++)
-    if (m != dir) {
-      lvalue[m] = sqrt(alpha/lvalue[dir]);
-      alpha=1.0/alpha;
-    }
-
-  // apply appropriate rescaling in each dimension
-
-  double **x = atom->x;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  for (m = 0; m < 3; m++) {
-    oldlo = *domainlo[m];
-    oldhi = *domainhi[m];
-    
-    newlo = domainloi[m] * lvalue[m]/lambdai[m];
-    newhi = domainhii[m] * lvalue[m]/lambdai[m];
-    ratio = (newhi - newlo) / *domainprd[m];
-
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	x[i][m] = newlo + (x[i][m] - oldlo) * ratio;
-    
-    *domainlo[m] = newlo;
-    *domainhi[m] = newhi;
-    domain->prd[m] = *domainprd[m] = newhi - newlo;
-  
-    if (nrigid)
-      for (i = 0; i < nrigid; i++)
-	modify->fix[rfix[i]]->dilate(m,oldlo,oldhi,newlo,newhi);
-  }
-
-  // redo KSpace coeffs since volume has changed
-
-  if (kspace_flag) force->kspace->setup();
-}

src/fix_volume_rescale.cpp

@@ -1,194 +0,0 @@
-/* ----------------------------------------------------------------------
-   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_volume_rescale.h"
-#include "update.h"
-#include "force.h"
-#include "modify.h"
-#include "kspace.h"
-#include "domain.h"
-#include "atom.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-/* ---------------------------------------------------------------------- */
-
-FixVolRescale::FixVolRescale(LAMMPS *lmp, int narg, char **arg) :
-  Fix(lmp, narg, arg)
-{
-  if (narg < 7) error->all("Illegal fix volume/rescale command");
-
-  box_change = 1;
-
-  nevery = atoi(arg[3]);
-  if (nevery <= 0) error->all("Illegal fix volume/rescale command");
-
-  xflag = yflag = zflag = 0;
-  int iarg = 4;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"x") == 0) {
-      if (iarg+3 > narg) error->all("Illegal fix volume/rescale command");
-      if (domain->xperiodic == 0) 
-	error->all("Cannot fix volume/rescale on a non-periodic boundary");
-      xflag = 1;
-      xlo_stop = atof(arg[iarg+1]);
-      xhi_stop = atof(arg[iarg+2]);
-      iarg += 3;
-    } else if (strcmp(arg[iarg],"y") == 0) {
-      if (iarg+3 > narg) error->all("Illegal fix volume/rescale command");
-      if (domain->yperiodic == 0) 
-	error->all("Cannot fix volume/rescale on a non-periodic boundary");
-      yflag = 1;
-      ylo_stop = atof(arg[iarg+1]);
-      yhi_stop = atof(arg[iarg+2]);
-      iarg += 3;
-    } else if (strcmp(arg[iarg],"z") == 0) {
-      if (iarg+3 > narg) error->all("Illegal fix volume/rescale command");
-      if (domain->zperiodic == 0) 
-	error->all("Cannot fix volume/rescale on a non-periodic boundary");
-      zflag = 1;
-      zlo_stop = atof(arg[iarg+1]);
-      zhi_stop = atof(arg[iarg+2]);
-      iarg += 3;
-    } else error->all("Illegal fix volume/rescale command");
-  }
-
-  if (domain->triclinic)
-    error->all("Cannot use fix volume/rescale with triclinic box");
-
-  nrigid = 0;
-  rfix = NULL;
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixVolRescale::~FixVolRescale()
-{
-  delete [] rfix;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixVolRescale::setmask()
-{
-  int mask = 0;
-  mask |= END_OF_STEP;
-  return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixVolRescale::init()
-{
-  xlo_start = domain->boxlo[0];
-  xhi_start = domain->boxhi[0];
-  ylo_start = domain->boxlo[1];
-  yhi_start = domain->boxhi[1];
-  zlo_start = domain->boxlo[2];
-  zhi_start = domain->boxhi[2];
-
-  if (force->kspace) kspace_flag = 1;
-  else kspace_flag = 0;
-
-  // detect if any fix rigid exist so rigid bodies can be rescaled
-  // rfix[] = indices to each fix rigid
-
-  delete [] rfix;
-  nrigid = 0;
-  rfix = NULL;
-
-  for (int i = 0; i < modify->nfix; i++)
-    if (strcmp(modify->fix[i]->style,"rigid") == 0 ||
-	strcmp(modify->fix[i]->style,"poems") == 0) nrigid++;
-  if (nrigid) {
-    rfix = new int[nrigid];
-    nrigid = 0;
-    for (int i = 0; i < modify->nfix; i++)
-      if (strcmp(modify->fix[i]->style,"rigid") == 0 ||
-	  strcmp(modify->fix[i]->style,"poems") == 0) rfix[nrigid++] = i;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixVolRescale::end_of_step()
-{
-  int i;
-  double oldlo,oldhi,newlo,newhi,ratio;
-
-  double delta = update->ntimestep - update->beginstep;
-  delta /= update->endstep - update->beginstep;
-
-  double **x = atom->x;
-  int *mask = atom->mask;
-  int nlocal = atom->nlocal;
-
-  if (xflag) {
-    oldlo = domain->boxlo[0];
-    oldhi = domain->boxhi[0];
-    newlo = xlo_start + delta * (xlo_stop-xlo_start);
-    newhi = xhi_start + delta * (xhi_stop-xhi_start);
-    ratio = (newhi - newlo) / domain->xprd;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	x[i][0] = newlo + (x[i][0] - oldlo) * ratio;
-    domain->boxlo[0] = newlo;
-    domain->boxhi[0] = newhi;
-    domain->prd[0] = domain->xprd = newhi - newlo;
-    if (nrigid)
-      for (i = 0; i < nrigid; i++)
-	modify->fix[rfix[i]]->dilate(0,oldlo,oldhi,newlo,newhi);
-  }
-
-  if (yflag) {
-    oldlo = domain->boxlo[1];
-    oldhi = domain->boxhi[1];
-    newlo = ylo_start + delta * (ylo_stop-ylo_start);
-    newhi = yhi_start + delta * (yhi_stop-yhi_start);
-    ratio = (newhi - newlo) / domain->yprd;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	x[i][1] = newlo + (x[i][1] - oldlo) * ratio;
-    domain->boxlo[1] = newlo;
-    domain->boxhi[1] = newhi;
-    domain->prd[1] = domain->yprd = newhi - newlo;
-    if (nrigid)
-      for (i = 0; i < nrigid; i++)
-	modify->fix[rfix[i]]->dilate(1,oldlo,oldhi,newlo,newhi);
-  }
-
-  if (zflag) {
-    oldlo = domain->boxlo[2];
-    oldhi = domain->boxhi[2];
-    newlo = zlo_start + delta * (zlo_stop-zlo_start);
-    newhi = zhi_start + delta * (zhi_stop-zhi_start);
-    ratio = (newhi - newlo) / domain->zprd;
-    for (i = 0; i < nlocal; i++)
-      if (mask[i] & groupbit)
-	x[i][2] = newlo + (x[i][2] - oldlo) * ratio;
-    domain->boxlo[2] = newlo;
-    domain->boxhi[2] = newhi;
-    domain->prd[2] = domain->zprd = newhi - newlo;
-    if (nrigid)
-      for (i = 0; i < nrigid; i++)
-	  modify->fix[rfix[i]]->dilate(2,oldlo,oldhi,newlo,newhi);
-  }
-
-  // redo KSpace coeffs since volume has changed
-
-  if (kspace_flag) force->kspace->setup();
-}

src/math_extra.h

@@ -0,0 +1,435 @@
+/* ----------------------------------------------------------------------
+   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)
+------------------------------------------------------------------------- */
+
+#ifndef MATH_EXTRA_H
+#define MATH_EXTRA_H
+
+#include "math.h"
+#include "stdio.h"
+#include "string.h"
+#include "error.h"
+
+namespace MathExtra {
+
+  // 3 vector operations
+
+  inline void normalize3(const double *v, double *ans);
+  inline double dot3(const double *v1, const double *v2);
+  inline void cross3(const double *v1, const double *v2, double *ans);
+
+  // 3x3 matrix operations
+
+  inline double det3(const double mat[3][3]);
+  inline void diag_times3(const double *diagonal, const double mat[3][3],
+                          double ans[3][3]);
+  inline void plus3(const double m[3][3], const double m2[3][3],
+                    double ans[3][3]);
+  inline void transpose_times3(const double mat1[3][3], 
+			       const double mat2[3][3],
+                               double ans[3][3]);
+  inline void invert3(const double mat[3][3], double ans[3][3]);
+  inline void times_column3(const double mat[3][3], const double*vec,
+                            double *ans);
+  inline void transpose_times_column3(const double mat[3][3],const double*vec,
+                                      double *ans);
+  inline void row_times3(const double *v, const double m[3][3],
+                         double *ans);
+  inline void mldivide3(const double mat[3][3], const double *vec,
+                        double *ans, LAMMPS_NS::Error *error);
+  inline void write3(const double mat[3][3]);
+      
+  // quaternion operations
+  
+  inline void normalize4(double *quat);
+  inline void quat_to_mat(const double *quat, double mat[3][3]);
+  inline void quat_to_mat_trans(const double *quat, double mat[3][3]);
+  inline void axisangle_to_quat(const double *v, const double angle,
+                                double *quat);
+  inline void multiply_quat_quat(const double *one, const double *two,
+				 double *ans);
+  inline void multiply_vec_quat(const double *one, const double *two,
+				double *ans);
+
+  // shape matrix operations
+
+  inline void multiply_shape_shape(const double *one, const double *two,
+				   double *ans);
+};
+
+/* ----------------------------------------------------------------------
+   normalize a vector
+------------------------------------------------------------------------- */
+
+void MathExtra::normalize3(const double *v, double *ans)
+{
+  double den = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+  ans[0] = v[0]/den;
+  ans[1] = v[1]/den;
+  ans[2] = v[2]/den;
+}
+
+/* ----------------------------------------------------------------------
+   dot product of 2 vectors
+------------------------------------------------------------------------- */
+
+double MathExtra::dot3(const double *v1, const double *v2)
+{
+  return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
+}
+
+/* ----------------------------------------------------------------------
+   cross product of 2 vectors
+------------------------------------------------------------------------- */
+
+void MathExtra::cross3(const double *v1, const double *v2, double *ans)
+{
+  ans[0] = v1[1]*v2[2]-v1[2]*v2[1];
+  ans[1] = v1[2]*v2[0]-v1[0]*v2[2];
+  ans[2] = v1[0]*v2[1]-v1[1]*v2[0];
+}
+
+/* ----------------------------------------------------------------------
+   determinant of a matrix
+------------------------------------------------------------------------- */
+
+double MathExtra::det3(const double m[3][3])
+{
+  double ans = m[0][0]*m[1][1]*m[2][2] - m[0][0]*m[1][2]*m[2][1] - 
+    m[1][0]*m[0][1]*m[2][2] + m[1][0]*m[0][2]*m[2][1] + 
+    m[2][0]*m[0][1]*m[1][2] - m[2][0]*m[0][2]*m[1][1];
+  return ans;
+}
+
+/* ----------------------------------------------------------------------
+   diagonal matrix times a full matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::diag_times3(const double *d, const double m[3][3],
+			    double ans[3][3])
+{
+  ans[0][0] = d[0]*m[0][0];
+  ans[0][1] = d[0]*m[0][1];
+  ans[0][2] = d[0]*m[0][2];
+  ans[1][0] = d[1]*m[1][0];
+  ans[1][1] = d[1]*m[1][1];
+  ans[1][2] = d[1]*m[1][2];
+  ans[2][0] = d[2]*m[2][0];
+  ans[2][1] = d[2]*m[2][1];
+  ans[2][2] = d[2]*m[2][2];
+}
+
+/* ----------------------------------------------------------------------
+   add two matrices
+------------------------------------------------------------------------- */
+
+void MathExtra::plus3(const double m[3][3], const double m2[3][3],
+		      double ans[3][3])
+{
+  ans[0][0] = m[0][0]+m2[0][0];
+  ans[0][1] = m[0][1]+m2[0][1];
+  ans[0][2] = m[0][2]+m2[0][2];
+  ans[1][0] = m[1][0]+m2[1][0];
+  ans[1][1] = m[1][1]+m2[1][1];
+  ans[1][2] = m[1][2]+m2[1][2];
+  ans[2][0] = m[2][0]+m2[2][0];
+  ans[2][1] = m[2][1]+m2[2][1];
+  ans[2][2] = m[2][2]+m2[2][2];
+}
+
+/* ----------------------------------------------------------------------
+   multiply the transpose of mat1 times mat2
+------------------------------------------------------------------------- */
+
+void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3],
+				 double ans[3][3])
+{
+  ans[0][0] = m[0][0]*m2[0][0]+m[1][0]*m2[1][0]+m[2][0]*m2[2][0];
+  ans[0][1] = m[0][0]*m2[0][1]+m[1][0]*m2[1][1]+m[2][0]*m2[2][1];
+  ans[0][2] = m[0][0]*m2[0][2]+m[1][0]*m2[1][2]+m[2][0]*m2[2][2];
+  ans[1][0] = m[0][1]*m2[0][0]+m[1][1]*m2[1][0]+m[2][1]*m2[2][0];
+  ans[1][1] = m[0][1]*m2[0][1]+m[1][1]*m2[1][1]+m[2][1]*m2[2][1];
+  ans[1][2] = m[0][1]*m2[0][2]+m[1][1]*m2[1][2]+m[2][1]*m2[2][2];
+  ans[2][0] = m[0][2]*m2[0][0]+m[1][2]*m2[1][0]+m[2][2]*m2[2][0];
+  ans[2][1] = m[0][2]*m2[0][1]+m[1][2]*m2[1][1]+m[2][2]*m2[2][1];
+  ans[2][2] = m[0][2]*m2[0][2]+m[1][2]*m2[1][2]+m[2][2]*m2[2][2];
+}
+
+/* ----------------------------------------------------------------------
+   invert a matrix
+   does NOT checks for singular or badly scaled matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::invert3(const double m[3][3], double ans[3][3])
+{
+  double den = m[0][0]*m[1][1]*m[2][2]-m[0][0]*m[1][2]*m[2][1];
+  den += -m[1][0]*m[0][1]*m[2][2]+m[1][0]*m[0][2]*m[2][1];
+  den += m[2][0]*m[0][1]*m[1][2]-m[2][0]*m[0][2]*m[1][1];
+
+  ans[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1]) / den;
+  ans[0][1] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]) / den;
+  ans[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1]) / den;
+  ans[1][0] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]) / den;
+  ans[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0]) / den;
+  ans[1][2] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]) / den;
+  ans[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0]) / den;
+  ans[2][1] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]) / den;
+  ans[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0]) / den;
+}
+
+/* ----------------------------------------------------------------------
+   matrix times vector
+------------------------------------------------------------------------- */
+
+void MathExtra::times_column3(const double m[3][3], const double *v,
+			      double *ans) 
+{
+  ans[0] = m[0][0]*v[0]+m[0][1]*v[1]+m[0][2]*v[2];
+  ans[1] = m[1][0]*v[0]+m[1][1]*v[1]+m[1][2]*v[2];
+  ans[2] = m[2][0]*v[0]+m[2][1]*v[1]+m[2][2]*v[2];
+}
+
+/* ----------------------------------------------------------------------
+   transposed matrix times vector
+------------------------------------------------------------------------- */
+
+void MathExtra::transpose_times_column3(const double m[3][3], const double *v,
+					double *ans)
+{
+  ans[0] = m[0][0]*v[0]+m[1][0]*v[1]+m[2][0]*v[2];
+  ans[1] = m[0][1]*v[0]+m[1][1]*v[1]+m[2][1]*v[2];
+  ans[2] = m[0][2]*v[0]+m[1][2]*v[1]+m[2][2]*v[2];
+}
+
+/* ----------------------------------------------------------------------
+   row vector times matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::row_times3(const double *v, const double m[3][3],
+			   double *ans)
+{
+  ans[0] = m[0][0]*v[0]+v[1]*m[1][0]+v[2]*m[2][0];
+  ans[1] = v[0]*m[0][1]+m[1][1]*v[1]+v[2]*m[2][1];
+  ans[2] = v[0]*m[0][2]+v[1]*m[1][2]+m[2][2]*v[2];
+}
+
+/* ----------------------------------------------------------------------
+   solve Ax = b or M ans = v
+   use gaussian elimination & partial pivoting on matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::mldivide3(const double m[3][3], const double *v, double *ans,
+			  LAMMPS_NS::Error *error)
+{
+  // create augmented matrix for pivoting
+
+  double aug[3][4];
+  for (unsigned i = 0; i < 3; i++) {
+    aug[i][3] = v[i];
+    for (unsigned j = 0; j < 3; j++) aug[i][j] = m[i][j];
+  }
+
+  for (unsigned i = 0; i < 2; i++) {
+    unsigned p = i;
+    for (unsigned j = i+1; j < 3; j++) {
+      if (fabs(aug[j][i]) > fabs(aug[i][i])) {
+        double tempv[4];
+        memcpy(tempv,aug[i],4*sizeof(double));
+        memcpy(aug[i],aug[j],4*sizeof(double));
+        memcpy(aug[j],tempv,4*sizeof(double));
+      }
+    }
+
+    while (aug[p][i] == 0 && p < 3) p++;
+
+    if (p == 3) error->all("Bad matrix inversion in MathExtra::mldivide3");
+    else
+      if (p != i) {
+        double tempv[4];
+        memcpy(tempv,aug[i],4*sizeof(double));
+        memcpy(aug[i],aug[p],4*sizeof(double));
+        memcpy(aug[p],tempv,4*sizeof(double));
+      }
+
+    for (unsigned j = i+1; j < 3; j++) {
+      double m = aug[j][i]/aug[i][i];
+      for (unsigned k=i+1; k<4; k++) aug[j][k]-=m*aug[i][k];
+    }
+  }
+
+  if (aug[2][2] == 0)
+    error->all("Bad matrix inversion in MathExtra::mldivide3");
+  
+  // back substitution
+
+  ans[2] = aug[2][3]/aug[2][2];
+  for (int i = 1; i >= 0; i--) {
+    double sumax = 0.0;
+    for (unsigned j = i+1; j < 3; j++) sumax += aug[i][j]*ans[j];
+    ans[i] = (aug[i][3]-sumax) / aug[i][i];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   output a matrix
+------------------------------------------------------------------------- */
+
+void MathExtra::write3(const double mat[3][3])
+{
+  for (unsigned i = 0; i < 3; i++) {
+    for (unsigned j = 0; j < 3; j++) printf("%g ",mat[i][j]);
+    printf("\n");
+  }
+}
+
+/* ----------------------------------------------------------------------
+   normalize a quaternion
+------------------------------------------------------------------------- */
+
+void MathExtra::normalize4(double *quat)
+{
+  double den = sqrt(quat[0]*quat[0]+quat[1]*quat[1] +
+		    quat[2]*quat[2]+quat[3]*quat[3]);
+  quat[0] /= den;
+  quat[1] /= den;
+  quat[2] /= den;
+  quat[3] /= den;
+}
+
+/* ----------------------------------------------------------------------
+   compute rotation matrix from quaternion
+   quat = [w i j k]
+------------------------------------------------------------------------- */
+
+void MathExtra::quat_to_mat(const double *quat, double mat[3][3])
+{
+  double w2 = quat[0]*quat[0];
+  double i2 = quat[1]*quat[1];
+  double j2 = quat[2]*quat[2];
+  double k2 = quat[3]*quat[3];
+  double twoij = 2.0*quat[1]*quat[2];
+  double twoik = 2.0*quat[1]*quat[3];
+  double twojk = 2.0*quat[2]*quat[3];
+  double twoiw = 2.0*quat[1]*quat[0];
+  double twojw = 2.0*quat[2]*quat[0];
+  double twokw = 2.0*quat[3]*quat[0];
+
+  mat[0][0] = w2+i2-j2-k2;
+  mat[0][1] = twoij-twokw;
+  mat[0][2] = twojw+twoik;
+
+  mat[1][0] = twoij+twokw;
+  mat[1][1] = w2-i2+j2-k2;
+  mat[1][2] = twojk-twoiw;
+	
+  mat[2][0] = twoik-twojw;
+  mat[2][1] = twojk+twoiw;
+  mat[2][2] = w2-i2-j2+k2;
+}
+
+/* ----------------------------------------------------------------------
+   compute rotation matrix from quaternion conjugate
+   quat = [w i j k]
+------------------------------------------------------------------------- */
+
+void MathExtra::quat_to_mat_trans(const double *quat, double mat[3][3])
+{
+  double w2 = quat[0]*quat[0];
+  double i2 = quat[1]*quat[1];
+  double j2 = quat[2]*quat[2];
+  double k2 = quat[3]*quat[3];
+  double twoij = 2.0*quat[1]*quat[2];
+  double twoik = 2.0*quat[1]*quat[3];
+  double twojk = 2.0*quat[2]*quat[3];
+  double twoiw = 2.0*quat[1]*quat[0];
+  double twojw = 2.0*quat[2]*quat[0];
+  double twokw = 2.0*quat[3]*quat[0];
+
+  mat[0][0] = w2+i2-j2-k2;
+  mat[1][0] = twoij-twokw;
+  mat[2][0] = twojw+twoik;
+
+  mat[0][1] = twoij+twokw;
+  mat[1][1] = w2-i2+j2-k2;
+  mat[2][1] = twojk-twoiw;
+	
+  mat[0][2] = twoik-twojw;
+  mat[1][2] = twojk+twoiw;
+  mat[2][2] = w2-i2-j2+k2;
+}
+
+/* ----------------------------------------------------------------------
+   compute quaternion from axis-angle rotation
+   v MUST be a unit vector
+------------------------------------------------------------------------- */
+
+void MathExtra::axisangle_to_quat(const double *v, const double angle,
+				  double *quat)
+{
+  double halfa = 0.5*angle;
+  double sina = sin(halfa);
+  quat[0] = cos(halfa);
+  quat[1] = v[0]*sina;
+  quat[2] = v[1]*sina;
+  quat[3] = v[2]*sina;
+}
+
+/* ----------------------------------------------------------------------
+   multiply 2 quaternions
+   effectively concatenates rotations
+   NOT a commutative operation
+------------------------------------------------------------------------- */
+
+void MathExtra::multiply_quat_quat(const double *one, const double *two,
+				   double *ans)
+{
+  ans[0] = one[0]*two[0]-one[1]*two[1]-one[2]*two[2]-one[3]*two[3];
+  ans[1] = one[0]*two[1]+one[1]*two[0]+one[2]*two[3]-one[3]*two[2];
+  ans[2] = one[0]*two[2]-one[1]*two[3]+one[2]*two[0]+one[3]*two[1];
+  ans[3] = one[0]*two[3]+one[1]*two[2]-one[2]*two[1]+one[3]*two[0];
+}
+
+/* ----------------------------------------------------------------------
+   multiply 3 vector times quaternion
+   3 vector one is treated as quaternion (0,one)
+------------------------------------------------------------------------- */
+
+void MathExtra::multiply_vec_quat(const double *one, const double *two,
+				  double *ans)
+{
+  ans[0] = -one[0]*two[1]-one[1]*two[2]-one[2]*two[3];
+  ans[1] =  one[0]*two[0]+one[1]*two[3]-one[2]*two[2];
+  ans[2] = -one[0]*two[3]+one[1]*two[0]+one[2]*two[1];
+  ans[3] =  one[0]*two[2]-one[1]*two[1]+one[2]*two[0];
+}
+
+/* ----------------------------------------------------------------------
+   multiply 2 shape matrices
+   upper-triangular 3x3, stored as 6-vector in Voigt notation
+------------------------------------------------------------------------- */
+
+void MathExtra::multiply_shape_shape(const double *one, const double *two,
+				     double *ans)
+{
+  ans[0] = one[0]*two[0];
+  ans[1] = one[1]*two[1];
+  ans[2] = one[2]*two[2];
+  ans[3] = one[1]*two[3] + one[3]*two[2];
+  ans[4] = one[0]*two[4] + one[5]*two[3] + one[4]*two[2];
+  ans[5] = one[0]*two[5] + one[5]*two[1];
+}
+
+#endif

src/style_class2.h

@@ -1,56 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 AngleInclude
-#include "angle_class2.h"
-#endif
-
-#ifdef AngleClass
-AngleStyle(class2,AngleClass2)
-#endif
-
-#ifdef BondInclude
-#include "bond_class2.h"
-#endif
-
-#ifdef BondClass
-BondStyle(class2,BondClass2)
-#endif
-
-#ifdef DihedralInclude
-#include "dihedral_class2.h"
-#endif
-
-#ifdef DihedralClass
-DihedralStyle(class2,DihedralClass2)
-#endif
-
-#ifdef ImproperInclude
-#include "improper_class2.h"
-#endif
-
-#ifdef ImproperClass
-ImproperStyle(class2,ImproperClass2)
-#endif
-
-#ifdef PairInclude
-#include "pair_lj_class2.h"
-#include "pair_lj_class2_coul_cut.h"
-#include "pair_lj_class2_coul_long.h"
-#endif
-
-#ifdef PairClass
-PairStyle(lj/class2,PairLJClass2)
-PairStyle(lj/class2/coul/cut,PairLJClass2CoulCut)
-PairStyle(lj/class2/coul/long,PairLJClass2CoulLong)
-#endif

src/style_dpd.h

@@ -1,28 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 AtomInclude
-#include "atom_vec_dpd.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(dpd,AtomVecDPD)
-#endif
-
-#ifdef PairInclude
-#include "pair_dpd.h"
-#endif
-
-#ifdef PairClass
-PairStyle(dpd,PairDPD)
-#endif

src/style_granular.h

@@ -1,50 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 AtomInclude
-#include "atom_vec_granular.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(granular,AtomVecGranular)
-# endif
-
-#ifdef FixInclude
-#include "fix_freeze.h"
-#include "fix_gran_diag.h"
-#include "fix_nve_gran.h"
-#include "fix_pour.h"
-#include "fix_shear_history.h"
-#include "fix_wall_gran.h"
-#endif
-
-#ifdef FixClass
-FixStyle(freeze,FixFreeze)
-FixStyle(gran/diag,FixGranDiag)
-FixStyle(nve/gran,FixNVEGran)
-FixStyle(pour,FixPour)
-FixStyle(SHEAR_HISTORY,FixShearHistory)
-FixStyle(wall/gran,FixWallGran)
-#endif
-
-#ifdef PairInclude
-#include "pair_gran_hertzian.h"
-#include "pair_gran_history.h"
-#include "pair_gran_no_history.h"
-#endif
-
-#ifdef PairClass
-PairStyle(gran/hertzian,PairGranHertzian)
-PairStyle(gran/history,PairGranHistory)
-PairStyle(gran/no_history,PairGranNoHistory)
-#endif

src/style_kspace.h

@@ -1,38 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 KSpaceInclude
-#include "ewald.h"
-#include "pppm.h"
-#include "pppm_tip4p.h"
-#endif
-
-#ifdef KSpaceClass
-KSpaceStyle(ewald,Ewald)
-KSpaceStyle(pppm,PPPM)
-KSpaceStyle(pppm/tip4p,PPPMTIP4P)
-#endif
-
-#ifdef PairInclude
-#include "pair_buck_coul_long.h"
-#include "pair_lj_cut_coul_long.h"
-#include "pair_lj_cut_coul_long_tip4p.h"
-#include "pair_lj_charmm_coul_long.h"
-#endif
-
-#ifdef PairClass
-PairStyle(buck/coul/long,PairBuckCoulLong)
-PairStyle(lj/cut/coul/long,PairLJCutCoulLong)
-PairStyle(lj/cut/coul/long/tip4p,PairLJCutCoulLongTIP4P)
-PairStyle(lj/charmm/coul/long,PairLJCharmmCoulLong)
-#endif

src/style_manybody.h

@@ -1,30 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 PairInclude
-#include "pair_airebo.h"
-#include "pair_eam.h"
-#include "pair_eam_alloy.h"
-#include "pair_eam_fs.h"
-#include "pair_sw.h"
-#include "pair_tersoff.h"
-#endif
-
-#ifdef PairClass
-PairStyle(airebo,PairAIREBO)
-PairStyle(eam,PairEAM)
-PairStyle(eam/alloy,PairEAMAlloy)
-PairStyle(eam/fs,PairEAMFS)
-PairStyle(sw,PairSW)
-PairStyle(tersoff,PairTersoff)
-#endif

src/style_molecule.h

@@ -1,116 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 AngleInclude
-#include "angle_charmm.h"
-#include "angle_cosine.h"
-#include "angle_cosine_squared.h"
-#include "angle_harmonic.h"
-#include "angle_hybrid.h"
-#endif
-
-#ifdef AngleClass
-AngleStyle(charmm,AngleCharmm)
-AngleStyle(cosine,AngleCosine)
-AngleStyle(cosine/squared,AngleCosineSquared)
-AngleStyle(harmonic,AngleHarmonic)
-AngleStyle(hybrid,AngleHybrid)
-#endif
-
-#ifdef AtomInclude
-#include "atom_vec_angle.h"
-#include "atom_vec_bond.h"
-#include "atom_vec_full.h"
-#include "atom_vec_molecular.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(angle,AtomVecAngle)
-AtomStyle(bond,AtomVecBond)
-AtomStyle(full,AtomVecFull)
-AtomStyle(molecular,AtomVecMolecular)
-#endif
-
-#ifdef BondInclude
-#include "bond_fene.h"
-#include "bond_fene_expand.h"
-#include "bond_harmonic.h"
-#include "bond_hybrid.h"
-#include "bond_morse.h"
-#include "bond_nonlinear.h"
-#include "bond_quartic.h"
-#endif
-
-#ifdef BondClass
-BondStyle(fene,BondFENE)
-BondStyle(fene/expand,BondFENEExpand)
-BondStyle(harmonic,BondHarmonic)
-BondStyle(hybrid,BondHybrid)
-BondStyle(morse,BondMorse)
-BondStyle(nonlinear,BondNonlinear)
-BondStyle(quartic,BondQuartic)
-#endif
-
-#ifdef DihedralInclude
-#include "dihedral_charmm.h"
-#include "dihedral_harmonic.h"
-#include "dihedral_helix.h"
-#include "dihedral_hybrid.h"
-#include "dihedral_multi_harmonic.h"
-#include "dihedral_opls.h"
-#endif
-
-#ifdef DihedralClass
-DihedralStyle(charmm,DihedralCharmm)
-DihedralStyle(harmonic,DihedralHarmonic)
-DihedralStyle(helix,DihedralHelix)
-DihedralStyle(hybrid,DihedralHybrid)
-DihedralStyle(multi/harmonic,DihedralMultiHarmonic)
-DihedralStyle(opls,DihedralOPLS)
-#endif
-
-#ifdef DumpInclude
-#include "dump_bond.h"
-#endif
-
-#ifdef DumpClass
-DumpStyle(bond,DumpBond)
-#endif
-
-#ifdef FixInclude
-#endif
-
-#ifdef FixClass
-#endif
-
-#ifdef ImproperInclude
-#include "improper_cvff.h"
-#include "improper_harmonic.h"
-#include "improper_hybrid.h"
-#endif
-
-#ifdef ImproperClass
-ImproperStyle(cvff,ImproperCvff)
-ImproperStyle(harmonic,ImproperHarmonic)
-ImproperStyle(hybrid,ImproperHybrid)
-#endif
-
-#ifdef PairInclude
-#include "pair_lj_charmm_coul_charmm.h"
-#include "pair_lj_charmm_coul_charmm_implicit.h"
-#endif
-
-#ifdef PairClass
-PairStyle(lj/charmm/coul/charmm,PairLJCharmmCoulCharmm)
-PairStyle(lj/charmm/coul/charmm/implicit,PairLJCharmmCoulCharmmImplicit)
-#endif

src/style_opt.h

@@ -1,30 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 PairInclude
-#include "pair_eam_opt.h"
-#include "pair_eam_alloy_opt.h"
-#include "pair_eam_fs_opt.h"
-#include "pair_lj_charmm_coul_long_opt.h"
-#include "pair_lj_cut_opt.h"
-#include "pair_morse_opt.h"
-#endif
-
-#ifdef PairClass
-PairStyle(eam/opt,PairEAMOpt)
-PairStyle(eam/alloy/opt,PairEAMAlloyOpt)
-PairStyle(eam/fs/opt,PairEAMFSOpt)
-PairStyle(lj/cut/opt,PairLJCutOpt)
-PairStyle(lj/charmm/coul/long/opt,PairLJCharmmCoulLongOpt)
-PairStyle(morse/opt,PairMorseOpt)
-#endif

src/style_xtc.h

@@ -1,20 +0,0 @@
-/* ----------------------------------------------------------------------
-   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 DumpInclude
-#include "dump_xtc.h"
-#endif
-
-#ifdef DumpClass
-DumpStyle(xtc,DumpXTC)
-#endif