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

src/compute_ti.cpp

@@ -0,0 +1,180 @@
+/* ----------------------------------------------------------------------
+   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: Sai Jayaraman (University of Notre Dame)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "string.h"
+#include "compute_ti.h"
+#include "update.h"
+#include "modify.h"
+#include "domain.h"
+#include "force.h"
+#include "pair.h"
+#include "kspace.h"
+#include "input.h"
+#include "variable.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{PAIR,TAIL,KSPACE};
+
+#define INVOKED_SCALAR 1
+
+/* ---------------------------------------------------------------------- */
+
+ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : 
+  Compute(lmp, narg, arg) 
+{ 
+  if (narg < 4) error->all("Illegal compute ti command");
+
+  peflag = 1;
+  scalar_flag = 1;
+  extscalar = 1;
+  timeflag = 1;
+
+  // terms come in triplets
+
+  nterms = (narg-3) / 3;
+  if (narg != 3*nterms + 3) error->all("Illegal compute ti command");
+
+  which = new int[nterms];
+  ivar1 = new int[nterms];
+  ivar2 = new int[nterms];
+  var1 = new char*[nterms];
+  var2 = new char*[nterms];
+  pptr = new Pair*[nterms];
+  pstyle = new char*[nterms];
+
+  for (int m = 0; m < nterms; m++) pstyle[m] = NULL;
+    
+  // parse keywords
+
+  nterms = 0;
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (iarg+3 > narg) error->all("Illegal fix adapt command");
+    if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE;
+    else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL;
+    else {
+      which[nterms] = PAIR;
+      int n = strlen(&arg[iarg+1][2]) + 1;
+      pstyle[nterms] = new char[n];
+      strcpy(pstyle[nterms],&arg[iarg+1][2]);
+    }
+
+    if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
+      int n = strlen(&arg[iarg+2][2]) + 1;
+      var1[nterms] = new char[n];
+      strcpy(var1[nterms],&arg[iarg+2][2]);
+    } else error->all("Illegal compute ti command");
+    if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) {
+      int n = strlen(&arg[iarg+3][2]) + 1;
+      var2[nterms] = new char[n];
+      strcpy(var2[nterms],&arg[iarg+3][2]);
+    } else error->all("Illegal compute ti command");
+
+    nterms++;
+    iarg += 3;
+  }
+}
+
+/* --------------------------------------------------------------------- */
+
+ComputeTI::~ComputeTI()
+{
+  for (int m = 0; m < nterms; m++) {
+    delete [] var1[m];
+    delete [] var2[m];
+    delete [] pstyle[m];
+  }
+  delete [] which;
+  delete [] ivar1;
+  delete [] ivar2;
+  delete [] var1;
+  delete [] var2;
+  delete [] pptr;
+  delete [] pstyle;
+}
+
+/* --------------------------------------------------------------------- */
+
+void ComputeTI::init()
+{
+  // setup and error checks
+
+  for (int m = 0; m < nterms; m++) {
+    ivar1[m] = input->variable->find(var1[m]);
+    ivar2[m] = input->variable->find(var2[m]);
+    if (ivar1[m] < 0 || ivar2 < 0)
+      error->all("Variable name for compute ti does not exist");
+    if (!input->variable->equalstyle(ivar1[m]) ||
+	!input->variable->equalstyle(ivar2[m]))
+      error->all("Variable for compute ti is invalid style");
+
+    if (which[m] == PAIR) {
+      pptr[m] = force->pair_match(pstyle[m],1);
+      if (pptr[m] == NULL) error->all("Compute ti pair style does not exist");
+
+    } else if (which[m] == TAIL) {
+      if (force->pair == NULL || force->pair->tail_flag == 0) 
+	error->all("Compute ti tail when pair style does not "
+		   "compute tail corrections");
+
+    } else if (which[m] == KSPACE) {
+      if (force->kspace == NULL) 
+	error->all("Compute ti is incompatible with KSpace style");
+    }
+  }
+}
+
+/* --------------------------------------------------------------------- */
+
+double ComputeTI::compute_scalar()
+{
+  double eng,engall,value1,value2;
+
+  invoked_scalar |= INVOKED_SCALAR;
+  if (update->eflag_global != invoked_scalar)
+    error->all("Energy was not tallied on needed timestep");
+
+  double dUdl = 0.0;
+
+  for (int m = 0; m < nterms; m++) {
+    value1 = input->variable->compute_equal(ivar1[m]);
+    value2 = input->variable->compute_equal(ivar2[m]);
+    if (value1 == 0.0) continue;
+
+    if (which[m] == PAIR) {
+      eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul;
+      MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
+      dUdl += engall/value1 * value2;
+
+    } else if (which[m] == TAIL) {
+      double vol = domain->xprd*domain->yprd*domain->zprd;
+      eng = force->pair->etail / vol;
+      dUdl += eng/value1 * value2;
+
+    } else if (which[m] == KSPACE) {
+      eng = force->kspace->energy;
+      dUdl += eng/value1 * value2;
+    }
+  }
+
+  scalar = dUdl;
+  return scalar;
+}

src/compute_ti.h

@@ -0,0 +1,46 @@
+/* ----------------------------------------------------------------------
+   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 COMPUTE_CLASS
+
+ComputeStyle(ti,ComputeTI)
+
+#else
+
+#ifndef COMPUTE_TI_H
+#define COMPUTE_TI_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+  
+class ComputeTI : public Compute {
+ public: 
+  ComputeTI(class LAMMPS *, int, char **);
+  ~ComputeTI();
+  void init();
+  double compute_scalar();
+
+ private:
+  int nterms;
+  int *which;
+  int *ivar1,*ivar2;
+  char **var1,**var2;
+  class Pair **pptr;
+  char **pstyle;
+};
+
+}
+
+#endif
+#endif

src/pair_born.cpp

@@ -403,3 +403,14 @@ double PairBorn::single(int i, int j, int itype, int jtype,
     d[itype][jtype]*r2inv*r6inv - offset[itype][jtype];
   return factor_lj*phiborn;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *PairBorn::extract(char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"a") == 0) return (void *) a;
+  if (strcmp(str,"c") == 0) return (void *) c;
+  if (strcmp(str,"d") == 0) return (void *) d;
+  return NULL;
+}

src/pair_born.h

@@ -37,6 +37,7 @@ class PairBorn : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
+  void *extract(char *, int &);
 
  private:
   double cut_global;

src/pair_buck.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "stdio.h"
 #include "stdlib.h"
+#include "string.h"
 #include "pair_buck.h"
 #include "atom.h"
 #include "comm.h"
@@ -342,3 +343,13 @@ double PairBuck::single(int i, int j, int itype, int jtype,
     offset[itype][jtype];
   return factor_lj*phibuck;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *PairBuck::extract(char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"a") == 0) return (void *) a;
+  if (strcmp(str,"c") == 0) return (void *) c;
+  return NULL;
+}

src/pair_buck.h

@@ -37,6 +37,7 @@ class PairBuck : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
+  void *extract(char *, int &);
 
  private:
   double cut_global;

src/pair_gauss.cpp

@@ -37,12 +37,18 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp) {}
+PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp)
+{
+  nextra = 1;
+  pvector = new double[1];
+}
 
 /* ---------------------------------------------------------------------- */
 
 PairGauss::~PairGauss()
 {
+  delete [] pvector;
+
   if (allocated) {
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
@@ -66,7 +72,8 @@ void PairGauss::compute(int eflag, int vflag)
   evdwl = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
-  
+  int occ = 0;
+
   double **x = atom->x;
   double **f = atom->f;
   int *type = atom->type;
@@ -98,7 +105,12 @@ void PairGauss::compute(int eflag, int vflag)
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
       jtype = type[j];
-  
+      
+      // define a Gaussian well to be occupied if
+      // the site it interacts with is within the force maximum    
+    
+      if (eflag_global && rsq < 0.5/b[itype][jtype]) occ++;
+
       if (rsq < cutsq[itype][jtype]) {
 	r2inv = 1.0/rsq;
 	r = sqrt(rsq);
@@ -120,11 +132,12 @@ void PairGauss::compute(int eflag, int vflag)
 		    offset[itype][jtype]);
 
 	if (evflag) ev_tally(i,j,nlocal,newton_pair,
-			     evdwl,0.0,fpair,delx,dely,delz); 
+			     evdwl,0.0,fpair,delx,dely,delz);
       }
     }    
   }
-  
+
+  if (eflag_global) pvector[0] = occ;
   if (vflag_fdotr) virial_compute();
 }
 
@@ -314,3 +327,12 @@ double PairGauss::single(int i, int j, int itype, int jtype, double rsq,
   fforce = forcelj*r2inv;
   return philj;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void *PairGauss::extract(char *str, int &dim)
+{
+  dim = 2;
+  if (strcmp(str,"a") == 0) return (void *) a;
+  return NULL;
+}

src/pair_gauss.h

@@ -24,28 +24,30 @@ PairStyle(gauss,PairGauss)
 
 namespace LAMMPS_NS {
 
-  class PairGauss : public Pair {
-  public:
-    PairGauss(class LAMMPS *);
-    ~PairGauss();
-    void compute(int,int);
-    void settings(int, char **);
-    void coeff(int, char **);
-    double init_one(int, int);  
-    void write_restart(FILE *);
-    void read_restart(FILE *);
-    void write_restart_settings(FILE *);
-    void read_restart_settings(FILE *);
-    double single(int, int, int, int, double, double, double, double &);
+class PairGauss : public Pair {
+ public:
+  PairGauss(class LAMMPS *);
+  ~PairGauss();
+  void compute(int,int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);  
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+  double single(int, int, int, int, double, double, double, double &);
+  void *extract(char *, int &);
  
-  private: 
-    double cut_global;
-    double **cut;
-    double **a,**b;
-    double **offset;
+ private: 
+  double cut_global;
+  double **cut;
+  double **a,**b;
+  double **offset;
+  
+  void allocate();
+};
 
-    void allocate();
-  };
 }
 
 #endif