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

src/fix_addforce.cpp

@@ -0,0 +1,175 @@
+/* ----------------------------------------------------------------------
+   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 "fix_addforce.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "region.h"
+#include "respa.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixAddForce::FixAddForce(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 6) error->all("Illegal fix addforce command");
+
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  global_freq = 1;
+  extscalar = 1;
+  extvector = 1;
+
+  xvalue = atof(arg[3]);
+  yvalue = atof(arg[4]);
+  zvalue = atof(arg[5]);
+
+  // optional args
+
+  iregion = -1;
+
+  int iarg = 6;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"region") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix addforce command");
+      iregion = domain->find_region(arg[iarg+1]);
+      if (iregion == -1) error->all("Fix addforce region ID does not exist");
+      iarg += 2;
+    } else error->all("Illegal fix addforce command");
+  }
+
+  force_flag = 0;
+  foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAddForce::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+  mask |= THERMO_ENERGY;
+  mask |= POST_FORCE_RESPA;
+  mask |= MIN_POST_FORCE;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::init()
+{
+  if (strcmp(update->integrate_style,"respa") == 0)
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::setup(int vflag)
+{
+  if (strcmp(update->integrate_style,"verlet") == 0)
+    post_force(vflag);
+  else {
+    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+    post_force_respa(vflag,nlevels_respa-1,0);
+    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::min_setup(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::post_force(int vflag)
+{
+  double **x = atom->x;
+  double **f = atom->f;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
+  force_flag = 0;
+
+  // foriginal[0] = - x dot f = "potential" for added force
+  // foriginal[123] = force on atoms before extra force added
+
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (iregion >= 0 && 
+          !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
+	continue;
+
+      foriginal[0] -= xvalue*x[i][0] + yvalue*x[i][1] + zvalue*x[i][2];
+      foriginal[1] += f[i][0];
+      foriginal[2] += f[i][1];
+      foriginal[3] += f[i][2];
+      f[i][0] += xvalue;
+      f[i][1] += yvalue;
+      f[i][2] += zvalue;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::post_force_respa(int vflag, int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAddForce::min_post_force(int vflag)
+{
+  post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   potential energy of added force
+------------------------------------------------------------------------- */
+
+double FixAddForce::compute_scalar()
+{
+  // only sum across procs one time
+
+  if (force_flag == 0) {
+    MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world);
+    force_flag = 1;
+  }
+  return foriginal_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   return components of total force on fix group before force was changed
+------------------------------------------------------------------------- */
+
+double FixAddForce::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (force_flag == 0) {
+    MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world);
+    force_flag = 1;
+  }
+  return foriginal_all[n+1];
+}