jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2015-09-08 23:55:07 +00:00
parent fafba07dff
commit 9aac00e7ed
6 changed files with 70 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
src/CORESHELL/pair_born_coul_long_cs.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Hendrik Heenen (hendrik.heenen@mytum.de)
Added epsilon for use with CORESHELL and USER-DRUDE
------------------------------------------------------------------------- */
#include "math.h"
@ -34,12 +34,13 @@ using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
#define EWALD_P 9.95473818e-1
#define B0 -0.1335096380159268
#define B1 -2.57839507e-1
#define B2 -1.37203639e-1
#define B3 -8.88822059e-3
#define B4 -5.80844129e-3
#define B5 1.14652755e-1
#define EPSILON 1.0e-20
#define EPS_EWALD 1.0e-6
@ -62,7 +63,7 @@ void PairBornCoulLongCS::compute(int eflag, int vflag)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double fraction,table;
double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
double grij,expm2,prefactor,t,erfc,u;
double r,rexp;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -123,7 +124,8 @@ void PairBornCoulLongCS::compute(int eflag, int vflag)
grij = g_ewald * (r+EPS_EWALD);
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= (r+EPS_EWALD);
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
// Additionally r2inv needs to be accordingly modified since the later
@ -133,7 +135,8 @@ void PairBornCoulLongCS::compute(int eflag, int vflag)
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
}

2
src/CORESHELL/pair_born_coul_long_cs.h
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Hendrik Heenen (hendrik.heenen@mytum.com)
Added epsilon for use with CORESHELL and USER-DRUDE
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS

23
src/CORESHELL/pair_buck_coul_long_cs.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Hendrik Heenen (hendrik.heenen@mytum.de)
Added epsilon for use with CORESHELL and USER-DRUDE
------------------------------------------------------------------------- */
#include "math.h"
@ -34,12 +34,13 @@ using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
#define EWALD_P 9.95473818e-1
#define B0 -0.1335096380159268
#define B1 -2.57839507e-1
#define B2 -1.37203639e-1
#define B3 -8.88822059e-3
#define B4 -5.80844129e-3
#define B5 1.14652755e-1
#define EPSILON 1.0e-20
#define EPS_EWALD 1.0e-6
@ -62,7 +63,7 @@ void PairBuckCoulLongCS::compute(int eflag, int vflag)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double fraction,table;
double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
double grij,expm2,prefactor,t,erfc,u;
double r,rexp;
int *ilist,*jlist,*numneigh,**firstneigh;
@ -123,7 +124,8 @@ void PairBuckCoulLongCS::compute(int eflag, int vflag)
grij = g_ewald * (r+EPS_EWALD);
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= (r+EPS_EWALD);
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
// Additionally r2inv needs to be accordingly modified since the later
@ -133,7 +135,8 @@ void PairBuckCoulLongCS::compute(int eflag, int vflag)
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
}

2
src/CORESHELL/pair_buck_coul_long_cs.h
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Hendrik Heenen (hendrik.heenen@mytum.com)
Added epsilon for use with CORESHELL and USER-DRUDE
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS

46
src/CORESHELL/pair_coul_long_cs.cpp
View File

@ -43,9 +43,18 @@ using namespace LAMMPS_NS;
#define B4 -5.80844129e-3
#define B5 1.14652755e-1
#define EPSILON 1.0e-20
#define EPS_EWALD 1.0e-6
#define EPS_EWALD_SQR 1.0e-12
/* ---------------------------------------------------------------------- */
PairCoulLongCS::PairCoulLongCS(LAMMPS *lmp) : PairCoulLong(lmp) {}
PairCoulLongCS::PairCoulLongCS(LAMMPS *lmp) : PairCoulLong(lmp)
{
ewaldflag = pppmflag = 1;
ftable = NULL;
qdist = 0.0;
}
/* ---------------------------------------------------------------------- */
@ -101,18 +110,34 @@ void PairCoulLongCS::compute(int eflag, int vflag)
jtype = type[j];
if (rsq < cut_coulsq) {
rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
r2inv = 1.0/rsq;
r = sqrt(rsq);
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j];
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
if (factor_coul < 1.0) {
// When bonded parts are being calculated a minimal distance (EPS_EWALD)
// has to be added to the prefactor and erfc in order to make the
// used approximation functions for the Ewald correction valid
grij = g_ewald * (r+EPS_EWALD);
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= (r+EPS_EWALD);
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
// Additionally r2inv needs to be accordingly modified since the later
// scaling of the overall force shall be consistent
r2inv = 1.0/(rsq + EPS_EWALD_SQR);
} else {
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
u = 1. - t;
erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
prefactor /= r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
}
} else {
union_int_float_t rsq_lookup;
rsq_lookup.f = rsq;
@ -157,3 +182,4 @@ void PairCoulLongCS::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
}

8
src/CORESHELL/pair_coul_long_cs.h
View File

@ -11,9 +11,13 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Added epsilon for use with CORESHELL and USER-DRUDE
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(coul/long/cs,PairCoulLongCS)
PairStyle(coul/long/cs/,PairCoulLongCS)
#else
@ -47,7 +51,7 @@ E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style lj/cut/coul/long/drude requires atom attribute q
E: Pair style lj/cut/coul/long requires atom attribute q
The atom style defined does not have this attribute.