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

src/MANYBODY/fix_qeq_comb.cpp

@@ -12,7 +12,8 @@
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
-   Contributing author: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
+
+   Contributing authors: Ray Shan (Sandia, tnshan@sandia.gov)
 ------------------------------------------------------------------------- */
 
 #include "lmptype.h"
@@ -21,6 +22,7 @@
 #include "stdlib.h"
 #include "string.h"
 #include "pair_comb.h"
+#include "pair_comb3.h"
 #include "fix_qeq_comb.h"
 #include "neighbor.h"
 #include "neigh_list.h"
@@ -88,6 +90,7 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
   for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
 
   comb = NULL;
+  comb3 = NULL;
 
   comm_forward = 1;
 }
@@ -121,8 +124,9 @@ void FixQEQComb::init()
     error->all(FLERR,"Fix qeq/comb requires atom attribute q");
 
   comb = (PairComb *) force->pair_match("comb",1);
-  if (comb == NULL)
+  comb3 = (PairComb3 *) force->pair_match("comb3",1);
-    error->all(FLERR,"Must use pair_style comb with fix qeq/comb");
+  if (comb == NULL && comb3 == NULL)
+    error->all(FLERR,"Must use pair_style comb or comb3 with fix qeq/comb");
 
   if (strstr(update->integrate_style,"respa"))
     nlevels_respa = ((Respa *) update->integrate)->nlevels;
@@ -188,8 +192,8 @@ void FixQEQComb::post_force(int vflag)
   // more loops for first-time charge equilibrium
 
   iloop = 0;
-  if (firstflag) loopmax = 500;
+  if (firstflag) loopmax = 200;
-  else loopmax = 200;
+  else loopmax = 100;
 
   // charge-equilibration loop
 
@@ -211,9 +215,14 @@ void FixQEQComb::post_force(int vflag)
   int *tag = atom->tag;
   int nlocal = atom->nlocal;
 
+ if (comb) {
     inum = comb->list->inum;
     ilist = comb->list->ilist;
-
+  }
+  if (comb3) {
+    inum = comb3->list->inum;
+    ilist = comb3->list->ilist;
+  }
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
     q1[i] = q2[i] = qf[i] = 0.0;
@@ -227,9 +236,11 @@ void FixQEQComb::post_force(int vflag)
         q[i]  += q1[i];
       }
     }
-    comm->forward_comm_fix(this);
 
+    comm->forward_comm_fix(this);
     if(comb) enegtot = comb->yasu_char(qf,igroup);
+    if(comb3) enegtot = comb3->combqeq(qf,igroup);
+
     enegtot /= ngroup;
     enegchk = enegmax = 0.0;
 

src/MANYBODY/fix_qeq_comb.h

@@ -48,6 +48,7 @@ class FixQEQComb : public Fix {
   FILE *fp;
 
   class PairComb *comb;
+  class PairComb3 *comb3;
   int nmax;
   double *qf,*q1,*q2;
 };

src/MANYBODY/pair_comb3.h

@@ -0,0 +1,333 @@
+/* ----------------------------------------------------------------------
+   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 PAIR_CLASS
+
+PairStyle(comb3,PairComb3)
+
+#else
+
+#ifndef LMP_PAIR_COMB3_H
+#define LMP_PAIR_COMB3_H
+
+#include "pair.h"
+#include "my_page.h"
+
+namespace LAMMPS_NS {
+
+class PairComb3 : public Pair {
+ public:
+  PairComb3(class LAMMPS *);
+  virtual ~PairComb3();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  double init_one(int, int);
+  double memory_usage();
+  virtual double combqeq(double *, int &);
+
+ // general potential parameters
+ protected:
+  struct Param {
+    int ielement,jelement,kelement,powermint;
+    int ielementgp,jelementgp,kelementgp;  	//element group
+    int ang_flag,pcn_flag,rad_flag,tor_flag;  	//angle, coordination,radical, torsion flag
+    double lami,lambda,alfi,alpha1,alpha2,alpha3,beta;
+    double pcos6,pcos5,pcos4,pcos3,pcos2,pcos1,pcos0;
+    double gamma,powerm,powern,bigA,bigB1,bigB2,bigB3;
+    double bigd,bigr,cut,cutsq,c1,c2,c3,c4;
+    double plp1,plp2,plp3,plp4,plp5,plp6;
+    double pbb1,pbb2,ptork1,ptork2;
+    double addrep, vdwflag;
+    double QU,QL,DU,DL,Qo,dQ,aB,bB,nD,bD,qmin,qmax;
+    double chi,dj,dk,dl,dm,esm,cmn1,cmn2,pcmn1,pcmn2;
+    double coulcut, lcut, lcutsq;
+    double veps, vsig, pcna, pcnb, pcnc, pcnd, polz, curl, pcross;
+    double paaa, pbbb; 
+    double curlcut1, curlcut2, curl0;
+  };
+
+  // general setups
+  int nelements;                	// # of unique elements
+  int ***elem2param;            	// mapping from element triplets to parameters
+  int *map;                     	// mapping from atom types to elements
+  int nparams;                  	// # of stored parameter sets
+  int maxparam;                 	// max # of parameter sets
+  double PI,PI2,PI4,PIsq;		// PIs
+  double cutmin;                	// min cutoff for all elements
+  double cutmax;                	// max cutoff for all elements
+  double precision;			// tolerance for QEq convergence
+  char **elements;              	// names of unique elements
+  Param *params;                	// parameter set for an I-J-K interaction
+  int debug_eng1, debug_eng2, debug_fq;	// logic controlling debugging outputs
+  int pack_flag;
+
+  // Short range neighbor list
+  void Short_neigh();
+  int pgsize, oneatom;
+  int *sht_num, **sht_first;	
+  MyPage<int> *ipage;
+
+  // loop up tables and flags
+  int nmax, **intype;
+  int  pol_flag, polar;
+  double *qf, **bbij, *charge, *NCo;
+  double *esm, **fafb, **dfafb, **ddfafb, **phin, **dphin, **erpaw;
+  double **vvdw, **vdvdw;
+  double **afb, **dafb;
+  double **dpl, bytes;
+  double *xcctmp, *xchtmp, *xcotmp;
+
+  // additional carbon parameters
+  int cflag;
+  int nsplpcn,nsplrad,nspltor;
+  int maxx,maxy,maxz,maxxc,maxyc,maxconj;
+  int maxxcn[4];            
+  double vmaxxcn[4],dvmaxxcn[4];
+  int ntab;
+  double iin2[16][2],iin3[64][3];
+  double brad[4], btor[4], bbtor, ptorr;
+  double fi_tor[3], fj_tor[3], fk_tor[3], fl_tor[3];
+  double radtmp, fi_rad[3], fj_rad[3], fk_rad[3];
+
+  double ccutoff[6],ch_a[7];
+
+  //COMB3-v18 arrays for CHO
+	// We wanna dynamic arrays
+	// C angle arrays, size = ntab+1
+  double pang[20001];
+  double dpang[20001];
+  double ddpang[20001];
+
+  //coordination spline arrays
+  double pcn_grid[4][5][5][5]; 
+  double pcn_gridx[4][5][5][5];
+  double pcn_gridy[4][5][5][5];
+  double pcn_gridz[4][5][5][5];
+  double pcn_cubs[4][4][4][4][64]; 
+
+  //coordination spline arrays
+  double rad_grid[3][5][5][11];
+  double rad_gridx[3][5][5][11];
+  double rad_gridy[3][5][5][11];
+  double rad_gridz[3][5][5][11];
+  double rad_spl[3][4][4][10][64];
+
+  //torsion spline arrays
+  double tor_grid[1][5][5][11]; 
+  double tor_gridx[1][5][5][11];
+  double tor_gridy[1][5][5][11];
+  double tor_gridz[1][5][5][11];
+  double tor_spl[1][4][4][10][64];
+
+  // initialization functions
+  void allocate();
+  void read_lib();
+  void setup();
+  virtual void read_file(char *);
+
+  // cutoff functions
+  double comb_fc(double, Param *);
+  double comb_fc_d(double, Param *);
+  double comb_fc_curl(double, Param *);
+  double comb_fc_curl_d(double, Param *);
+  double comb_fccc(double);
+  double comb_fccc_d(double);
+  double comb_fcch(double);
+  double comb_fcch_d(double);
+  double comb_fccch(double);
+  double comb_fccch_d(double);
+  double comb_fcsw(double);
+
+  // short range terms
+  void attractive(Param *, Param *, Param *, double, double, double, double,
+	double, double, double, double *, double *, double *, 
+	double *, double *, int, double);
+  virtual void comb_fa(double, Param *, Param *, double, double, 
+	  double &, double &);
+  virtual void repulsive(Param *, Param *,double, double &, int, 
+	 double &, double, double);
+
+  // bond order terms
+  double comb_bij(double, Param *, double, int, double);
+  double comb_gijk(double, Param *, double);
+  void comb_gijk_d(double, Param *, double, double &, double &);
+  double zeta(Param *, Param *, double, double, double *, double *, int, double);
+  void comb_bij_d(double, Param *, double, int, double &,
+	  double &, double &, double &, double &, double &, double);
+  void coord(Param *, double, int, double &, double &, 
+	  double &, double &, double &, double);
+  void comb_zetaterm_d(double, double, double, double, double, 
+	double *, double, double *, double, double *, double *, 
+	double *, Param *, Param *, Param *, double);
+  void costheta_d(double *, double, double *, double,
+	  double *, double *, double *);
+  void force_zeta(Param *, Param *, double, double, double, double &, 
+	double &, double &, double &, double &, double &, double &, 
+	double &, double &, double &, double &, double &, double &, 
+	int, double &, double,double, int, int, int, 
+	double , double , double);
+  void cntri_int(int, double, double, double, int, int, int,
+        double &, double &, double &, double &, Param *);
+
+  // Legendre polynomials
+  void selflp(Param *, Param *, double, double &, double &);
+  double elp(Param *, Param *, double, double, double *, double * ,double &);
+  void flp(Param *, Param *, double, double, double *, double *, double *, 
+	  double *, double *);
+
+  // long range q-dependent terms
+  double self(Param *, double);
+  void tables();
+  void potal_calc(double &, double &, double &);
+  void tri_point(double, int &, int &, int &, double &, double &, 
+	 double &);
+  void vdwaals(int,int,int,int,double,double,double,double,
+	  double &, double &);
+  void direct(Param *, Param *, int,int,int,double,double,
+	double,double,double,double, double,double,double &,double &,
+         int, int);
+  void field(Param *, Param *,double,double,double,double &,double &);
+  int heaviside(double);
+  double switching(double);
+  double switching_d(double);
+  double chicut1, chicut2;
+
+  // radical terms
+  double rad_init(double, Param *, int, double &, double);
+  void rad_calc(double, Param *, Param *, double, double, int, 
+		  int, double, double);
+  void rad_int(int , double, double, double, int, int, int,
+	double &, double &, double &, double &);
+  void rad_forceik(Param *,  double, double *, double, double);
+  void rad_force(Param *,  double, double *,  double);
+
+  // torsion terms
+  double bbtor1(int, Param *, Param *, double, double, double, 
+        double *, double *, double *, double);                     //modified by TAO
+  void tor_calc(double, Param *, Param *, double, double, int, 
+		  int, double, double);
+  void tor_int(int , double, double, double, int, int, int,
+        double &, double &, double &, double &);
+  void tor_force(int, Param *, Param *, double, double, double, 
+        double, double *, double *, double *);                //modified by TAO
+
+  // charge force terms
+  double qfo_self(Param *, double);
+  void qfo_short(Param *, Param *, double, double, double, 
+	double &, double &, int, int, int);
+  void qfo_direct(Param *, Param *, int, int, int, double, 
+	double, double, double, double, double &, double &, 
+	double, double, int, int);
+  void qfo_field(Param *, Param *,double,double ,double ,double &, double &);
+  void qfo_dipole(double, int, int, int, int, double, double *, double, 
+	double, double, double &, double &, int, int);
+  void qsolve(double *);
+  
+  // dipole - polarization terms
+  double dipole_self(Param *, int);
+  void dipole_init(Param *, Param *, double, double *, double,
+	int, int, int, double, double, double, double, double, int , int);
+  void dipole_calc(Param *, Param *, double, double, double, double, double,
+	int, int, int, double, double, double, double, double, int , int,
+	double &, double &, double *);
+
+  // communication functions
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  int pack_comm(int , int *, double *, int, int *);
+  void unpack_comm(int , int , double *);
+
+  // vector functions, inline for efficiency
+  inline double vec3_dot(double *x, double *y) {
+    return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
+  }
+  inline void vec3_add(double *x, double *y, double *z) {
+    z[0] = x[0]+y[0];  z[1] = x[1]+y[1];  z[2] = x[2]+y[2];
+  }
+  inline void vec3_scale(double k, double *x, double *y) {
+    y[0] = k*x[0];  y[1] = k*x[1];  y[2] = k*x[2];
+  }
+  inline void vec3_scaleadd(double k, double *x, double *y, double *z) {
+    z[0] = k*x[0]+y[0];  z[1] = k*x[1]+y[1];  z[2] = k*x[2]+y[2];
+  }
+ };
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair style COMB requires atom IDs
+
+This is a requirement to use the AIREBO potential.
+
+E: Pair style COMB requires newton pair on
+
+See the newton command.  This is a restriction to use the COMB
+potential.
+
+E: Pair style COMB requires atom attribute q
+
+Self-explanatory.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Cannot open COMB potential file %s
+
+The specified COMB potential file cannot be opened.  Check that the
+path and name are correct.
+
+E: Incorrect format in COMB potential file
+
+Incorrect number of words per line in the potential file.
+
+E: Illegal COMB parameter
+
+One or more of the coefficients defined in the potential file is
+invalid.
+
+E: Potential file has duplicate entry
+
+The potential file for a SW or Tersoff potential has more than
+one entry for the same 3 ordered elements.
+
+E: Potential file is missing an entry
+
+The potential file for a SW or Tersoff potential does not have a
+needed entry.
+
+W: Pair COMB charge %.10f with force %.10f hit min barrier
+
+Something is possibly wrong with your model.
+
+W: Pair COMB charge %.10f with force %.10f hit max barrier
+
+Something is possibly wrong with your model.
+
+*/