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

src/fix_efield.cpp

@@ -15,6 +15,7 @@
    Contributing author: Christina Payne (Vanderbilt U)
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "string.h"
 #include "stdlib.h"
 #include "fix_efield.h"
@@ -22,11 +23,15 @@
 #include "update.h"
 #include "force.h"
 #include "respa.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
 #include "error.h"
-#include "math.h"
 
 using namespace LAMMPS_NS;
 
+enum{NONE,EQUAL,ATOM};
+
 /* ---------------------------------------------------------------------- */
 
 FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
@@ -34,10 +39,39 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 6) error->all("Illegal fix efield command");
 
-  double factor = force->qe2f;
-  ex = factor * atof(arg[3]);
-  ey = factor * atof(arg[4]);
-  ez = factor * atof(arg[5]);
+  efactor = force->qe2f;
+  xstr = ystr = zstr = NULL;
+
+  if (strstr(arg[3],"v_") == arg[3]) {
+    int n = strlen(&arg[3][2]) + 1;
+    xstr = new char[n];
+    strcpy(xstr,&arg[3][2]);
+  } else ex = efactor * atof(arg[3]);
+
+  if (strstr(arg[4],"v_") == arg[4]) {
+    int n = strlen(&arg[4][2]) + 1;
+    xstr = new char[n];
+    strcpy(xstr,&arg[4][2]);
+  } else ey = efactor * atof(arg[4]);
+
+  if (strstr(arg[5],"v_") == arg[5]) {
+    int n = strlen(&arg[5][2]) + 1;
+    xstr = new char[n];
+    strcpy(xstr,&arg[5][2]);
+  } else ez = efactor * atof(arg[5]);
+
+  maxatom = 0;
+  efield = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixEfield::~FixEfield()
+{
+  delete [] xstr;
+  delete [] ystr;
+  delete [] zstr;
+  memory->destroy_2d_double_array(efield);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -59,6 +93,36 @@ void FixEfield::init()
   if (atom->q == NULL)
     error->all("Must use charged atom style with fix efield");
 
+  // check variables
+
+  if (xstr) {
+    xvar = input->variable->find(xstr);
+    if (xvar < 0) error->all("Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(xvar)) xvarstyle = EQUAL;
+    else if (input->variable->atomstyle(xvar)) xvarstyle = ATOM;
+    else error->all("Variable for fix efield is invalid style");
+  } else xvarstyle = NONE;
+  if (ystr) {
+    yvar = input->variable->find(ystr);
+    if (yvar < 0) error->all("Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(yvar)) yvarstyle = EQUAL;
+    else if (input->variable->atomstyle(yvar)) yvarstyle = ATOM;
+    else error->all("Variable for fix efield is invalid style");
+  } else yvarstyle = NONE;
+  if (zstr) {
+    zvar = input->variable->find(zstr);
+    if (zvar < 0) error->all("Variable name for fix efield does not exist");
+    if (input->variable->equalstyle(zvar)) zvarstyle = EQUAL;
+    else if (input->variable->atomstyle(zvar)) zvarstyle = ATOM;
+    else error->all("Variable for fix efield is invalid style");
+  } else zvarstyle = NONE;
+
+  if (xvarstyle == ATOM || yvarstyle == ATOM || zvarstyle == ATOM)
+    varflag = ATOM;
+  else if (xvarstyle == EQUAL || yvarstyle == EQUAL || zvarstyle == EQUAL)
+    varflag = EQUAL;
+  else varflag = NONE;
+
   if (strcmp(update->integrate_style,"respa") == 0)
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
 }
@@ -87,12 +151,52 @@ void FixEfield::post_force(int vflag)
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  for (int i = 0; i < nlocal; i++)
-    if (mask[i] & groupbit) {
-      f[i][0] += q[i]*ex;
-      f[i][1] += q[i]*ey;
-      f[i][2] += q[i]*ez;
+  // reallocate efield array if necessary
+
+  if (varflag == ATOM && nlocal > maxatom) {
+    maxatom = atom->nmax;
+    memory->destroy_2d_double_array(efield);
+    efield = memory->create_2d_double_array(maxatom,3,"efield:efield");
+  }
+
+  if (varflag == NONE) {
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	f[i][0] += q[i]*ex;
+	f[i][1] += q[i]*ey;
+	f[i][2] += q[i]*ez;
+      }
+
+  } else {
+    if (xstr) {
+      if (xvarstyle == EQUAL)
+	ex = efactor * input->variable->compute_equal(xvar);
+      else 
+	input->variable->compute_atom(xvar,igroup,&efield[0][0],3,0);
     }
+    if (ystr) {
+      if (yvarstyle == EQUAL)
+	ey = efactor * input->variable->compute_equal(yvar);
+      else 
+	input->variable->compute_atom(yvar,igroup,&efield[0][1],3,0);
+    }
+    if (zstr) {
+      if (zvarstyle == EQUAL)
+	ez = efactor * input->variable->compute_equal(zvar);
+      else 
+	input->variable->compute_atom(zvar,igroup,&efield[0][2],3,0);
+    }
+
+    for (int i = 0; i < nlocal; i++)
+      if (mask[i] & groupbit) {
+	if (xvarstyle == ATOM) f[i][0] += q[i]*efield[i][0];
+	else f[i][0] += q[i]*ex;
+	if (yvarstyle == ATOM) f[i][1] += q[i]*efield[i][1];
+	else f[i][1] += q[i]*ey;
+	if (zvarstyle == ATOM) f[i][2] += q[i]*efield[i][2];
+	else f[i][2] += q[i]*ez;
+      }
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -101,3 +205,14 @@ void FixEfield::post_force_respa(int vflag, int ilevel, int iloop)
 {
   if (ilevel == nlevels_respa-1) post_force(vflag);
 }
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+double FixEfield::memory_usage()
+{
+  double bytes = 0.0;
+  if (varflag == ATOM) bytes = atom->nmax*3 * sizeof(double);
+  return bytes;
+}

src/fix_efield.h

@@ -27,15 +27,24 @@ namespace LAMMPS_NS {
 class FixEfield : public Fix {
  public:
   FixEfield(class LAMMPS *, int, char **);
+  ~FixEfield();
   int setmask();
   void init();
   void setup(int);
   void post_force(int);
   void post_force_respa(int, int, int);
+  double memory_usage();
 
  private:
   double ex,ey,ez;
+  int varflag;
+  char *xstr,*ystr,*zstr;
+  int xvar,yvar,zvar,xvarstyle,yvarstyle,zvarstyle;
   int nlevels_respa;
+  double efactor;
+
+  int maxatom;
+  double **efield;
 };
 
 }

src/fix_move.cpp

@@ -791,6 +791,8 @@ void FixMove::final_integrate_respa(int ilevel, int iloop)
 double FixMove::memory_usage()
 {
   double bytes = atom->nmax*3 * sizeof(double);
+  if (displaceflag) bytes += atom->nmax*3 * sizeof(double);
+  if (velocityflag) bytes += atom->nmax*3 * sizeof(double);
   return bytes;
 }