From f62f8fc6e6a0ba26619eb1599ba76f90860a21e4 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 1 Oct 2008 15:02:36 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2121 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/MANYBODY/Install.csh | 4 + src/MANYBODY/pair_tersoff.cpp | 1 + src/MANYBODY/pair_tersoff.h | 18 +- src/MANYBODY/pair_tersoff_zbl.cpp | 285 ++++++++++++++++++++++++++++++ src/MANYBODY/pair_tersoff_zbl.h | 43 +++++ src/MANYBODY/style_manybody.h | 2 + src/style_manybody.h | 2 + 7 files changed, 347 insertions(+), 8 deletions(-) create mode 100644 src/MANYBODY/pair_tersoff_zbl.cpp create mode 100755 src/MANYBODY/pair_tersoff_zbl.h diff --git a/src/MANYBODY/Install.csh b/src/MANYBODY/Install.csh index 5cfed7a791..3857346082 100644 --- a/src/MANYBODY/Install.csh +++ b/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 diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 6e2d75ebf0..3a70e6cd63 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/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); diff --git a/src/MANYBODY/pair_tersoff.h b/src/MANYBODY/pair_tersoff.h index bb291bce8a..7769ddcecc 100755 --- a/src/MANYBODY/pair_tersoff.h +++ b/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 *); - double ters_fa_d(double, Param *); + virtual double ters_fa(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 *); diff --git a/src/MANYBODY/pair_tersoff_zbl.cpp b/src/MANYBODY/pair_tersoff_zbl.cpp new file mode 100644 index 0000000000..fcc707815c --- /dev/null +++ b/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); +} diff --git a/src/MANYBODY/pair_tersoff_zbl.h b/src/MANYBODY/pair_tersoff_zbl.h new file mode 100755 index 0000000000..dd3f110886 --- /dev/null +++ b/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 diff --git a/src/MANYBODY/style_manybody.h b/src/MANYBODY/style_manybody.h index c564993e3e..18cd907008 100644 --- a/src/MANYBODY/style_manybody.h +++ b/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 diff --git a/src/style_manybody.h b/src/style_manybody.h index c564993e3e..18cd907008 100644 --- a/src/style_manybody.h +++ b/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