From 0dacc993ae27c509f8389eee5774c2ab5b7febb5 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 1 Dec 2011 14:50:34 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7254 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/pair_coul_diel.cpp | 348 +++++++++ src/USER-MISC/pair_coul_diel.h | 57 ++ src/USER-MISC/pair_gauss_cut.cpp | 371 ++++++++++ src/USER-MISC/pair_gauss_cut.h | 60 ++ src/USER-MISC/pair_tersoff_table.cpp | 1007 ++++++++++++++++++++++++++ src/USER-MISC/pair_tersoff_table.h | 96 +++ 6 files changed, 1939 insertions(+) create mode 100644 src/USER-MISC/pair_coul_diel.cpp create mode 100644 src/USER-MISC/pair_coul_diel.h create mode 100644 src/USER-MISC/pair_gauss_cut.cpp create mode 100644 src/USER-MISC/pair_gauss_cut.h create mode 100644 src/USER-MISC/pair_tersoff_table.cpp create mode 100644 src/USER-MISC/pair_tersoff_table.h diff --git a/src/USER-MISC/pair_coul_diel.cpp b/src/USER-MISC/pair_coul_diel.cpp new file mode 100644 index 0000000000..7f63b1a201 --- /dev/null +++ b/src/USER-MISC/pair_coul_diel.cpp @@ -0,0 +1,348 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributiong authors: Arben Jusufi, Axel Kohlmeyer (Temple U.) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_diel.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairCoulDiel::PairCoulDiel(LAMMPS *lmp) : Pair(lmp) {} + +/* ---------------------------------------------------------------------- */ + +PairCoulDiel::~PairCoulDiel() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(sigmae); + memory->destroy(rme); + memory->destroy(offset); + memory->destroy(cutsq); + memory->destroy(cut); + allocated = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulDiel::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair; + double rsq,r,rarg,th,depsdr,epsr,forcecoul,factor_coul; + int *ilist,*jlist,*numneigh,**firstneigh; + + ecoul = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + int newton_pair = force->newton_pair; + double qqrd2e = force->qqrd2e; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + qtmp = q[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + rarg = (r-rme[itype][jtype])/sigmae[itype][jtype]; + th=tanh(rarg); + epsr=a_eps+b_eps*th; + depsdr=b_eps * (1.0 - th*th) / sigmae[itype][jtype]; + + forcecoul = qqrd2e*qtmp*q[j]*((eps_s*(epsr+r*depsdr)/epsr/epsr) -1.)/rsq; + fpair = factor_coul*forcecoul/r; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + ecoul = (qqrd2e*qtmp*q[j]*((eps_s/epsr) -1.)/r) - offset[itype][jtype]; + ecoul *= factor_coul; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair,0.0, + ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairCoulDiel::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(sigmae,n+1,n+1,"pair:sigmae"); + memory->create(rme,n+1,n+1,"pair:rme"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCoulDiel::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + + cut_global = force->numeric(arg[0]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulDiel::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) error->all(FLERR,"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); + + eps_s = force->numeric(arg[2]); + double rme_one =force->numeric(arg[3]); + double sigmae_one = force->numeric(arg[4]); + + double cut_one = cut_global; + if (narg == 6) cut_one = force->numeric(arg[5]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + sigmae[i][j] = sigmae_one; + rme[i][j] = rme_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + a_eps = 0.5*(5.2+eps_s); + b_eps = 0.5*(eps_s-5.2); + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairCoulDiel::init_style() +{ + if (!atom->q_flag) + error->all(FLERR,"Pair style coul/diel requires atom attribute q"); + + int irequest = neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairCoulDiel::init_one(int i, int j) +{ + if (setflag[i][j] == 0) { + error->all(FLERR,"for pair style coul/diel, parameters need to be set explicitly for all pairs."); + } + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + if (offset_flag) { + double rarg = (cut[i][j]-rme[i][j])/sigmae[i][j]; + double epsr=a_eps+b_eps*tanh(rarg); + offset[i][j] = qqrd2e*q[i]*q[j]*((eps_s/epsr) -1.)/cut[i][j]; + } else offset[i][j] = 0.0; + + + sigmae[j][i] = sigmae[i][j]; + rme[j][i] = rme[i][j]; + offset[j][i] = offset[i][j]; + cut[j][i] = cut[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulDiel::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&rme[i][j],sizeof(double),1,fp); + fwrite(&sigmae[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairCoulDiel::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (setflag[i][j]) { + if (me == 0) { + fread(&rme[i][j],sizeof(double),1,fp); + fread(&sigmae[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&rme[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigmae[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulDiel::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairCoulDiel::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairCoulDiel::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, + double &fforce) +{ + double r, rarg,forcedielec,phidielec; + double th,epsr,depsdr; + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r=sqrt(rsq); + rarg = (r-rme[itype][jtype])/sigmae[itype][jtype]; + th = tanh(rarg); + epsr=a_eps+b_eps*th; + depsdr=b_eps*(1.-th*th)/sigmae[itype][jtype]; + + forcedielec = qqrd2e*q[i]*q[j]*((eps_s*(epsr+r*depsdr)/epsr/epsr) -1.)/rsq; + fforce = factor_coul*forcedielec/r; + + phidielec = (qqrd2e*q[i]*q[j]*((eps_s/epsr) -1.)/r)- offset[itype][jtype]; + return factor_coul*phidielec; +} diff --git a/src/USER-MISC/pair_coul_diel.h b/src/USER-MISC/pair_coul_diel.h new file mode 100644 index 0000000000..14d7b8b143 --- /dev/null +++ b/src/USER-MISC/pair_coul_diel.h @@ -0,0 +1,57 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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(coul/diel,PairCoulDiel) + +#else + +#ifndef LMP_PAIR_COUL_DIEL_H +#define LMP_PAIR_COUL_DIEL_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairCoulDiel : public Pair { + public: + PairCoulDiel(class LAMMPS *); + virtual ~PairCoulDiel(); + + virtual void compute(int, int); + + virtual void settings(int, char **); + virtual void coeff(int, char **); + virtual void init_style(); + virtual double init_one(int, int); + virtual void write_restart(FILE *); + virtual void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double cut_global; + double **cut; + double **sigmae, **rme, **offset; + double a_eps, b_eps, eps_s; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/pair_gauss_cut.cpp b/src/USER-MISC/pair_gauss_cut.cpp new file mode 100644 index 0000000000..0568c40ea0 --- /dev/null +++ b/src/USER-MISC/pair_gauss_cut.cpp @@ -0,0 +1,371 @@ +/* ---------------------------------------------------------------------- + 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 authors: Arben Jusufi, Axel Kohlmeyer (Temple U.) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_gauss_cut.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "update.h" +#include "integrate.h" +#include "memory.h" +#include "error.h" +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +PairGaussCut::PairGaussCut(LAMMPS *lmp) : Pair(lmp) +{ + respa_enable = 0; +} + +/* ---------------------------------------------------------------------- */ + +PairGaussCut::~PairGaussCut() +{ + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(hgauss); + memory->destroy(sigmah); + memory->destroy(rmh); + memory->destroy(pgauss); + memory->destroy(offset); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairGaussCut::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,r,rexp,ugauss,factor_lj; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + jtype = type[j]; + + if (rsq < cutsq[itype][jtype]) { + r = sqrt(rsq); + rexp = (r-rmh[itype][jtype])/sigmah[itype][jtype]; + ugauss = pgauss[itype][jtype]*exp(-0.5*rexp*rexp); + fpair = factor_lj*rexp/r*ugauss/sigmah[itype][jtype]; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) { + evdwl = ugauss - offset[itype][jtype]; + evdwl *= factor_lj; + } + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairGaussCut::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(hgauss,n+1,n+1,"pair:hgauss"); + memory->create(sigmah,n+1,n+1,"pair:sigmah"); + memory->create(rmh,n+1,n+1,"pair:rmh"); + memory->create(pgauss,n+1,n+1,"pair:pgauss"); + memory->create(offset,n+1,n+1,"pair:offset"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairGaussCut::settings(int narg, char **arg) +{ + if (narg != 1) error->all(FLERR,"Illegal pair_style command"); + + cut_global = force->numeric(arg[0]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGaussCut::coeff(int narg, char **arg) +{ + if (narg < 5 || narg > 6) error->all(FLERR,"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 hgauss_one = force->numeric(arg[2]); + double rmh_one = force->numeric(arg[3]); + double sigmah_one = force->numeric(arg[4]); + + double cut_one = cut_global; + if (narg == 6) cut_one = force->numeric(arg[5]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + hgauss[i][j] = hgauss_one; + sigmah[i][j] = sigmah_one; + rmh[i][j] = rmh_one; + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairGaussCut::init_one(int i, int j) +{ + + if (setflag[i][j] == 0) { + error->all(FLERR,"for gauss/cut pair style, parameters need to be set explicitly for all pairs."); + } + pgauss[i][j] = hgauss[i][j] / sqrt(MY_2PI)/ sigmah[i][j]; + + + if (offset_flag) { + double rexp = (cut[i][j]-rmh[i][j])/sigmah[i][j]; + offset[i][j] = pgauss[i][j] * exp(-0.5*rexp*rexp); + } else offset[i][j] = 0.0; + + hgauss[j][i] = hgauss[i][j]; + sigmah[j][i] = sigmah[i][j]; + rmh[j][i] = rmh[i][j]; + pgauss[j][i] = pgauss[i][j]; + offset[j][i] = offset[i][j]; + cut[j][i] = cut[i][j]; + + // compute I,J contribution to long-range tail correction + // count total # of atoms of type I and J via Allreduce + + if (tail_flag) { + int *type = atom->type; + int nlocal = atom->nlocal; + + double count[2],all[2]; + count[0] = count[1] = 0.0; + for (int k = 0; k < nlocal; k++) { + if (type[k] == i) count[0] += 1.0; + if (type[k] == j) count[1] += 1.0; + } + MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); + } + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairGaussCut::write_restart(FILE *fp) +{ + write_restart_settings(fp); + + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&hgauss[i][j],sizeof(double),1,fp); + fwrite(&rmh[i][j],sizeof(double),1,fp); + fwrite(&sigmah[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairGaussCut::read_restart(FILE *fp) +{ + read_restart_settings(fp); + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&hgauss[i][j],sizeof(double),1,fp); + fread(&rmh[i][j],sizeof(double),1,fp); + fread(&sigmah[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&hgauss[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&rmh[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&sigmah[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairGaussCut::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&offset_flag,sizeof(int),1,fp); + fwrite(&mix_flag,sizeof(int),1,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairGaussCut::read_restart_settings(FILE *fp) +{ + int me = comm->me; + if (me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairGaussCut::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &fforce) +{ + double r, rexp,ugauss,phigauss; + + r=sqrt(rsq); + rexp = (r-rmh[itype][jtype])/sigmah[itype][jtype]; + ugauss = pgauss[itype][jtype]*exp(-0.5*rexp*rexp); + + fforce = factor_lj*rexp/r*ugauss/sigmah[itype][jtype]; + + phigauss = ugauss - offset[itype][jtype]; + return factor_lj*phigauss; +} + +/* ---------------------------------------------------------------------- */ +double PairGaussCut::memory_usage() +{ + const int n=atom->ntypes; + + double bytes = Pair::memory_usage(); + + bytes += 7*((n+1)*(n+1) * sizeof(double) + (n+1)*sizeof(double *)); + bytes += 1*((n+1)*(n+1) * sizeof(int) + (n+1)*sizeof(int *)); + + return bytes; +} + diff --git a/src/USER-MISC/pair_gauss_cut.h b/src/USER-MISC/pair_gauss_cut.h new file mode 100644 index 0000000000..8a10308278 --- /dev/null +++ b/src/USER-MISC/pair_gauss_cut.h @@ -0,0 +1,60 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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(gauss/cut,PairGaussCut) + +#else + +#ifndef LMP_PAIR_GAUSS_CUT_H +#define LMP_PAIR_GAUSS_CUT_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairGaussCut : public Pair { + public: + PairGaussCut(class LAMMPS *); + ~PairGaussCut(); + + virtual void compute(int, int); + + virtual double single(int, int, int, int, double, double, double, double &); + + virtual void settings(int, char **); + virtual void coeff(int, char **); + + virtual double init_one(int, int); + + virtual void write_restart(FILE *); + virtual void read_restart(FILE *); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + + virtual double memory_usage(); + + protected: + double cut_global; + double **cut; + double **hgauss,**sigmah,**rmh; + double **pgauss,**offset; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp new file mode 100644 index 0000000000..25a8f32fa5 --- /dev/null +++ b/src/USER-MISC/pair_tersoff_table.cpp @@ -0,0 +1,1007 @@ +/* ---------------------------------------------------------------------- + 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: Luca Ferraro (CASPUR) + email: luca.ferraro@caspur.it + + Tersoff Potential + References: (referenced as tersoff_2 functional form in LAMMPS manual) + 1) Tersoff, Phys. Rev. B 39, 5566 (1988) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_tersoff_table.h" +#include "atom.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 + +/* ---------------------------------------------------------------------- */ + +PairTersoffTable::PairTersoffTable(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + one_coeff = 1; + + nelements = 0; + elements = NULL; + nparams = maxparam = 0; + params = NULL; + elem2param = NULL; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairTersoffTable::~PairTersoffTable() +{ + if (elements) + for (int i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + memory->destroy(params); + memory->destroy(elem2param); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + delete [] map; + + deallocateGrids(); + deallocatePreLoops(); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffTable::compute(int eflag, int vflag) +{ + int i,j,k,ii,inum,jnum; + int itype,jtype,ktype,ijparam,ikparam,ijkparam; + double xtmp,ytmp,ztmp; + double fxtmp,fytmp,fztmp; + int *ilist,*jlist,*numneigh,**firstneigh; + + + int interpolIDX; + double r_ik_x, r_ik_y, r_ik_z; + double directorCos_ij_x, directorCos_ij_y, directorCos_ij_z, directorCos_ik_x, directorCos_ik_y, directorCos_ik_z; + double invR_ij, invR_ik, cosTeta; + double repulsivePotential, attractivePotential; + double exponentRepulsivePotential, exponentAttractivePotential,interpolTMP,interpolDeltaX,interpolY1; + double interpolY2, cutoffFunctionIJ, attractiveExponential, repulsiveExponential, cutoffFunctionDerivedIJ,zeta; + double gtetaFunctionIJK,gtetaFunctionDerivedIJK,cutoffFunctionIK; + double cutoffFunctionDerivedIK,factor_force3_ij,factor_1_force3_ik; + double factor_2_force3_ik,betaZetaPowerIJK,betaZetaPowerDerivedIJK,factor_force_tot; + double factor_force_ij; + double gtetaFunctionDerived_temp,gtetaFunction_temp; + + double evdwl = 0.0; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over full neighbor list of my atoms + for (ii = 0; ii < inum; ii++) { + + i = ilist[ii]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + fxtmp = fytmp = fztmp = 0.0; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // Pre-calculate gteta and cutoff function + for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { + + double dr_ij[3], r_ij; + + j = jlist[neighbor_j]; + j &= NEIGHMASK; + + dr_ij[0] = xtmp - x[j][0]; + dr_ij[1] = ytmp - x[j][1]; + dr_ij[2] = ztmp - x[j][2]; + r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + + jtype = map[type[j]]; + ijparam = elem2param[itype][jtype][jtype]; + + if (r_ij > params[ijparam].cutsq) continue; + + r_ij = sqrt(r_ij); + + invR_ij = 1.0 / r_ij; + + directorCos_ij_x = invR_ij * dr_ij[0]; + directorCos_ij_y = invR_ij * dr_ij[1]; + directorCos_ij_z = invR_ij * dr_ij[2]; + + // preCutoffFunction + interpolDeltaX = r_ij - GRIDSTART; + interpolTMP = (interpolDeltaX * GRIDDENSITY_FCUTOFF); + interpolIDX = (int) interpolTMP; + interpolY1 = cutoffFunction[itype][jtype][interpolIDX]; + interpolY2 = cutoffFunction[itype][jtype][interpolIDX+1]; + preCutoffFunction[neighbor_j] = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + // preCutoffFunctionDerived + interpolY1 = cutoffFunctionDerived[itype][jtype][interpolIDX]; + interpolY2 = cutoffFunctionDerived[itype][jtype][interpolIDX+1]; + preCutoffFunctionDerived[neighbor_j] = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + + + for (int neighbor_k = neighbor_j + 1; neighbor_k < jnum; neighbor_k++) { + double dr_ik[3], r_ik; + + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + dr_ik[0] = xtmp -x[k][0]; + dr_ik[1] = ytmp -x[k][1]; + dr_ik[2] = ztmp -x[k][2]; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + + if (r_ik > params[ikparam].cutsq) continue; + + r_ik = sqrt(r_ik); + + invR_ik = 1.0 / r_ik; + + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; + + cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + directorCos_ij_z * directorCos_ik_z; + + // preGtetaFunction + interpolDeltaX=cosTeta+1.0; + interpolTMP = (interpolDeltaX * GRIDDENSITY_GTETA); + interpolIDX = (int) interpolTMP; + interpolY1 = gtetaFunction[itype][interpolIDX]; + interpolY2 = gtetaFunction[itype][interpolIDX+1]; + gtetaFunction_temp = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + // preGtetaFunctionDerived + interpolY1 = gtetaFunctionDerived[itype][interpolIDX]; + interpolY2 = gtetaFunctionDerived[itype][interpolIDX+1]; + gtetaFunctionDerived_temp = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + + preGtetaFunction[neighbor_j][neighbor_k]=params[ijkparam].gamma*gtetaFunction_temp; + preGtetaFunctionDerived[neighbor_j][neighbor_k]=params[ijkparam].gamma*gtetaFunctionDerived_temp; + preGtetaFunction[neighbor_k][neighbor_j]=params[ijkparam].gamma*gtetaFunction_temp; + preGtetaFunctionDerived[neighbor_k][neighbor_j]=params[ijkparam].gamma*gtetaFunctionDerived_temp; + + } // loop on K + + } // loop on J + + + // loop over neighbours of atom i + for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { + + double dr_ij[3], r_ij, f_ij[3]; + + j = jlist[neighbor_j]; + j &= NEIGHMASK; + + dr_ij[0] = xtmp - x[j][0]; + dr_ij[1] = ytmp - x[j][1]; + dr_ij[2] = ztmp - x[j][2]; + r_ij = dr_ij[0]*dr_ij[0] + dr_ij[1]*dr_ij[1] + dr_ij[2]*dr_ij[2]; + + jtype = map[type[j]]; + ijparam = elem2param[itype][jtype][jtype]; + + if (r_ij > params[ijparam].cutsq) continue; + + r_ij = sqrt(r_ij); + invR_ij = 1.0 / r_ij; + + directorCos_ij_x = invR_ij * dr_ij[0]; + directorCos_ij_y = invR_ij * dr_ij[1]; + directorCos_ij_z = invR_ij * dr_ij[2]; + + exponentRepulsivePotential = params[ijparam].lam1 * r_ij; + exponentAttractivePotential = params[ijparam].lam2 * r_ij; + + // repulsiveExponential + interpolDeltaX = exponentRepulsivePotential - minArgumentExponential; + interpolTMP = (interpolDeltaX * GRIDDENSITY_EXP); + interpolIDX = (int) interpolTMP; + interpolY1 = exponential[interpolIDX]; + interpolY2 = exponential[interpolIDX+1]; + repulsiveExponential = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + // attractiveExponential + interpolDeltaX = exponentAttractivePotential - minArgumentExponential; + interpolTMP = (interpolDeltaX * GRIDDENSITY_EXP); + interpolIDX = (int) interpolTMP; + interpolY1 = exponential[interpolIDX]; + interpolY2 = exponential[interpolIDX+1]; + attractiveExponential = interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX); + + repulsivePotential = params[ijparam].biga * repulsiveExponential; + attractivePotential = -params[ijparam].bigb * attractiveExponential; + + cutoffFunctionIJ = preCutoffFunction[neighbor_j]; + cutoffFunctionDerivedIJ = preCutoffFunctionDerived[neighbor_j]; + + zeta = 0.0; + + // first loop over neighbours of atom i except j - part 1/2 + for (int neighbor_k = 0; neighbor_k < neighbor_j; neighbor_k++) { + double dr_ik[3], r_ik; + + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + dr_ik[0] = xtmp -x[k][0]; + dr_ik[1] = ytmp -x[k][1]; + dr_ik[2] = ztmp -x[k][2]; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + + if (r_ik > params[ikparam].cutsq) continue; + + r_ik = sqrt(r_ik); + + invR_ik = 1.0 / r_ik; + + directorCos_ik_x = invR_ik * r_ik_x; + directorCos_ik_y = invR_ik * r_ik_y; + directorCos_ik_z = invR_ik * r_ik_z; + + gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k]; + + cutoffFunctionIK = preCutoffFunction[neighbor_k]; + + zeta += cutoffFunctionIK * gtetaFunctionIJK; + + } + + // first loop over neighbours of atom i except j - part 2/2 + for (int neighbor_k = neighbor_j+1; neighbor_k < jnum; neighbor_k++) { + double dr_ik[3], r_ik; + + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + dr_ik[0] = xtmp -x[k][0]; + dr_ik[1] = ytmp -x[k][1]; + dr_ik[2] = ztmp -x[k][2]; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + + if (r_ik > params[ikparam].cutsq) continue; + + r_ik = sqrt(r_ik); + invR_ik = 1.0 / r_ik; + + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; + + gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k]; + + cutoffFunctionIK = preCutoffFunction[neighbor_k]; + + zeta += cutoffFunctionIK * gtetaFunctionIJK; + } + + // betaZetaPowerIJK + interpolDeltaX= params[ijparam].beta * zeta; + interpolTMP = (interpolDeltaX * GRIDDENSITY_BIJ); + interpolIDX = (int) interpolTMP; + interpolY1 = betaZetaPower[itype][interpolIDX]; + interpolY2 = betaZetaPower[itype][interpolIDX+1]; + betaZetaPowerIJK = (interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX)); + // betaZetaPowerDerivedIJK + interpolY1 = betaZetaPowerDerived[itype][interpolIDX]; + interpolY2 = betaZetaPowerDerived[itype][interpolIDX+1]; + betaZetaPowerDerivedIJK = params[ijparam].beta*(interpolY1 + (interpolY2 - interpolY1) * (interpolTMP - interpolIDX)); + + // Forces and virial + factor_force_ij = 0.5*cutoffFunctionDerivedIJ*(repulsivePotential + attractivePotential * betaZetaPowerIJK)+0.5*cutoffFunctionIJ*(-repulsivePotential*params[ijparam].lam1-betaZetaPowerIJK*attractivePotential*params[ijparam].lam2); + + f_ij[0] = factor_force_ij * directorCos_ij_x; + f_ij[1] = factor_force_ij * directorCos_ij_y; + f_ij[2] = factor_force_ij * directorCos_ij_z; + + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; + + fxtmp -= f_ij[0]; + fytmp -= f_ij[1]; + fztmp -= f_ij[2]; + + // potential energy + evdwl = cutoffFunctionIJ * repulsivePotential + cutoffFunctionIJ * attractivePotential * betaZetaPowerIJK; + + if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.5 * evdwl, 0.0, + -factor_force_ij*invR_ij, dr_ij[0], dr_ij[1], dr_ij[2]); + + factor_force_tot= 0.5*cutoffFunctionIJ*attractivePotential*betaZetaPowerDerivedIJK; + + // second loop over neighbours of atom i except j, forces and virial only - part 1/2 + for (int neighbor_k = 0; neighbor_k < neighbor_j; neighbor_k++) { + double dr_ik[3], r_ik, f_ik[3]; + + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + dr_ik[0] = xtmp -x[k][0]; + dr_ik[1] = ytmp -x[k][1]; + dr_ik[2] = ztmp -x[k][2]; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + + if (r_ik > params[ikparam].cutsq) continue; + + r_ik = sqrt(r_ik); + invR_ik = 1.0 / r_ik; + + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; + + cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + directorCos_ij_z * directorCos_ik_z; + + gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k]; + + gtetaFunctionDerivedIJK = preGtetaFunctionDerived[neighbor_j][neighbor_k]; + + cutoffFunctionIK = preCutoffFunction[neighbor_k]; + + cutoffFunctionDerivedIK = preCutoffFunctionDerived[neighbor_k]; + + factor_force3_ij= cutoffFunctionIK * gtetaFunctionDerivedIJK * invR_ij *factor_force_tot; + + f_ij[0] = factor_force3_ij * (directorCos_ij_x*cosTeta - directorCos_ik_x); + f_ij[1] = factor_force3_ij * (directorCos_ij_y*cosTeta - directorCos_ik_y); + f_ij[2] = factor_force3_ij * (directorCos_ij_z*cosTeta - directorCos_ik_z); + + factor_1_force3_ik = (cutoffFunctionIK * gtetaFunctionDerivedIJK * invR_ik)*factor_force_tot; + factor_2_force3_ik = -(cutoffFunctionDerivedIK * gtetaFunctionIJK)*factor_force_tot; + + f_ik[0] = factor_1_force3_ik * (directorCos_ik_x*cosTeta - directorCos_ij_x) + factor_2_force3_ik * directorCos_ik_x; + f_ik[1] = factor_1_force3_ik * (directorCos_ik_y*cosTeta - directorCos_ij_y) + factor_2_force3_ik * directorCos_ik_y; + f_ik[2] = factor_1_force3_ik * (directorCos_ik_z*cosTeta - directorCos_ij_z) + factor_2_force3_ik * directorCos_ik_z; + + f[j][0] -= f_ij[0]; + f[j][1] -= f_ij[1]; + f[j][2] -= f_ij[2]; + + f[k][0] -= f_ik[0]; + f[k][1] -= f_ik[1]; + f[k][2] -= f_ik[2]; + + fxtmp += f_ij[0] + f_ik[0]; + fytmp += f_ij[1] + f_ik[1]; + fztmp += f_ij[2] + f_ik[2]; + + // potential energy + evdwl = 0.0; + + if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + } + + // second loop over neighbours of atom i except j, forces and virial only - part 2/2 + for (int neighbor_k = neighbor_j+1; neighbor_k < jnum; neighbor_k++) { + double dr_ik[3], r_ik, f_ik[3]; + + k = jlist[neighbor_k]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + dr_ik[0] = xtmp -x[k][0]; + dr_ik[1] = ytmp -x[k][1]; + dr_ik[2] = ztmp -x[k][2]; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + + if (r_ik > params[ikparam].cutsq) continue; + + r_ik = sqrt(r_ik); + invR_ik = 1.0 / r_ik; + + directorCos_ik_x = invR_ik * dr_ik[0]; + directorCos_ik_y = invR_ik * dr_ik[1]; + directorCos_ik_z = invR_ik * dr_ik[2]; + + cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + directorCos_ij_z * directorCos_ik_z; + + gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k]; + + gtetaFunctionDerivedIJK = preGtetaFunctionDerived[neighbor_j][neighbor_k]; + + cutoffFunctionIK = preCutoffFunction[neighbor_k]; + + cutoffFunctionDerivedIK = preCutoffFunctionDerived[neighbor_k]; + + factor_force3_ij= cutoffFunctionIK * gtetaFunctionDerivedIJK * invR_ij *factor_force_tot; + + f_ij[0] = factor_force3_ij * (directorCos_ij_x*cosTeta - directorCos_ik_x); + f_ij[1] = factor_force3_ij * (directorCos_ij_y*cosTeta - directorCos_ik_y); + f_ij[2] = factor_force3_ij * (directorCos_ij_z*cosTeta - directorCos_ik_z); + + factor_1_force3_ik = (cutoffFunctionIK * gtetaFunctionDerivedIJK * invR_ik)*factor_force_tot; + factor_2_force3_ik = -(cutoffFunctionDerivedIK * gtetaFunctionIJK)*factor_force_tot; + + f_ik[0] = factor_1_force3_ik * (directorCos_ik_x*cosTeta - directorCos_ij_x) + factor_2_force3_ik * directorCos_ik_x; + f_ik[1] = factor_1_force3_ik * (directorCos_ik_y*cosTeta - directorCos_ij_y) + factor_2_force3_ik * directorCos_ik_y; + f_ik[2] = factor_1_force3_ik * (directorCos_ik_z*cosTeta - directorCos_ij_z) + factor_2_force3_ik * directorCos_ik_z; + + f[j][0] -= f_ij[0]; + f[j][1] -= f_ij[1]; + f[j][2] -= f_ij[2]; + + f[k][0] -= f_ik[0]; + f[k][1] -= f_ik[1]; + f[k][2] -= f_ik[2]; + + fxtmp += f_ij[0] + f_ik[0]; + fytmp += f_ij[1] + f_ik[1]; + fztmp += f_ij[2] + f_ik[2]; + + // potential energy + evdwl = 0.0; + + if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + + } + } // loop on J + f[i][0] += fxtmp; + f[i][1] += fytmp; + f[i][2] += fztmp; + } // loop on I + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffTable::deallocatePreLoops(void) +{ + memory->destroy (preGtetaFunction); + memory->destroy (preGtetaFunctionDerived); + memory->destroy(preCutoffFunction); + memory->destroy(preCutoffFunctionDerived); +} + +void PairTersoffTable::allocatePreLoops(void) +{ + memory->create(preGtetaFunction,leadingDimensionInteractionList,leadingDimensionInteractionList,"tersofftable:preGtetaFunction"); + + memory->create(preGtetaFunctionDerived,leadingDimensionInteractionList,leadingDimensionInteractionList,"tersofftable:preGtetaFunctionDerived"); + + memory->create(preCutoffFunction,leadingDimensionInteractionList,"tersofftable:preCutoffFunction"); + + memory->create(preCutoffFunctionDerived,leadingDimensionInteractionList,"tersofftable:preCutoffFunctionDerived"); +} + +void PairTersoffTable::deallocateGrids() +{ + int i,j; + + memory->destroy(exponential); + memory->destroy(gtetaFunction); + memory->destroy(gtetaFunctionDerived); + memory->destroy(cutoffFunction); + memory->destroy(cutoffFunctionDerived); + memory->destroy(betaZetaPower); + memory->destroy(betaZetaPowerDerived); +} + +void PairTersoffTable::allocateGrids(void) +{ + int i, j, l; + + int numGridPointsExponential, numGridPointsGtetaFunction, numGridPointsOneCutoffFunction; + int numGridPointsNotOneCutoffFunction, numGridPointsCutoffFunction, numGridPointsBetaZetaPower; + // double minArgumentExponential; + double deltaArgumentCutoffFunction, deltaArgumentExponential, deltaArgumentBetaZetaPower; + double deltaArgumentGtetaFunction; + double r, minMu, maxLambda, maxCutoff; + double const PI=acos(-1.0); + + // exponential + + // find min and max argument + minMu=params[0].lam2; + maxLambda=params[0].lam1; + for (i=1; i maxLambda) maxLambda = params[i].lam1; + } + maxCutoff=cutmax; + + minArgumentExponential=minMu*GRIDSTART; + + numGridPointsExponential=(int)((maxLambda*maxCutoff-minArgumentExponential)*GRIDDENSITY_EXP)+2; + + memory->create(exponential,numGridPointsExponential,"tersofftable:exponential"); + + r = minArgumentExponential; + deltaArgumentExponential = 1.0 / GRIDDENSITY_EXP; + for (i = 0; i < numGridPointsExponential; i++) + { + exponential[i] = exp(-r); + r += deltaArgumentExponential; + } + + + // gtetaFunction + + numGridPointsGtetaFunction=(int)(2.0*GRIDDENSITY_GTETA)+2; + + memory->create(gtetaFunction,nelements,numGridPointsGtetaFunction,"tersofftable:gtetaFunction"); + memory->create(gtetaFunctionDerived,nelements,numGridPointsGtetaFunction,"tersofftable:gtetaFunctionDerived"); + + r = minArgumentExponential; + for (i=0; icreate(cutoffFunction,nelements,nelements,ngrid_max,"tersoff:cutfunc"); + memory->create(cutoffFunctionDerived,nelements,nelements,ngrid_max,"tersoff:cutfuncD"); + + // cutoffFunction, compute. + + for (i=0; icreate(betaZetaPower,nelements,zeta_max,"tersoff:zetafunc"); + memory->create(betaZetaPowerDerived,nelements,zeta_max,"tersoff:zetafuncD"); + + for (i=0; intypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + + map = new int[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairTersoffTable::settings(int narg, char **arg) +{ + if (narg != 0) error->all(FLERR,"Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairTersoffTable::coeff(int narg, char **arg) +{ + int i,j,n; + + if (!allocated) allocate(); + + if (narg != 3 + atom->ntypes) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // insure I,J args are * * + + if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // read args that map atom types to elements in potential file + // map[i] = which element the Ith atom type is, -1 if NULL + // nelements = # of unique elements + // elements = list of element names + + if (elements) { + for (i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + } + elements = new char*[atom->ntypes]; + for (i = 0; i < atom->ntypes; i++) elements[i] = NULL; + + nelements = 0; + for (i = 3; i < narg; i++) { + if (strcmp(arg[i],"NULL") == 0) { + map[i-2] = -1; + continue; + } + for (j = 0; j < nelements; j++) + if (strcmp(arg[i],elements[j]) == 0) break; + map[i-2] = j; + if (j == nelements) { + n = strlen(arg[i]) + 1; + elements[j] = new char[n]; + strcpy(elements[j],arg[i]); + nelements++; + } + } + + // read potential file and initialize potential parameters + + read_file(arg[2]); + setup(); + + // clear setflag since coeff() called once with I,J = * * + + n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + count++; + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + + // allocate tables and internal structures + allocatePreLoops(); + allocateGrids(); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairTersoffTable::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style Tersoff requires atom IDs"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style Tersoff requires newton pair on"); + + // need a full neighbor list + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairTersoffTable::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + return cutmax; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffTable::read_file(char *file) +{ + int params_per_line = 17; + char **words = new char*[params_per_line+1]; + + memory->sfree(params); + params = NULL; + nparams = maxparam = 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(FLERR,str); + } + } + + // read each set of params from potential file + // one set of params can span multiple lines + // store 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(FLERR,"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 entry in file + + 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]); // not used (only tersoff_2 is implemented) + params[nparams].gamma = atof(words[4]); // not used (only tersoff_2 is implemented) + params[nparams].lam3 = atof(words[5]); // not used (only tersoff_2 is implemented) + 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]); + // current implementation is based on functional form + // of tersoff_2 as reported in the reference paper + double bigr = atof(words[13]); + double bigd = atof(words[14]); + params[nparams].cutoffR = bigr - bigd; + params[nparams].cutoffS = bigr + bigd; + params[nparams].lam1 = atof(words[15]); + params[nparams].biga = atof(words[16]); + + // currently only allow m exponent of 1 or 3 + params[nparams].powermint = int(params[nparams].powerm); + + if (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].cutoffR < 0.0 ||params[nparams].cutoffS < 0.0 || + params[nparams].cutoffR > params[nparams].cutoffS || + params[nparams].lam1 < 0.0 || params[nparams].biga < 0.0 + ) error->all(FLERR,"Illegal Tersoff parameter"); + + // only tersoff_2 parametrization is implemented + if (params[nparams].gamma != 1.0 || params[nparams].lam3 != 0.0) + error->all(FLERR,"Current tersoff/table pair_style implements only tersoff_2 parametrization"); + nparams++; + } + + delete [] words; +} + +/* ---------------------------------------------------------------------- */ + +void PairTersoffTable::setup() +{ + int i,j,k,m,n; + + // set elem2param for all triplet combinations + // must be a single exact match to lines read from file + // do not allow for ACB in place of ABC + + memory->destroy(elem2param); + memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param"); + + for (i = 0; i < nelements; i++) + for (j = 0; j < nelements; j++) + for (k = 0; k < nelements; k++) { + n = -1; + for (m = 0; m < nparams; m++) { + if (i == params[m].ielement && j == params[m].jelement && + k == params[m].kelement) { + if (n >= 0) error->all(FLERR,"Potential file has duplicate entry"); + n = m; + } + } + if (n < 0) error->all(FLERR,"Potential file is missing an entry"); + elem2param[i][j][k] = n; + } + + // set cutoff square + for (m = 0; m < nparams; m++) { + params[m].cut = params[m].cutoffS; + params[m].cutsq = params[m].cut*params[m].cut; + } + + // set cutmax to max of all params + cutmax = 0.0; + for (m = 0; m < nparams; m++) { + if (params[m].cut > cutmax) cutmax = params[m].cut; + } +} diff --git a/src/USER-MISC/pair_tersoff_table.h b/src/USER-MISC/pair_tersoff_table.h new file mode 100644 index 0000000000..3be63cde0f --- /dev/null +++ b/src/USER-MISC/pair_tersoff_table.h @@ -0,0 +1,96 @@ +/* ---------------------------------------------------------------------- + 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: Luca Ferraro (CASPUR) + email: luca.ferraro@caspur.it + + Tersoff Potential + References: + 1) Tersoff, Phys. Rev. B 39, 5566 (1988) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(tersoff/table,PairTersoffTable) + +#else + +#ifndef LMP_PAIR_Tersoff_H +#define LMP_PAIR_Tersoff_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairTersoffTable : public Pair { + public: + PairTersoffTable(class LAMMPS *); + virtual ~PairTersoffTable(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + + protected: + struct Param { + double lam1,lam2,lam3; + double c,d,h; + double gamma,powerm; + double powern,beta; + double biga,bigb,cutoffR, cutoffS; + double cut,cutsq; + int ielement,jelement,kelement; + int powermint; + }; + + + double cutmax; // max cutoff for all elements + int nelements; // # of unique elements + char **elements; // names 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 + Param *params; // parameter set for an I-J-K interaction + + void allocate(); + + void read_file(char *); + void setup(); + + // max number of interaction per atom for environment potential + static const int leadingDimensionInteractionList = 64; + // pre-loop coordination functions + double **preGtetaFunction, **preGtetaFunctionDerived; + double *preCutoffFunction, *preCutoffFunctionDerived; + void allocatePreLoops(void); + void deallocatePreLoops(void); + // grids + static const double GRIDSTART = 0.1; + static const int GRIDDENSITY_FCUTOFF = 5000; + static const int GRIDDENSITY_EXP = 12000; + static const int GRIDDENSITY_GTETA = 12000; + static const int GRIDDENSITY_BIJ = 7500; + double minArgumentExponential; + double *exponential, ***cutoffFunction, ***cutoffFunctionDerived; + double **gtetaFunction, **gtetaFunctionDerived, **betaZetaPower, **betaZetaPowerDerived; + void allocateGrids(void); + void deallocateGrids(void); +}; + +} + +#endif +#endif