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

src/MANYBODY/Install.csh

@@ -10,6 +10,7 @@ if ($1 == 1) then
   cp pair_eam_fs.cpp ..
   cp pair_sw.cpp ..
   cp pair_tersoff.cpp ..
+  cp pair_tersoff_zbl.cpp ..
 
   cp pair_airebo.h ..
   cp pair_eam.h ..
@@ -17,6 +18,7 @@ if ($1 == 1) then
   cp pair_eam_fs.h ..
   cp pair_sw.h ..
   cp pair_tersoff.h ..
+  cp pair_tersoff_zbl.h ..
 
 else if ($1 == 0) then
 
@@ -29,6 +31,7 @@ else if ($1 == 0) then
   rm ../pair_eam_fs.cpp
   rm ../pair_sw.cpp
   rm ../pair_tersoff.cpp
+  rm ../pair_tersoff_zbl.cpp
 
   rm ../pair_airebo.h
   rm ../pair_eam.h
@@ -36,5 +39,6 @@ else if ($1 == 0) then
   rm ../pair_eam_fs.h
   rm ../pair_sw.h
   rm ../pair_tersoff.h
+  rm ../pair_tersoff_zbl.h
 
 endif

src/MANYBODY/pair_tersoff.cpp

@@ -41,6 +41,7 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
   single_enable = 0;
   one_coeff = 1;
 
+  PI = 4.0*atan(1.0);
   PI2 = 2.0*atan(1.0);
   PI4 = atan(1.0);
 

src/MANYBODY/pair_tersoff.h

@@ -21,14 +21,14 @@ namespace LAMMPS_NS {
 class PairTersoff : public Pair {
  public:
   PairTersoff(class LAMMPS *);
-  ~PairTersoff();
+  virtual ~PairTersoff();
   void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
   double init_one(int, int);
 
- private:
+ protected:
   struct Param {
     double lam1,lam2,lam3;
     double c,d,h;
@@ -39,9 +39,11 @@ class PairTersoff : public Pair {
     double c1,c2,c3,c4;
     int ielement,jelement,kelement;
     int powermint;
+    double Z_i,Z_j;
+    double ZBLcut,ZBLexpscale;
   };
   
-  double PI2,PI4;
+  double PI,PI2,PI4;
   double cutmax;                // max cutoff for all elements
   int nelements;                // # of unique elements
   char **elements;              // names of unique elements
@@ -52,9 +54,9 @@ class PairTersoff : public Pair {
   Param *params;                // parameter set for an I-J-K interaction
 
   void allocate();
-  void read_file(char *);
+  virtual void read_file(char *);
   void setup();
-  void repulsive(Param *, double, double &, int, double &);
+  virtual void repulsive(Param *, double, double &, int, double &);
   double zeta(Param *, double, double, double *, double *);
   void force_zeta(Param *, double, double, double &, double &, int, double &);
   void attractive(Param *, double, double, double, double *, double *,
@@ -62,14 +64,14 @@ class PairTersoff : public Pair {
 
   double ters_fc(double, Param *);
   double ters_fc_d(double, Param *);
-  double ters_fa(double, Param *);
+  virtual double ters_fa(double, Param *);
-  double ters_fa_d(double, Param *);
+  virtual double ters_fa_d(double, Param *);
   double ters_bij(double, Param *);
   double ters_bij_d(double, Param *);
   double ters_gijk(double, Param *);
   double ters_gijk_d(double, Param *);
   void ters_zetaterm_d(double, double *, double, double *, double,
-		       double *, double *, double *, Param *);
+			       double *, double *, double *, Param *);
   void costheta_d(double *, double, double *, double,
 		  double *, double *, double *);
 

src/MANYBODY/pair_tersoff_zbl.cpp

@@ -0,0 +1,285 @@
+/* ----------------------------------------------------------------------
+   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: Aidan Thompson (SNL) - original Tersoff implementation
+                        David Farrell (NWU) - ZBL addition
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_tersoff_zbl.h"
+#include "atom.h"
+#include "update.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "force.h"
+#include "comm.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024
+#define DELTA 4
+
+/* ---------------------------------------------------------------------- */
+
+PairTersoffZBL::PairTersoffZBL(LAMMPS *lmp) : PairTersoff(lmp)
+{
+  // hard-wired constants in metal or real units
+  // a0 = Bohr radius
+  // epsilon0 = permittivity of vacuum = q / energy-distance units
+  // e = unit charge
+  // 1 Kcal/mole = 0.043365121 eV
+
+  if (strcmp(update->unit_style,"metal") == 0) {
+    global_a_0 = 0.529;
+    global_epsilon_0 = 0.00552635;
+    global_e = 1.0;
+  } else if (strcmp(update->unit_style,"real") == 0) {
+    global_a_0 = 0.529;
+    global_epsilon_0 = 0.00552635 * 0.043365121;
+    global_e = 1.0;
+  } else error->all("Pair tersoff/zbl requires metal or real units");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffZBL::read_file(char *file)
+{
+  int params_per_line = 21;
+  char **words = new char*[params_per_line+1];
+
+  if (params) delete [] params;
+  params = NULL;
+  nparams = 0;
+
+  // open file on proc 0
+
+  FILE *fp;
+  if (comm->me == 0) {
+    fp = fopen(file,"r");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open Tersoff potential file %s",file);
+      error->one(str);
+    }
+  }
+
+  // read each line out of file, skipping blank lines or leading '#'
+  // store line of params if all 3 element tags are in element list
+
+  int n,nwords,ielement,jelement,kelement;
+  char line[MAXLINE],*ptr;
+  int eof = 0;
+
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+	eof = 1;
+	fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if (ptr = strchr(line,'#')) *ptr = '\0';
+    nwords = atom->count_words(line);
+    if (nwords == 0) continue;
+
+    // concatenate additional lines until have params_per_line words
+
+    while (nwords < params_per_line) {
+      n = strlen(line);
+      if (comm->me == 0) {
+        ptr = fgets(&line[n],MAXLINE-n,fp);
+        if (ptr == NULL) {
+	  eof = 1;
+	  fclose(fp);
+        } else n = strlen(line) + 1;
+      }
+      MPI_Bcast(&eof,1,MPI_INT,0,world);
+      if (eof) break;
+      MPI_Bcast(&n,1,MPI_INT,0,world);
+      MPI_Bcast(line,n,MPI_CHAR,0,world);
+      if (ptr = strchr(line,'#')) *ptr = '\0';
+      nwords = atom->count_words(line);
+    }
+
+    if (nwords != params_per_line)
+      error->all("Incorrect format in Tersoff potential file");
+
+    // words = ptrs to all words in line
+
+    nwords = 0;
+    words[nwords++] = strtok(line," \t\n\r\f");
+    while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
+
+    // ielement,jelement,kelement = 1st args
+    // if all 3 args are in element list, then parse this line
+    // else skip to next line
+
+    for (ielement = 0; ielement < nelements; ielement++)
+      if (strcmp(words[0],elements[ielement]) == 0) break;
+    if (ielement == nelements) continue;
+    for (jelement = 0; jelement < nelements; jelement++)
+      if (strcmp(words[1],elements[jelement]) == 0) break;
+    if (jelement == nelements) continue;
+    for (kelement = 0; kelement < nelements; kelement++)
+      if (strcmp(words[2],elements[kelement]) == 0) break;
+    if (kelement == nelements) continue;
+
+    // load up parameter settings and error check their values
+
+    if (nparams == maxparam) {
+      maxparam += DELTA;
+      params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+					  "pair:params");
+    }
+
+    params[nparams].ielement = ielement;
+    params[nparams].jelement = jelement;
+    params[nparams].kelement = kelement;
+    params[nparams].powerm = atof(words[3]);
+    params[nparams].gamma = atof(words[4]);
+    params[nparams].lam3 = atof(words[5]);
+    params[nparams].c = atof(words[6]);
+    params[nparams].d = atof(words[7]);
+    params[nparams].h = atof(words[8]);
+    params[nparams].powern = atof(words[9]);
+    params[nparams].beta = atof(words[10]);
+    params[nparams].lam2 = atof(words[11]);
+    params[nparams].bigb = atof(words[12]);
+    params[nparams].bigr = atof(words[13]);
+    params[nparams].bigd = atof(words[14]);
+    params[nparams].lam1 = atof(words[15]);
+    params[nparams].biga = atof(words[16]);
+    params[nparams].Z_i = atof(words[17]);
+    params[nparams].Z_j = atof(words[18]);
+    params[nparams].ZBLcut = atof(words[19]);
+    params[nparams].ZBLexpscale = atof(words[20]);
+
+    // currently only allow m exponent of 1 or 3
+
+    params[nparams].powermint = int(params[nparams].powerm);
+
+    if (
+	params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 || 
+	params[nparams].d < 0.0 || params[nparams].powern < 0.0 || 
+	params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 || 
+	params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
+	params[nparams].bigd < 0.0 ||
+	params[nparams].bigd > params[nparams].bigr ||
+	params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
+	params[nparams].powerm - params[nparams].powermint != 0.0 ||
+        (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
+	params[nparams].gamma < 0.0 ||
+	params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 ||
+	params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0)
+      error->all("Illegal Tersoff parameter");
+
+    nparams++;
+  }
+
+  delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce,
+			       int eflag, double &eng)
+{
+  double r,tmp_fc,tmp_fc_d,tmp_exp;
+
+  // Tersoff repulsive portion
+
+  r = sqrt(rsq);
+  tmp_fc = ters_fc(r,param);
+  tmp_fc_d = ters_fc_d(r,param);
+  tmp_exp = exp(-param->lam1 * r);
+  double fforce_ters = param->biga * tmp_exp * (tmp_fc_d - tmp_fc*param->lam1);
+  double eng_ters = tmp_fc * param->biga * tmp_exp;
+	
+  // ZBL repulsive portion
+
+  double esq = pow(global_e,2.0);
+  double a_ij = (0.8854*global_a_0) / 
+    (pow(param->Z_i,0.23) + pow(param->Z_j,0.23));
+  double premult = (param->Z_i * param->Z_j * esq)/(4.0*PI*global_epsilon_0);
+  double r_ov_a = r/a_ij;
+  double phi = 0.1818*exp(-3.2*r_ov_a) + 0.5099*exp(-0.9423*r_ov_a) + 
+    0.2802*exp(-0.4029*r_ov_a) + 0.02817*exp(-0.2016*r_ov_a);
+  double dphi = (1.0/a_ij) * (-3.2*0.1818*exp(-3.2*r_ov_a) - 
+			      0.9423*0.5099*exp(-0.9423*r_ov_a) - 
+			      0.4029*0.2802*exp(-0.4029*r_ov_a) - 
+			      0.2016*0.02817*exp(-0.2016*r_ov_a));
+  double fforce_ZBL = premult*-pow(r,-2.0)* phi + premult*pow(r,-1.0)*dphi;
+  double eng_ZBL = premult*(1.0/r)*phi;
+  
+  // combine two parts with smoothing by Fermi-like function
+
+  fforce = -(-F_fermi_d(r,param) * eng_ZBL + 
+	     (1.0 - F_fermi(r,param))*fforce_ZBL + 
+	     F_fermi_d(r,param)*eng_ters + F_fermi(r,param)*fforce_ters) / r;
+  
+  if (eflag)
+    eng = (1.0 - F_fermi(r,param))*eng_ZBL + F_fermi(r,param)*eng_ters;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairTersoffZBL::ters_fa(double r, Param *param)
+{
+  if (r > param->bigr + param->bigd) return 0.0;
+  return -param->bigb * exp(-param->lam2 * r) * ters_fc(r,param) * 
+    F_fermi(r,param);
+}   
+
+/* ---------------------------------------------------------------------- */
+
+double PairTersoffZBL::ters_fa_d(double r, Param *param)
+{
+  if (r > param->bigr + param->bigd) return 0.0;
+  return param->bigb * exp(-param->lam2 * r) *
+    (param->lam2 * ters_fc(r,param) * F_fermi(r,param) - 
+     ters_fc_d(r,param) * F_fermi(r,param) - ters_fc(r,param) * 
+     F_fermi_d(r,param));
+}
+
+/* ----------------------------------------------------------------------
+   Fermi-like smoothing function
+------------------------------------------------------------------------- */
+
+double PairTersoffZBL::F_fermi(double r, Param *param)
+{
+  return 1.0 / (1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)));
+}
+
+/* ----------------------------------------------------------------------
+   Fermi-like smoothing function derivative with respect to r
+------------------------------------------------------------------------- */
+
+double PairTersoffZBL::F_fermi_d(double r, Param *param)
+{
+  return param->ZBLexpscale*exp(-param->ZBLexpscale*(r-param->ZBLcut)) / 
+    pow(1.0 + exp(-param->ZBLexpscale*(r-param->ZBLcut)),2.0);
+}

src/MANYBODY/pair_tersoff_zbl.h

@@ -0,0 +1,43 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef PAIR_TERSOFF_ZBL_H
+#define PAIR_TERSOFF_ZBL_H
+
+#include "pair_tersoff.h"
+
+namespace LAMMPS_NS {
+
+class PairTersoffZBL : public PairTersoff {
+ public:
+  PairTersoffZBL(class LAMMPS *);
+  ~PairTersoffZBL() {}
+
+ private:
+  double global_a_0;		// Bohr radius for Coulomb repulsion
+  double global_epsilon_0;	// permittivity of vacuum for Coulomb repulsion
+  double global_e;		// proton charge (negative of electron charge)
+
+  void read_file(char *);
+  void repulsive(Param *, double, double &, int, double &);
+
+  double ters_fa(double, Param *);
+  double ters_fa_d(double, Param *);
+	
+  double F_fermi(double, Param *);
+  double F_fermi_d(double, Param *);
+};
+
+}
+
+#endif

src/MANYBODY/style_manybody.h

@@ -18,6 +18,7 @@
 #include "pair_eam_fs.h"
 #include "pair_sw.h"
 #include "pair_tersoff.h"
+#include "pair_tersoff_zbl.h"
 #endif
 
 #ifdef PairClass
@@ -27,4 +28,5 @@ PairStyle(eam/alloy,PairEAMAlloy)
 PairStyle(eam/fs,PairEAMFS)
 PairStyle(sw,PairSW)
 PairStyle(tersoff,PairTersoff)
+PairStyle(tersoff/zbl,PairTersoffZBL)
 #endif

src/style_manybody.h

@@ -18,6 +18,7 @@
 #include "pair_eam_fs.h"
 #include "pair_sw.h"
 #include "pair_tersoff.h"
+#include "pair_tersoff_zbl.h"
 #endif
 
 #ifdef PairClass
@@ -27,4 +28,5 @@ PairStyle(eam/alloy,PairEAMAlloy)
 PairStyle(eam/fs,PairEAMFS)
 PairStyle(sw,PairSW)
 PairStyle(tersoff,PairTersoff)
+PairStyle(tersoff/zbl,PairTersoffZBL)
 #endif