ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2011-10-10 17:03:52 +00:00
parent 2e2abb3a72
commit e6a0e39b50
7 changed files with 70 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/KSPACE/pppm.cpp
View File

@ -35,7 +35,10 @@
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXORDER 7
#define OFFSET 16384
@ -58,7 +61,6 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command");
precision = atof(arg[0]);
PI = 4.0*atan(1.0);
nfactors = 3;
factors = new int[nfactors];
@ -520,9 +522,9 @@ void PPPM::setup()
delvolinv = delxinv*delyinv*delzinv;
double unitkx = (2.0*PI/xprd);
double unitky = (2.0*PI/yprd);
double unitkz = (2.0*PI/zprd_slab);
double unitkx = (2.0*MY_PI/xprd);
double unitky = (2.0*MY_PI/yprd);
double unitkz = (2.0*MY_PI/zprd_slab);
// fkx,fky,fkz for my FFT grid pts
@ -581,11 +583,11 @@ void PPPM::setup()
double sum1,dot1,dot2;
double numerator,denominator;
int nbx = static_cast<int> ((g_ewald*xprd/(PI*nx_pppm)) *
int nbx = static_cast<int> ((g_ewald*xprd/(MY_PI*nx_pppm)) *
pow(-log(EPS_HOC),0.25));
int nby = static_cast<int> ((g_ewald*yprd/(PI*ny_pppm)) *
int nby = static_cast<int> ((g_ewald*yprd/(MY_PI*ny_pppm)) *
pow(-log(EPS_HOC),0.25));
int nbz = static_cast<int> ((g_ewald*zprd_slab/(PI*nz_pppm)) *
int nbz = static_cast<int> ((g_ewald*zprd_slab/(MY_PI*nz_pppm)) *
pow(-log(EPS_HOC),0.25));
double form = 1.0;
@ -708,7 +710,7 @@ void PPPM::compute(int eflag, int vflag)
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 +
0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e*scale;
}
@ -1027,7 +1029,7 @@ double PPPM::rms(double h, double prd, bigint natoms,
for (int m = 0; m < order; m++)
sum += acons[order][m] * pow(h*g_ewald,2.0*m);
double value = q2 * pow(h*g_ewald,order) *
sqrt(g_ewald*prd*sqrt(2.0*PI)*sum/natoms) / (prd*prd);
sqrt(g_ewald*prd*sqrt(2.0*MY_PI)*sum/natoms) / (prd*prd);
return value;
}
@ -1873,13 +1875,13 @@ void PPPM::slabcorr(int eflag)
// compute corrections
double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
if (eflag) energy += qqrd2e*scale * e_slabcorr;
// add on force corrections
double ffact = -4.0*PI*dipole_all/volume;
double ffact = -4.0*MY_PI*dipole_all/volume;
double **f = atom->f;
for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact;

1
src/KSPACE/pppm.h
View File

@ -47,7 +47,6 @@ class PPPM : public KSpace {
protected:
int me,nprocs;
double PI;
double precision;
int nfactors;
int *factors;

9
src/KSPACE/pppm_cg.cpp
View File

@ -26,7 +26,10 @@
#include "memory.h"
#include "pppm_cg.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define OFFSET 16384
#define SMALLQ 0.00001
@ -176,7 +179,7 @@ void PPPMCG::compute(int eflag, int vflag)
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 +
0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e*scale;
}
@ -372,13 +375,13 @@ void PPPMCG::slabcorr(int eflag)
// compute corrections
double e_slabcorr = 2.0*PI*dipole_all*dipole_all/volume;
double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume;
if (eflag) energy += qqrd2e*scale * e_slabcorr;
// add on force corrections
double ffact = -4.0*PI*dipole_all/volume * qqrd2e * scale;
double ffact = -4.0*MY_PI*dipole_all/volume * qqrd2e * scale;
double **f = atom->f;
for (int j = 0; j < num_charged; j++) {

34
src/MANYBODY/pair_tersoff.cpp
View File

@ -29,7 +29,10 @@
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
@ -41,10 +44,6 @@ 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);
nelements = 0;
elements = NULL;
nparams = maxparam = 0;
@ -634,7 +633,7 @@ double PairTersoff::ters_fc(double r, Param *param)
if (r < ters_R-ters_D) return 1.0;
if (r > ters_R+ters_D) return 0.0;
return 0.5*(1.0 - sin(PI2*(r - ters_R)/ters_D));
return 0.5*(1.0 - sin(MY_PI2*(r - ters_R)/ters_D));
}
/* ---------------------------------------------------------------------- */
@ -646,7 +645,7 @@ double PairTersoff::ters_fc_d(double r, Param *param)
if (r < ters_R-ters_D) return 0.0;
if (r > ters_R+ters_D) return 0.0;
return -(PI4/ters_D) * cos(PI2*(r - ters_R)/ters_D);
return -(MY_PI4/ters_D) * cos(MY_PI2*(r - ters_R)/ters_D);
}
/* ---------------------------------------------------------------------- */
@ -700,27 +699,6 @@ double PairTersoff::ters_bij_d(double zeta, Param *param)
/* ---------------------------------------------------------------------- */
double PairTersoff::ters_gijk(double costheta, Param *param)
{
double ters_c = param->c;
double ters_d = param->d;
return param->gamma*(1.0 + pow(ters_c/ters_d,2.0) -
pow(ters_c,2.0) / (pow(ters_d,2.0) + pow(param->h - costheta,2.0)));
};
/* ---------------------------------------------------------------------- */
double PairTersoff::ters_gijk_d(double costheta, Param *param)
{
double numerator = -2.0 * pow(param->c,2) * (param->h - costheta);
double denominator = pow(pow(param->d,2.0) +
pow(param->h - costheta,2.0),2.0);
return param->gamma*numerator/denominator;
}
/* ---------------------------------------------------------------------- */
void PairTersoff::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
@ -794,4 +772,4 @@ void PairTersoff::costheta_d(double *rij_hat, double rij,
vec3_scale(1.0/rik,drk,drk);
vec3_add(drj,drk,dri);
vec3_scale(-1.0,dri,dri);
};
}

54
src/MANYBODY/pair_tersoff.h
View File

@ -28,7 +28,7 @@ class PairTersoff : public Pair {
public:
PairTersoff(class LAMMPS *);
virtual ~PairTersoff();
void compute(int, int);
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
@ -49,22 +49,22 @@ class PairTersoff : public Pair {
double ZBLcut,ZBLexpscale;
};
double PI,PI2,PI4;
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
Param *params; // parameter set for an I-J-K interaction
char **elements; // names of unique elements
int ***elem2param; // mapping from element triplets to parameters
int *map; // mapping from atom types to elements
double cutmax; // max cutoff for all elements
int nelements; // # of unique 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();
virtual void read_file(char *);
void setup();
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 &);
virtual void force_zeta(Param *, double, double, double &,
double &, int, double &);
void attractive(Param *, double, double, double, double *, double *,
double *, double *, double *);
@ -74,26 +74,52 @@ class PairTersoff : public Pair {
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 *);
void costheta_d(double *, double, double *, double,
double *, double *, double *);
// vector functions, inline for efficiency
// inlined functions for efficiency
inline double vec3_dot(double *x, double *y) {
inline double ters_gijk(const double costheta,
const Param * const param) const {
const double ters_c = param->c * param->c;
const double ters_d = param->d * param->d;
const double hcth = param->h - costheta;
return param->gamma*(1.0 + ters_c/ters_d - ters_c / (ters_d + hcth*hcth));
}
inline double ters_gijk_d(const double costheta,
const Param * const param) const {
const double ters_c = param->c * param->c;
const double ters_d = param->d * param->d;
const double hcth = param->h - costheta;
const double numerator = -2.0 * ters_c * hcth;
const double denominator = 1.0/(ters_d + hcth*hcth);
return param->gamma*numerator*denominator*denominator;
}
inline double vec3_dot(const double x[3], const double y[3]) const {
return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
}
inline void vec3_add(double *x, double *y, double *z) {
inline void vec3_add(const double x[3], const double y[3],
double * const z) const {
z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2];
}
inline void vec3_scale(double k, double *x, double *y) {
inline void vec3_scale(const double k, const double x[3],
double y[3]) const {
y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2];
}
inline void vec3_scaleadd(double k, double *x, double *y, double *z) {
z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2];
inline void vec3_scaleadd(const double k, const double x[3],
const double y[3], double * const z) const {
z[0] = k*x[0]+y[0];
z[1] = k*x[1]+y[1];
z[2] = k*x[2]+y[2];
}
};

6
src/MANYBODY/pair_tersoff_zbl.cpp
View File

@ -31,7 +31,9 @@
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
@ -224,7 +226,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce,
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 premult = (param->Z_i * param->Z_j * esq)/(4.0*MY_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);
@ -232,7 +234,7 @@ void PairTersoffZBL::repulsive(Param *param, double rsq, double &fforce,
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 fforce_ZBL = premult*-phi/rsq + premult*dphi/r;
double eng_ZBL = premult*(1.0/r)*phi;
// combine two parts with smoothing by Fermi-like function

1
src/USER-EWALDN/math_vector.h
View File

@ -19,6 +19,7 @@
#define LMP_MATH_VECTOR_H
#include "math.h"
#include "string.h"
#define VECTOR_NULL {0, 0, 0}
#define SHAPE_NULL {0, 0, 0, 0, 0, 0}