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

src/fix_adapt.cpp

@@ -11,6 +11,7 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "math.h"
 #include "string.h"
 #include "stdlib.h"
 #include "fix_adapt.h"
@@ -148,14 +149,17 @@ void FixAdapt::init()
 	error->all("Fix adapt pair style does not exist");
       pairindex[m] = 
 	pairptr[m]->pre_adapt(param[m],ilo[m],ihi[m],jlo[m],jhi[m]);
-      if (pairindex[m] < 0)
+      if (pairindex[m] == -1)
-	error->all("Fix adapt pair style is not compatible");
+	error->all("Fix adapt pair parameter is not recognized");
+      if (pairindex[m] == -2)
+	error->all("Fix adapt pair types are not valid");
+
     } else if (which[m] == ATOM) {
       if (strcmp(param[m],"diameter") == 0) {
 	awhich[m] = DIAMETER;
 	if (!atom->radius_flag)
 	  error->all("Fix adapt requires atom attribute radius");
-      } else error->all("Invalid fix adapt atom attribute");
+      } else error->all("Fix adapt atom attribute is not recognized");
     }
 
     ivar[m] = input->variable->find(var[m]);
@@ -182,14 +186,26 @@ void FixAdapt::pre_force(int vflag)
       pairptr[m]->adapt(pairindex[m],ilo[m],ihi[m],jlo[m],jhi[m],value);
 
     else if (which[m] == ATOM) {
-      double *radius = atom->radius;
+
-      int *mask = atom->mask;
+      // set radius from diameter
-      int nlocal = atom->nlocal;
+      // set rmass if both rmass and density are defined
 
       if (awhich[m] == DIAMETER) {
+	int mflag = 0;
+	if (atom->rmass_flag && atom->density_flag) mflag = 1;
+	double PI = 4.0*atan(1.0);
+
+	double *radius = atom->radius;
+	double *rmass = atom->rmass;
+	double *density = atom->density;
+	int *mask = atom->mask;
+	int nlocal = atom->nlocal;
+
 	for (int i = 0; i < nlocal; i++)
-	  if (mask[i] & groupbit)
+	  if (mask[i] & groupbit) {
 	    radius[i] = 0.5*value;
+	    rmass[i] = 4.0*PI/3.0 * radius[i]*radius[i]*radius[i] * density[i];
+	  }
       }
     }
   }

src/pair.cpp

@@ -23,7 +23,6 @@
 #include "stdlib.h"
 #include "string.h"
 #include "pair.h"
-#include "pair_soft.h"
 #include "atom.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -895,13 +894,6 @@ void Pair::write_file(int narg, char **arg)
 
   force->init();
 
-  // if pair style = soft, set prefactor to final value
-
-  Pair *spair = force->pair_match("soft",1);
-  if (spair)
-    ((PairSoft *) spair)->prefactor[itype][jtype] =
-      ((PairSoft *) spair)->prestop[itype][jtype];
-
   // if pair style = any of EAM, swap in dummy fp vector
 
   double eamfp[2];

src/pair_hybrid.h

@@ -32,7 +32,7 @@ class PairHybrid : public Pair {
   char **keywords;              // style name of each Pair style
 
   PairHybrid(class LAMMPS *);
-  ~PairHybrid();
+  virtual ~PairHybrid();
   void compute(int, int);
   void settings(int, char **);
   virtual void coeff(int, char **);

src/pair_hybrid_overlay.h

@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class PairHybridOverlay : public PairHybrid {
  public:
   PairHybridOverlay(class LAMMPS *);
+  ~PairHybridOverlay() {}
   void coeff(int, char **);
 
  private:

src/pair_soft.cpp

@@ -14,6 +14,7 @@
 #include "math.h"
 #include "stdio.h"
 #include "stdlib.h"
+#include "string.h"
 #include "pair_soft.h"
 #include "atom.h"
 #include "comm.h"
@@ -43,8 +44,6 @@ PairSoft::~PairSoft()
     memory->destroy_2d_int_array(setflag);
     memory->destroy_2d_double_array(cutsq);
 
-    memory->destroy_2d_double_array(prestart);
-    memory->destroy_2d_double_array(prestop);
     memory->destroy_2d_double_array(prefactor);
     memory->destroy_2d_double_array(cut);
   }
@@ -63,21 +62,6 @@ void PairSoft::compute(int eflag, int vflag)
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
-  // set current prefactor
-  // for minimization, set to prestop
-  // for dynamics, ramp from prestart to prestop
-  // for 0-step dynamics, set to prestart
-
-  double delta = update->ntimestep - update->beginstep;
-  if (update->whichflag == 2) delta = 1.0;
-  else if (update->nsteps) delta /= update->endstep - update->beginstep;
-  else delta = 0.0;
-  int ntypes = atom->ntypes;
-  for (i = 1; i <= ntypes; i++)
-    for (j = 1; j <= ntypes; j++)
-      prefactor[i][j] = prestart[i][j] + 
-	delta * (prestop[i][j] - prestart[i][j]);
-
   double **x = atom->x;
   double **f = atom->f;
   int *type = atom->type;
@@ -161,20 +145,8 @@ void PairSoft::allocate()
 
   cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
 
-  prestart = memory->create_2d_double_array(n+1,n+1,"pair:prestart");
-  prestop = memory->create_2d_double_array(n+1,n+1,"pair:prestop");
   prefactor = memory->create_2d_double_array(n+1,n+1,"pair:prefactor");
   cut = memory->create_2d_double_array(n+1,n+1,"pair:cut");
-
-  // init prestart and prestop to 0.0
-  // since pair_hybrid can use all types even if pair_soft sub-class
-  //   never sets them
-
-  for (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++) {
-      prestart[i][j] = 0.0;
-      prestop[i][j] = 0.0;
-    }
 }
 
 /* ----------------------------------------------------------------------
@@ -203,24 +175,22 @@ void PairSoft::settings(int narg, char **arg)
 
 void PairSoft::coeff(int narg, char **arg)
 {
-  if (narg < 4 || narg > 5) error->all("Incorrect args for pair coefficients");
+  if (narg < 3 || narg > 4) 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 prestart_one = force->numeric(arg[2]);
+  double prefactor_one = force->numeric(arg[2]);
-  double prestop_one = force->numeric(arg[3]);
 
   double cut_one = cut_global;
-  if (narg == 5) cut_one = force->numeric(arg[4]);
+  if (narg == 4) cut_one = force->numeric(arg[3]);
 
   int count = 0;
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
-      prestart[i][j] = prestart_one;
+      prefactor[i][j] = prefactor_one;
-      prestop[i][j] = prestop_one;
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
       count++;
@@ -239,13 +209,11 @@ double PairSoft::init_one(int i, int j)
   // always mix prefactors geometrically
 
   if (setflag[i][j] == 0) {
-    prestart[i][j] = sqrt(prestart[i][i]*prestart[j][j]);
+    prefactor[i][j] = sqrt(prefactor[i][i]*prefactor[j][j]);
-    prestop[i][j] = sqrt(prestop[i][i]*prestop[j][j]);
     cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
   }
 
-  prestart[j][i] = prestart[i][j];
+  prefactor[j][i] = prefactor[i][j];
-  prestop[j][i] = prestop[i][j];
   cut[j][i] = cut[i][j];
 
   return cut[i][j];
@@ -264,8 +232,7 @@ void PairSoft::write_restart(FILE *fp)
     for (j = i; j <= atom->ntypes; j++) {
       fwrite(&setflag[i][j],sizeof(int),1,fp);
       if (setflag[i][j]) {
-	fwrite(&prestart[i][j],sizeof(double),1,fp);
+	fwrite(&prefactor[i][j],sizeof(double),1,fp);
-	fwrite(&prestop[i][j],sizeof(double),1,fp);
 	fwrite(&cut[i][j],sizeof(double),1,fp);
       }
     }
@@ -289,12 +256,10 @@ void PairSoft::read_restart(FILE *fp)
       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
       if (setflag[i][j]) {
 	if (me == 0) {
-	  fread(&prestart[i][j],sizeof(double),1,fp);
+	  fread(&prefactor[i][j],sizeof(double),1,fp);
-	  fread(&prestop[i][j],sizeof(double),1,fp);
 	  fread(&cut[i][j],sizeof(double),1,fp);
 	}
-	MPI_Bcast(&prestart[i][j],1,MPI_DOUBLE,0,world);
+	MPI_Bcast(&prefactor[i][j],1,MPI_DOUBLE,0,world);
-	MPI_Bcast(&prestop[i][j],1,MPI_DOUBLE,0,world);
 	MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
       }
     }
@@ -324,6 +289,38 @@ void PairSoft::read_restart_settings(FILE *fp)
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }
 
+/* ----------------------------------------------------------------------
+   check if name is recognized, return integer index for that name
+   if name not recognized, return -1
+   if type pair setting, return -2 if no type pairs are set
+------------------------------------------------------------------------- */
+
+int PairSoft::pre_adapt(char *name, int ilo, int ihi, int jlo, int jhi)
+{
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      count++;
+  if (count == 0) return -2;
+
+  if (strcmp(name,"a") == 0) return 0;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   adapt parameter indexed by which
+   change all pair variables affected by the reset parameter
+   if type pair setting, set I-J and J-I coeffs
+------------------------------------------------------------------------- */
+
+void PairSoft::adapt(int which, int ilo, int ihi, int jlo, int jhi,
+		     double value)
+{
+  for (int i = ilo; i <= ihi; i++)
+    for (int j = MAX(jlo,i); j <= jhi; j++)
+      prefactor[i][j] = prefactor[j][i] = value;
+}
+
 /* ---------------------------------------------------------------------- */
 
 double PairSoft::single(int i, int j, int itype, int jtype, double rsq,

src/pair_soft.h

@@ -38,12 +38,13 @@ class PairSoft : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  int pre_adapt(char *, int, int, int, int);
+  void adapt(int, int, int, int, int, double);
   double single(int, int, int, int, double, double, double, double &);
 
  private:
   double PI;
   double cut_global;
-  double **prestart,**prestop;
   double **prefactor;
   double **cut;