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

src/compute_bond_local.cpp

@@ -21,14 +21,17 @@
 #include "domain.h"
 #include "force.h"
 #include "bond.h"
+#include "math_extra.h"
+#include "comm.h"
 #include "memory.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
 
 #define DELTA 10000
+#define EPSILON 1.0e-12
 
-enum{DIST,ENG,FORCE};
+enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -42,6 +45,8 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Compute bond/local used when bonds are not allowed");
 
   local_flag = 1;
+  comm_forward = 3;
+
   nvalues = narg - 3;
   if (nvalues == 1) size_local_cols = 0;
   else size_local_cols = nvalues;
@@ -51,16 +56,26 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
   nvalues = 0;
   for (int iarg = 3; iarg < narg; iarg++) {
     if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
-    else if (strcmp(arg[iarg],"eng") == 0) bstyle[nvalues++] = ENG;
+    else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
     else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
+    else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB;
+    else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT;
+    else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS;
+    else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA;
+    else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB;
     else error->all(FLERR,"Invalid keyword in compute bond/local command");
   }
 
   // set singleflag if need to call bond->single()
+  // set velflag if compute any quantities based on velocities
 
   singleflag = 0;
-  for (int i = 0; i < nvalues; i++)
+  ghostvelflag = 0;
-    if (bstyle[i] != DIST) singleflag = 1;
+  for (int i = 0; i < nvalues; i++) {
+    if (bstyle[i] == ENGPOT || bstyle[i] == FORCE) singleflag = 1;
+    if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS ||
+        bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1;
+  }
 
   nmax = 0;
 }
@@ -81,9 +96,17 @@ void ComputeBondLocal::init()
   if (force->bond == NULL)
     error->all(FLERR,"No bond style is defined for compute bond/local");
 
+  // set ghostvelflag if need to acquire ghost atom velocities
+
+  if (velflag && !comm->ghost_velocity) ghostvelflag = 1;
+  else ghostvelflag = 0;
+
   // do initial memory allocation so that memory_usage() is correct
 
+  initflag = 1;
   ncount = compute_bonds(0);
+  initflag = 0;
+
   if (ncount > nmax) reallocate(ncount);
   size_local_rows = ncount;
 }
@@ -117,10 +140,22 @@ int ComputeBondLocal::compute_bonds(int flag)
 {
   int i,m,n,nb,atom1,atom2,imol,iatom,btype;
   tagint tagprev;
-  double delx,dely,delz,rsq;
+  double dx,dy,dz,rsq;
+  double mass1,mass2,masstotal,invmasstotal;
+  double xcm[3],vcm[3];
+  double delr1[3],delr2[3],delv1[3],delv2[3];
+  double r12[3],vpar1,vpar2;
+  double vvib,vrotsq;
+  double inertia,omegasq;
+  double mvv2e;
+  double engpot,engtrans,engvib,engrot,engtot,fbond;
   double *ptr;
 
   double **x = atom->x;
+  double **v = atom->v;
+  int *type = atom->type;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
   tagint *tag = atom->tag;
   int *num_bond = atom->num_bond;
   tagint **bond_atom = atom->bond_atom;
@@ -136,7 +171,12 @@ int ComputeBondLocal::compute_bonds(int flag)
   int molecular = atom->molecular;
 
   Bond *bond = force->bond;
-  double eng,fbond;
+
+  // communicate ghost velocities if needed
+
+  if (ghostvelflag && !initflag) comm->forward_comm_compute(this);
+
+  // loop over all atoms and their bonds
 
   m = n = 0;
   for (atom1 = 0; atom1 < nlocal; atom1++) {
@@ -164,16 +204,93 @@ int ComputeBondLocal::compute_bonds(int flag)
       if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue;
       if (btype == 0) continue;
 
-      if (flag) {
+      if (!flag) {
-        delx = x[atom1][0] - x[atom2][0];
+        m++;
-        dely = x[atom1][1] - x[atom2][1];
+        continue;
-        delz = x[atom1][2] - x[atom2][2];
+      }
-        domain->minimum_image(delx,dely,delz);
-        rsq = delx*delx + dely*dely + delz*delz;
 
-        if (singleflag && (btype > 0))
+      dx = x[atom1][0] - x[atom2][0];
-          eng = bond->single(btype,rsq,atom1,atom2,fbond);
+      dy = x[atom1][1] - x[atom2][1];
-        else eng = fbond = 0.0;
+      dz = x[atom1][2] - x[atom2][2];
+      domain->minimum_image(dx,dy,dz);
+      rsq = dx*dx + dy*dy + dz*dz;
+      
+      if (btype == 0) {
+        engpot = fbond = 0.0;
+        engvib = engrot = engtrans = omegasq = vvib = 0.0;
+      } else {
+
+        if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond);
+        
+        if (velflag) {
+          if (rmass) {
+            mass1 = rmass[atom1];
+            mass2 = rmass[atom2];
+          }
+          else {
+            mass1 = mass[type[atom1]];
+            mass2 = mass[type[atom2]];
+          }
+          masstotal = mass1+mass2;
+          invmasstotal = 1.0 / (masstotal);
+          xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal;
+          xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal;
+          xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal;
+          vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal;
+          vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal;
+          vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal;
+          
+          engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm);
+          
+          // r12 = unit bond vector from atom1 to atom2
+          
+          MathExtra::sub3(x[atom2],x[atom1],r12);
+          MathExtra::norm3(r12);
+          
+          // delr = vector from COM to each atom
+          // delv = velocity of each atom relative to COM
+          
+          MathExtra::sub3(x[atom1],xcm,delr1);
+          MathExtra::sub3(x[atom2],xcm,delr2);
+          MathExtra::sub3(v[atom1],vcm,delv1);
+          MathExtra::sub3(v[atom2],vcm,delv2);
+          
+          // vpar = component of delv parallel to bond vector
+          
+          vpar1 = MathExtra::dot3(delv1,r12);
+          vpar2 = MathExtra::dot3(delv2,r12);
+          engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2);
+          
+          // vvib = relative velocity of 2 atoms along bond direction
+          // vvib < 0 for 2 atoms moving towards each other
+          // vvib > 0 for 2 atoms moving apart
+          
+          vvib = vpar2 - vpar1;
+          
+          // vrotsq = tangential speed squared of atom1 only
+          // omegasq = omega squared, and is the same for atom1 and atom2
+          
+          inertia = mass1*MathExtra::lensq3(delr1) + 
+            mass2*MathExtra::lensq3(delr2);
+          vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1;
+          omegasq = vrotsq / MathExtra::lensq3(delr1);
+          
+          engrot = 0.5 * inertia * omegasq;
+          
+          // sanity check: engtotal = engtrans + engvib + engrot
+          
+          //engtot = 0.5 * (mass1*MathExtra::lensq3(v[atom1]) + 
+          //                mass2*MathExtra::lensq3(v[atom2]));
+          //if (fabs(engtot-engtrans-engvib-engrot) > EPSILON)
+          //  error->one(FLERR,"Sanity check on 3 energy components failed");
+          
+          // scale energies by units
+
+          mvv2e = force->mvv2e;
+          engtrans *= mvv2e;
+          engvib *= mvv2e;
+          engrot *= mvv2e;
+        }
 
         if (nvalues == 1) ptr = &vector[m];
         else ptr = array[m];
@@ -183,12 +300,27 @@ int ComputeBondLocal::compute_bonds(int flag)
           case DIST:
             ptr[n] = sqrt(rsq);
             break;
-          case ENG:
+          case ENGPOT:
-            ptr[n] = eng;
+            ptr[n] = engpot;
             break;
           case FORCE:
             ptr[n] = sqrt(rsq)*fbond;
             break;
+          case ENGVIB:
+            ptr[n] = engvib;
+            break;
+          case ENGROT:
+            ptr[n] = engrot;
+            break;
+          case ENGTRANS:
+            ptr[n] = engtrans;
+            break;
+          case OMEGA:
+            ptr[n] = sqrt(omegasq);
+            break;
+          case VELVIB:
+            ptr[n] = vvib;
+            break;
           }
         }
       }
@@ -202,6 +334,43 @@ int ComputeBondLocal::compute_bonds(int flag)
 
 /* ---------------------------------------------------------------------- */
 
+int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf,
+                                        int pbc_flag, int *pbc)
+{
+  int i,j,m;
+
+  double **v = atom->v;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = v[j][0];
+    buf[m++] = v[j][1];
+    buf[m++] = v[j][2];
+  }
+
+  return m;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  double **v = atom->v;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    v[i][0] = buf[m++];
+    v[i][1] = buf[m++];
+    v[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void ComputeBondLocal::reallocate(int n)
 {
   // grow vector or array and indices array

src/compute_bond_local.h

@@ -30,13 +30,15 @@ class ComputeBondLocal : public Compute {
   ~ComputeBondLocal();
   void init();
   void compute_local();
+  int pack_forward_comm(int, int *, double *, int, int *);
+  void unpack_forward_comm(int, int, double *);
   double memory_usage();
 
  private:
   int nvalues;
   int ncount;
   int *bstyle;
-  int singleflag;
+  int singleflag,velflag,ghostvelflag,initflag;
 
   int nmax;
 

src/dump_local.cpp

@@ -21,6 +21,7 @@
 #include "compute.h"
 #include "domain.h"
 #include "update.h"
+#include "input.h"
 #include "memory.h"
 #include "error.h"
 #include "force.h"
@@ -45,7 +46,18 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
   nevery = force->inumeric(FLERR,arg[3]);
   if (nevery <= 0) error->all(FLERR,"Illegal dump local command");
 
-  size_one = nfield = narg-5;
+  nfield = narg - 5;
+
+  // expand args if any have wildcard character "*"
+
+  int expand = 0;
+  char **earg;
+  nfield = input->expand_args(nfield,&arg[5],1,earg);
+
+  if (earg != &arg[5]) expand = 1;
+
+  // allocate field vectors
+
   pack_choice = new FnPtrPack[nfield];
   vtype = new int[nfield];
 
@@ -67,7 +79,8 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
 
   // process attributes
 
-  parse_fields(narg,arg);
+  parse_fields(nfield,earg);
+  size_one = nfield;
 
   // setup format strings
 
@@ -88,11 +101,11 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
   // setup column string
 
   int n = 0;
-  for (int iarg = 5; iarg < narg; iarg++) n += strlen(arg[iarg]) + 2;
+  for (int iarg = 0; iarg < nfield; iarg++) n += strlen(earg[iarg]) + 2;
   columns = new char[n];
   columns[0] = '\0';
-  for (int iarg = 5; iarg < narg; iarg++) {
+  for (int iarg = 0; iarg < nfield; iarg++) {
-    strcat(columns,arg[iarg]);
+    strcat(columns,earg[iarg]);
     strcat(columns," ");
   }
 
@@ -102,6 +115,13 @@ DumpLocal::DumpLocal(LAMMPS *lmp, int narg, char **arg) :
   n = strlen(str) + 1;
   label = new char[n];
   strcpy(label,str);
+
+  // if wildcard expansion occurred, free earg memory from exapnd_args()
+
+  if (expand) {
+    for (int i = 0; i < nfield; i++) delete [] earg[i];
+    memory->sfree(earg);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -376,8 +396,8 @@ void DumpLocal::parse_fields(int narg, char **arg)
   // customize by adding to if statement
 
   int i;
-  for (int iarg = 5; iarg < narg; iarg++) {
+  for (int iarg = 0; iarg < narg; iarg++) {
-    i = iarg-5;
+    i = iarg;
 
     if (strcmp(arg[iarg],"index") == 0) {
       pack_choice[i] = &DumpLocal::pack_index;

src/fix_ave_atom.cpp

@@ -36,7 +36,8 @@ enum{X,V,F,COMPUTE,FIX,VARIABLE};
 
 FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg),
-  nvalues(0), which(NULL), argindex(NULL), value2index(NULL), ids(NULL), array(NULL)
+  nvalues(0), which(NULL), argindex(NULL), value2index(NULL), 
+  ids(NULL), array(NULL)
 {
   if (narg < 7) error->all(FLERR,"Illegal fix ave/atom command");