port code to current LAMMPS style and make it compatible with OpenMP 4.x compilers
This commit is contained in:
Axel Kohlmeyer
@ -13,6 +13,7 @@
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "atom_vec_dielectric.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "citeme.h"
|
||||
|
||||
|
||||
@ -25,6 +25,9 @@ AtomStyle(dielectric,AtomVecDielectric);
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class AtomVecDielectric : public AtomVec {
|
||||
friend class PairLJCutCoulDebyeDielectric;
|
||||
friend class PairLJLongCoulLongDielectric;
|
||||
|
||||
public:
|
||||
AtomVecDielectric(class LAMMPS *);
|
||||
|
||||
@ -32,19 +35,21 @@ class AtomVecDielectric : public AtomVec {
|
||||
void create_atom_post(int);
|
||||
void data_atom_post(int);
|
||||
void unpack_restart_init(int);
|
||||
int property_atom(char *);
|
||||
void pack_property_atom(int, double *, int, int);
|
||||
|
||||
private:
|
||||
int *num_bond,*num_angle,*num_dihedral,*num_improper;
|
||||
int **bond_type,**angle_type,**dihedral_type,**improper_type;
|
||||
protected:
|
||||
int *num_bond, *num_angle, *num_dihedral, *num_improper;
|
||||
int **bond_type, **angle_type, **dihedral_type, **improper_type;
|
||||
int **nspecial;
|
||||
|
||||
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
|
||||
int bond_per_atom, angle_per_atom, dihedral_per_atom, improper_per_atom;
|
||||
|
||||
double **mu;
|
||||
double *area,*ed,*em,*epsilon,*curvature,*q_unscaled;
|
||||
double *area, *ed, *em, *epsilon, *curvature, *q_unscaled;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -44,7 +44,7 @@ class ComputeEfieldAtom : public Compute {
|
||||
double **efield;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -32,13 +32,13 @@ class FixPolarizeBEMGMRES : public Fix {
|
||||
virtual void init();
|
||||
virtual void setup(int);
|
||||
virtual void pre_force(int);
|
||||
int pack_forward_comm(int, int*, double*, int, int*);
|
||||
void unpack_forward_comm(int, int, double*);
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
virtual double compute_vector(int);
|
||||
|
||||
int modify_param(int, char**);
|
||||
int modify_param(int, char **);
|
||||
double memory_usage();
|
||||
void grow_arrays(int);
|
||||
void copy_arrays(int, int, int);
|
||||
@ -53,53 +53,58 @@ class FixPolarizeBEMGMRES : public Fix {
|
||||
|
||||
protected:
|
||||
int nmax;
|
||||
int nevery; // to be invoked every time steps
|
||||
int* tag2mat; // tag2mat[atom->tag[i]] = the index of the atoms in the induced charge arrays from atom tags
|
||||
int* mat2tag; // mat2tag[idx] = the atom tag of the induced charge idx in the induced charge arrays
|
||||
int* induced_charge_idx; // return the index of the atoms in the induced charge arrays
|
||||
int num_induced_charges; // total number of induced charges
|
||||
double* induced_charges; // values of induced charges
|
||||
double* buffer; // buffer of size num_induced_charges
|
||||
double* q_backup; // backup for the real charges
|
||||
int nevery; // to be invoked every time steps
|
||||
int *
|
||||
tag2mat; // tag2mat[atom->tag[i]] = the index of the atoms in the induced charge arrays from atom tags
|
||||
int *
|
||||
mat2tag; // mat2tag[idx] = the atom tag of the induced charge idx in the induced charge arrays
|
||||
int *induced_charge_idx; // return the index of the atoms in the induced charge arrays
|
||||
int num_induced_charges; // total number of induced charges
|
||||
double *induced_charges; // values of induced charges
|
||||
double *buffer; // buffer of size num_induced_charges
|
||||
double *q_backup; // backup for the real charges
|
||||
int allocated;
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
int torqueflag,extraflag;
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
int torqueflag, extraflag;
|
||||
|
||||
void force_clear();
|
||||
double vec_dot(const double*, const double*, int); // dot product between two vectors of length n
|
||||
double vec_dot(const double *, const double *,
|
||||
int); // dot product between two vectors of length n
|
||||
|
||||
private:
|
||||
int mat_dim; // matrix dimension = total number of induced charges
|
||||
int mr; // number of vectors used to span the Krylov space
|
||||
int iterations; // actual number of iterations
|
||||
int itr_max; // maximum number of outer iterations
|
||||
int randomized; // 1 if generating random induced charges, 0 otherwise
|
||||
double ave_charge; // average random charge
|
||||
int mat_dim; // matrix dimension = total number of induced charges
|
||||
int mr; // number of vectors used to span the Krylov space
|
||||
int iterations; // actual number of iterations
|
||||
int itr_max; // maximum number of outer iterations
|
||||
int randomized; // 1 if generating random induced charges, 0 otherwise
|
||||
double ave_charge; // average random charge
|
||||
int seed_charge;
|
||||
|
||||
double *c, *g, *h, *r, *s, *v, *y; // vectors used by the solver
|
||||
double* rhs; // right-hand side vector of the equation Ax = b
|
||||
double tol_abs, tol_rel; // tolerance for convergence
|
||||
double normb; // norm of the rhs vector b
|
||||
double rho; // norm of (b - Ax)
|
||||
int first; // 1 if first time invoked (initializing induced charges with zero)
|
||||
double *c, *g, *h, *r, *s, *v, *y; // vectors used by the solver
|
||||
double *rhs; // right-hand side vector of the equation Ax = b
|
||||
double tol_abs, tol_rel; // tolerance for convergence
|
||||
double normb; // norm of the rhs vector b
|
||||
double rho; // norm of (b - Ax)
|
||||
int first; // 1 if first time invoked (initializing induced charges with zero)
|
||||
|
||||
void gmres_solve(double*, double*); // GMRES workhorse
|
||||
void apply_operator(double*, double*, int); // compute Ax without explicitly storing A
|
||||
void update_residual(double*, double*, int); // compute (b - Ax) directly (without computing b and Ax explcitly)
|
||||
void gmres_solve(double *, double *); // GMRES workhorse
|
||||
void apply_operator(double *, double *, int); // compute Ax without explicitly storing A
|
||||
void update_residual(double *, double *,
|
||||
int); // compute (b - Ax) directly (without computing b and Ax explcitly)
|
||||
|
||||
// Givens rotations
|
||||
inline void mult_givens(double c, double s, int k, double* g) {
|
||||
double g1 = c * g[k] - s * g[k+1];
|
||||
double g2 = s * g[k] + c * g[k+1];
|
||||
g[k] = g1;
|
||||
g[k+1] = g2;
|
||||
inline void mult_givens(double c, double s, int k, double *g)
|
||||
{
|
||||
double g1 = c * g[k] - s * g[k + 1];
|
||||
double g2 = s * g[k] + c * g[k + 1];
|
||||
g[k] = g1;
|
||||
g[k + 1] = g2;
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -33,9 +33,9 @@ class FixPolarizeBEMICC : public Fix {
|
||||
virtual void setup(int);
|
||||
virtual void pre_force(int);
|
||||
virtual double compute_vector(int);
|
||||
int modify_param(int, char**);
|
||||
int pack_forward_comm(int, int*, double*, int, int*);
|
||||
void unpack_forward_comm(int, int, double*);
|
||||
int modify_param(int, char **);
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
|
||||
virtual void compute_induced_charges();
|
||||
void set_dielectric_params(double, double, double, double, int, double);
|
||||
@ -43,27 +43,26 @@ class FixPolarizeBEMICC : public Fix {
|
||||
class AtomVecDielectric *avec;
|
||||
|
||||
protected:
|
||||
int nevery; // to be invoked every time steps
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
int torqueflag,extraflag;
|
||||
int nevery; // to be invoked every time steps
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
int torqueflag, extraflag;
|
||||
|
||||
void force_clear();
|
||||
|
||||
private:
|
||||
int iterations; // actual number of iterations
|
||||
int itr_max; // maximum number of outer iterations
|
||||
double tol_abs, tol_rel; // tolerance for convergence
|
||||
double rho; // current error
|
||||
double omega; // iterative weight
|
||||
int randomized; // 1 if generating random induced charges, 0 otherwise
|
||||
double ave_charge; // average random charge
|
||||
int iterations; // actual number of iterations
|
||||
int itr_max; // maximum number of outer iterations
|
||||
double tol_abs, tol_rel; // tolerance for convergence
|
||||
double rho; // current error
|
||||
double omega; // iterative weight
|
||||
int randomized; // 1 if generating random induced charges, 0 otherwise
|
||||
double ave_charge; // average random charge
|
||||
int seed_charge;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -30,14 +30,14 @@ class FixPolarizeFunctional : public Fix {
|
||||
~FixPolarizeFunctional();
|
||||
int setmask();
|
||||
void init();
|
||||
void init_list(int,class NeighList *);
|
||||
void init_list(int, class NeighList *);
|
||||
void setup(int);
|
||||
void setup_pre_force(int vflag);
|
||||
void pre_force(int);
|
||||
int pack_exchange(int, double *);
|
||||
int unpack_exchange(int, double *);
|
||||
int pack_forward_comm(int, int*, double*, int, int*);
|
||||
void unpack_forward_comm(int, int, double*);
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
|
||||
int modify_param(int, char **);
|
||||
double memory_usage();
|
||||
@ -49,7 +49,7 @@ class FixPolarizeFunctional : public Fix {
|
||||
|
||||
protected:
|
||||
int nmax;
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
class NeighList *list;
|
||||
|
||||
void set_dielectric_params(double, double, double, double, int, double);
|
||||
@ -59,25 +59,28 @@ class FixPolarizeFunctional : public Fix {
|
||||
double **inverse_matrix;
|
||||
double **G1ww, **ndotGww, **G2ww, **G3ww, **Rww;
|
||||
|
||||
int* tag2mat; // tag2mat[atom->tag[i]] = the index of the atoms in the induced charge arrays from atom tags
|
||||
int* mat2tag; // mat2tag[idx] = the atom tag of the induced charge idx in the induced charge arrays
|
||||
int* induced_charge_idx; // return the index of the atoms in the induced charge arrays
|
||||
int num_induced_charges; // total number of induced charges
|
||||
double* induced_charges; // values of induced charges
|
||||
int* tag2mat_ions; // tag2mat_ions[atom->tag[i]] returns the index of the atoms in the ion arrays from atom tags
|
||||
int* mat2tag_ions; // mat2tag_ions[idx] returns the atom tag of the ion idx in the ion arrays
|
||||
int* ion_idx; // return the index of the atoms in the ion arrays
|
||||
int num_ions; // total number of ions
|
||||
double* rhs1;
|
||||
double* rhs2;
|
||||
double** buffer1;
|
||||
double** buffer2;
|
||||
int *
|
||||
tag2mat; // tag2mat[atom->tag[i]] = the index of the atoms in the induced charge arrays from atom tags
|
||||
int *
|
||||
mat2tag; // mat2tag[idx] = the atom tag of the induced charge idx in the induced charge arrays
|
||||
int *induced_charge_idx; // return the index of the atoms in the induced charge arrays
|
||||
int num_induced_charges; // total number of induced charges
|
||||
double *induced_charges; // values of induced charges
|
||||
int *
|
||||
tag2mat_ions; // tag2mat_ions[atom->tag[i]] returns the index of the atoms in the ion arrays from atom tags
|
||||
int *mat2tag_ions; // mat2tag_ions[idx] returns the atom tag of the ion idx in the ion arrays
|
||||
int *ion_idx; // return the index of the atoms in the ion arrays
|
||||
int num_ions; // total number of ions
|
||||
double *rhs1;
|
||||
double *rhs2;
|
||||
double **buffer1;
|
||||
double **buffer2;
|
||||
|
||||
int allocated;
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int torqueflag,extraflag;
|
||||
int kspaceflag; // 1 if kspace is used for the induced charge computation
|
||||
double **efield_pair; // electrical field at position of atom i due to pair contribution
|
||||
double **efield_kspace; // electrical field at position of atom i due to kspace contribution
|
||||
int torqueflag, extraflag;
|
||||
double g_ewald;
|
||||
int includingG3ww;
|
||||
|
||||
@ -102,8 +105,8 @@ class FixPolarizeFunctional : public Fix {
|
||||
double tolerance;
|
||||
|
||||
void calculate_matrix_multiply_vector(double **, double *, double *, int);
|
||||
double inner_product(double*, double*, int);
|
||||
void cg_solver(double**, double*, double*, int);
|
||||
double inner_product(double *, double *, int);
|
||||
void cg_solver(double **, double *, double *, int);
|
||||
|
||||
inline double greens_real(double);
|
||||
inline double grad_greens_real_factor(double);
|
||||
@ -114,7 +117,7 @@ class FixPolarizeFunctional : public Fix {
|
||||
inline double calculate_ndotgreens_ewald_self_vertex(double, double);
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -32,14 +32,14 @@ class MSMDielectric : public MSM {
|
||||
virtual void compute(int, int);
|
||||
void fieldforce();
|
||||
|
||||
double** efield;
|
||||
double* phi;
|
||||
double **efield;
|
||||
double *phi;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -69,8 +68,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
|
||||
}
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag,vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
|
||||
@ -32,14 +32,14 @@ class PairCoulCutDielectric : public PairCoulCut {
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void init_style();
|
||||
|
||||
double** efield;
|
||||
double **efield;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -24,11 +23,11 @@
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
@ -36,13 +35,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_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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -63,31 +62,30 @@ PairCoulLongDielectric::~PairCoulLongDielectric()
|
||||
|
||||
void PairCoulLongDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itable,itype,jtype;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul;
|
||||
double fpair_i,fpair_j;
|
||||
double fraction,table;
|
||||
double r,r2inv,forcecoul,factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc,prefactorE,efield_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, inum, jnum, itable, itype, jtype;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul;
|
||||
double fpair_i, fpair_j;
|
||||
double fraction, table;
|
||||
double r, r2inv, forcecoul, factor_coul;
|
||||
double grij, expm2, prefactor, t, erfc, prefactorE, efield_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
double rsq;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
}
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
double** norm = atom->mu;
|
||||
double* curvature = atom->curvature;
|
||||
double* area = atom->area;
|
||||
double **norm = atom->mu;
|
||||
double *curvature = atom->curvature;
|
||||
double *area = atom->area;
|
||||
double *eps = atom->epsilon;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
@ -116,10 +114,10 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -132,24 +130,24 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
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;
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
expm2 = exp(-grij * grij);
|
||||
t = 1.0 / (1.0 + EWALD_P * grij);
|
||||
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp * q[j] / r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = qqrd2e * scale[itype][jtype] * q[j]/r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = qqrd2e * scale[itype][jtype] * q[j] / r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
@ -157,48 +155,48 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp * q[j] * table;
|
||||
efield_i = scale[itype][jtype] * q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = scale[itype][jtype] * qtmp * q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = scale[itype][jtype] * q[j] * table;
|
||||
efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
}
|
||||
}
|
||||
|
||||
fpair_i = etmp * forcecoul * r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
fpair_j = eps[j] * forcecoul * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(etmp+eps[j])*erfc;
|
||||
ecoul = prefactor * (etmp + eps[j]) * erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = scale[itype][jtype] * qtmp*q[j]*(etmp+eps[j]) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
ecoul = scale[itype][jtype] * qtmp * q[j] * (etmp + eps[j]) * table;
|
||||
}
|
||||
ecoul *= 0.5;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,0.0,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, 0.0, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -213,9 +211,9 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
|
||||
void PairCoulLongDielectric::init_style()
|
||||
{
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair coul/long/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair coul/long/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
|
||||
@ -223,12 +221,10 @@ void PairCoulLongDielectric::init_style()
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
if (force->kspace == NULL) error->all(FLERR, "Pair style requires a KSpace style");
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
// setup force tables
|
||||
|
||||
if (ncoultablebits) init_tables(cut_coul,NULL);
|
||||
if (ncoultablebits) init_tables(cut_coul, NULL);
|
||||
}
|
||||
|
||||
|
||||
@ -34,12 +34,11 @@ class PairCoulLongDielectric : public PairCoulLong {
|
||||
double **efield;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -23,11 +22,11 @@
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -38,8 +37,7 @@ using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulCutDielectric::PairLJCutCoulCutDielectric(LAMMPS *lmp) :
|
||||
PairLJCutCoulCut(lmp)
|
||||
PairLJCutCoulCutDielectric::PairLJCutCoulCutDielectric(LAMMPS *lmp) : PairLJCutCoulCut(lmp)
|
||||
{
|
||||
efield = nullptr;
|
||||
epot = nullptr;
|
||||
@ -58,32 +56,31 @@ PairLJCutCoulCutDielectric::~PairLJCutCoulCutDielectric()
|
||||
|
||||
void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fpair_i,fpair_j;
|
||||
double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj,efield_i,epot_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
|
||||
double fpair_i, fpair_j;
|
||||
double rsq, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj, efield_i, epot_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
double *q_real = atom->q_unscaled;
|
||||
double* eps = atom->epsilon;
|
||||
double** norm = atom->mu;
|
||||
double* curvature = atom->curvature;
|
||||
double* area = atom->area;
|
||||
double *eps = atom->epsilon;
|
||||
double **norm = atom->mu;
|
||||
double *curvature = atom->curvature;
|
||||
double *area = atom->area;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
@ -112,10 +109,10 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -131,55 +128,58 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
|
||||
if (rsq < cut_coulsq[itype][jtype] && rsq > EPSILON) {
|
||||
efield_i = q[j]*sqrt(r2inv);
|
||||
forcecoul = qqrd2e * qtmp*efield_i;
|
||||
efield_i = q[j] * sqrt(r2inv);
|
||||
forcecoul = qqrd2e * qtmp * efield_i;
|
||||
epot_i = efield_i;
|
||||
} else epot_i = efield_i = forcecoul = 0.0;
|
||||
} else
|
||||
epot_i = efield_i = forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fpair_i = (factor_coul*etmp*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
fpair_i = (factor_coul * etmp * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (factor_coul*etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (factor_coul * etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
fpair_j = (factor_coul*eps[j]*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
fpair_j = (factor_coul * eps[j] * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
ecoul = factor_coul * qqrd2e * qtmp*q[j]*(etmp+eps[j])*sqrt(r2inv);
|
||||
} else ecoul = 0.0;
|
||||
ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) * sqrt(r2inv);
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
ecoul *= 0.5;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,evdwl,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -194,47 +194,50 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
|
||||
void PairLJCutCoulCutDielectric::init_style()
|
||||
{
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair lj/cut/coul/cut/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair lj/cut/coul/cut/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj, double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,forcecoul,forcelj,phicoul,ei,ej,philj;
|
||||
double* eps = atom->epsilon;
|
||||
double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj;
|
||||
double *eps = atom->epsilon;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
if (rsq < cut_coulsq[itype][jtype])
|
||||
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv)*eps[i];
|
||||
else forcecoul = 0.0;
|
||||
forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv) * eps[i];
|
||||
else
|
||||
forcecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
fforce = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (eps[i] == 1) ei = 0;
|
||||
else ei = eps[i];
|
||||
if (eps[j] == 1) ej = 0;
|
||||
else ej = eps[j];
|
||||
if (eps[i] == 1)
|
||||
ei = 0;
|
||||
else
|
||||
ei = eps[i];
|
||||
if (eps[j] == 1)
|
||||
ej = 0;
|
||||
else
|
||||
ej = eps[j];
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
|
||||
phicoul *= 0.5*(ei+ej);
|
||||
eng += factor_coul*phicoul;
|
||||
phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * sqrt(r2inv);
|
||||
phicoul *= 0.5 * (ei + ej);
|
||||
eng += factor_coul * phicoul;
|
||||
}
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
eng += factor_lj * philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
|
||||
@ -32,15 +32,15 @@ class PairLJCutCoulCutDielectric : public PairLJCutCoulCut {
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void init_style();
|
||||
|
||||
double** efield;
|
||||
double* epot;
|
||||
double **efield;
|
||||
double *epot;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -16,21 +15,21 @@
|
||||
Contributing author: Trung Nguyen (Northwestern)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "pair_lj_cut_coul_debye_dielectric.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "atom_vec_dielectric.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
@ -58,33 +57,31 @@ PairLJCutCoulDebyeDielectric::~PairLJCutCoulDebyeDielectric()
|
||||
|
||||
void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fpair_i,fpair_j;
|
||||
double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj,efield_i,epot_i;
|
||||
double r,rinv,screening;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
|
||||
double fpair_i, fpair_j;
|
||||
double rsq, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj, efield_i, epot_i;
|
||||
double r, rinv, screening;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
double *q_real = avec->q_real;
|
||||
double* eps = avec->epsilon;
|
||||
double** norm = avec->mu;
|
||||
double* curvature = avec->curvature;
|
||||
double* area = avec->area;
|
||||
double *eps = avec->epsilon;
|
||||
double **norm = avec->mu;
|
||||
double *curvature = avec->curvature;
|
||||
double *area = avec->area;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
@ -113,10 +110,10 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag)
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -132,58 +129,61 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
|
||||
if (rsq < cut_coulsq[itype][jtype] && rsq > EPSILON) {
|
||||
r = sqrt(rsq);
|
||||
rinv = 1.0/r;
|
||||
screening = exp(-kappa*r);
|
||||
rinv = 1.0 / r;
|
||||
screening = exp(-kappa * r);
|
||||
efield_i = qqrd2e * q[j] * screening * (kappa + rinv);
|
||||
forcecoul = qtmp*efield_i;
|
||||
forcecoul = qtmp * efield_i;
|
||||
epot_i = efield_i;
|
||||
} else efield_i = forcecoul = 0.0;
|
||||
} else
|
||||
efield_i = forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fpair_i = (factor_coul*etmp*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
fpair_i = (factor_coul * etmp * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (factor_coul*etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (factor_coul * etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
fpair_j = (factor_coul*eps[j]*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
fpair_j = (factor_coul * eps[j] * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
ecoul = factor_coul * qqrd2e * qtmp*q[j]*(etmp+eps[j]) * rinv * screening;
|
||||
} else ecoul = 0.0;
|
||||
ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) * rinv * screening;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
ecoul *= 0.5;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,evdwl,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -198,52 +198,54 @@ void PairLJCutCoulDebyeDielectric::compute(int eflag, int vflag)
|
||||
void PairLJCutCoulDebyeDielectric::init_style()
|
||||
{
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair lj/cut/coul/debye/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair lj/cut/coul/debye/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
double PairLJCutCoulDebyeDielectric::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj, double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,forcecoul,forcelj,phicoul,ei,ej,philj;
|
||||
double r,rinv,screening;
|
||||
double* eps = avec->epsilon;
|
||||
double r2inv, r6inv, forcecoul, forcelj, phicoul, ei, ej, philj;
|
||||
double r, rinv, screening;
|
||||
double *eps = avec->epsilon;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
r = sqrt(rsq);
|
||||
rinv = 1.0/r;
|
||||
screening = exp(-kappa*r);
|
||||
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
|
||||
screening * (kappa + rinv) * eps[i];
|
||||
} else forcecoul = 0.0;
|
||||
rinv = 1.0 / r;
|
||||
screening = exp(-kappa * r);
|
||||
forcecoul = force->qqrd2e * atom->q[i] * atom->q[j] * screening * (kappa + rinv) * eps[i];
|
||||
} else
|
||||
forcecoul = 0.0;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
fforce = (factor_coul * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (eps[i] == 1) ei = 0;
|
||||
else ei = eps[i];
|
||||
if (eps[j] == 1) ej = 0;
|
||||
else ej = eps[j];
|
||||
if (eps[i] == 1)
|
||||
ei = 0;
|
||||
else
|
||||
ei = eps[i];
|
||||
if (eps[j] == 1)
|
||||
ej = 0;
|
||||
else
|
||||
ej = eps[j];
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * screening;
|
||||
phicoul *= 0.5*(ei+ej);
|
||||
eng += factor_coul*phicoul;
|
||||
phicoul = force->qqrd2e * atom->q[i] * atom->q[j] * rinv * screening;
|
||||
phicoul *= 0.5 * (ei + ej);
|
||||
eng += factor_coul * phicoul;
|
||||
}
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
eng += factor_lj * philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
|
||||
@ -32,15 +32,15 @@ class PairLJCutCoulDebyeDielectric : public PairLJCutCoulDebye {
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
void init_style();
|
||||
|
||||
double** efield;
|
||||
double* epot;
|
||||
double **efield;
|
||||
double *epot;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -24,31 +23,30 @@
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
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_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 EPSILON 1e-6
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLongDielectric::PairLJCutCoulLongDielectric(LAMMPS *lmp) :
|
||||
PairLJCutCoulLong(lmp)
|
||||
PairLJCutCoulLongDielectric::PairLJCutCoulLongDielectric(LAMMPS *lmp) : PairLJCutCoulLong(lmp)
|
||||
{
|
||||
respa_enable = 0;
|
||||
cut_respa = nullptr;
|
||||
@ -69,34 +67,33 @@ PairLJCutCoulLongDielectric::~PairLJCutCoulLongDielectric()
|
||||
|
||||
void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,ii,j,jj,inum,jnum,itype,jtype,itable;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
|
||||
double fpair_i,fpair_j;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc,prefactorE,efield_i,epot_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, ii, j, jj, inum, jnum, itype, jtype, itable;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul;
|
||||
double fpair_i, fpair_j;
|
||||
double fraction, table;
|
||||
double r, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj;
|
||||
double grij, expm2, prefactor, t, erfc, prefactorE, efield_i, epot_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
double rsq;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
double *eps = atom->epsilon;
|
||||
double** norm = atom->mu;
|
||||
double* curvature = atom->curvature;
|
||||
double* area = atom->area;
|
||||
double **norm = atom->mu;
|
||||
double *curvature = atom->curvature;
|
||||
double *area = atom->area;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
@ -126,10 +123,10 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
|
||||
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -145,27 +142,27 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (rsq < cut_coulsq && rsq > EPSILON) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
|
||||
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;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
expm2 = exp(-grij * grij);
|
||||
t = 1.0 / (1.0 + EWALD_P * grij);
|
||||
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
|
||||
prefactor = qqrd2e * qtmp * q[j] / r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j]/r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] / r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
epot_i = efield_i;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
@ -173,66 +170,69 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
efield_i = q[j] * table/qqrd2e;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = qtmp * q[j] * table;
|
||||
efield_i = q[j] * table / qqrd2e;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = qtmp * q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j] * table/qqrd2e;
|
||||
efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] * table / qqrd2e;
|
||||
efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
}
|
||||
epot_i = efield_i;
|
||||
}
|
||||
} else epot_i = efield_i = forcecoul = 0.0;
|
||||
} else
|
||||
epot_i = efield_i = forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fpair_i = (forcecoul*etmp + factor_lj*forcelj) * r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
fpair_i = (forcecoul * etmp + factor_lj * forcelj) * r2inv;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
|
||||
fpair_j = (forcecoul*eps[j] + factor_lj*forcelj) * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
fpair_j = (forcecoul * eps[j] + factor_lj * forcelj) * r2inv;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(etmp+eps[j])*erfc;
|
||||
ecoul = prefactor * (etmp + eps[j]) * erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j]*(etmp+eps[j]) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
ecoul = qtmp * q[j] * (etmp + eps[j]) * table;
|
||||
}
|
||||
ecoul *= 0.5;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,evdwl,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -247,9 +247,9 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
|
||||
void PairLJCutCoulLongDielectric::init_style()
|
||||
{
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair lj/cut/coul/long/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair lj/cut/coul/long/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
|
||||
@ -257,82 +257,84 @@ void PairLJCutCoulLongDielectric::init_style()
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
if (force->kspace == NULL) error->all(FLERR, "Pair style requires a KSpace style");
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
// setup force tables
|
||||
|
||||
if (ncoultablebits) init_tables(cut_coul,cut_respa);
|
||||
if (ncoultablebits) init_tables(cut_coul, cut_respa);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj, double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,grij,expm2,t,erfc,ei,ej,prefactor;
|
||||
double fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
double r2inv, r6inv, r, grij, expm2, t, erfc, ei, ej, prefactor;
|
||||
double fraction, table, forcecoul, forcelj, phicoul, philj;
|
||||
int itable;
|
||||
double *eps = atom->epsilon;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
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;
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
expm2 = exp(-grij * grij);
|
||||
t = 1.0 / (1.0 + EWALD_P * grij);
|
||||
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
|
||||
prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup_single;
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = atom->q[i] * atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = atom->q[i]*atom->q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = atom->q[i] * atom->q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
} else
|
||||
forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fforce = (forcecoul*eps[i] + factor_lj*forcelj) * r2inv;
|
||||
fforce = (forcecoul * eps[i] + factor_lj * forcelj) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (eps[i] == 1) ei = 0;
|
||||
else ei = eps[i];
|
||||
if (eps[j] == 1) ej = 0;
|
||||
else ej = eps[j];
|
||||
if (eps[i] == 1)
|
||||
ei = 0;
|
||||
else
|
||||
ei = eps[i];
|
||||
if (eps[j] == 1)
|
||||
ej = 0;
|
||||
else
|
||||
ej = eps[j];
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
phicoul = prefactor*(ei+ej)*erfc;
|
||||
phicoul = prefactor * (ei + ej) * erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
phicoul = atom->q[i]*atom->q[j]*(ei+ej) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
phicoul = atom->q[i] * atom->q[j] * (ei + ej) * table;
|
||||
}
|
||||
phicoul *= 0.5;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
|
||||
eng += phicoul;
|
||||
}
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
eng += factor_lj * philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
|
||||
@ -33,15 +33,15 @@ class PairLJCutCoulLongDielectric : public PairLJCutCoulLong {
|
||||
virtual void init_style();
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
|
||||
double** efield;
|
||||
double* epot;
|
||||
double **efield;
|
||||
double *epot;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -23,13 +22,13 @@
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "integrate.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "kspace.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -40,8 +39,7 @@ using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulMSMDielectric::PairLJCutCoulMSMDielectric(LAMMPS *lmp) :
|
||||
PairLJCutCoulLong(lmp)
|
||||
PairLJCutCoulMSMDielectric::PairLJCutCoulMSMDielectric(LAMMPS *lmp) : PairLJCutCoulLong(lmp)
|
||||
{
|
||||
ewaldflag = pppmflag = 0;
|
||||
msmflag = 1;
|
||||
@ -65,52 +63,50 @@ PairLJCutCoulMSMDielectric::~PairLJCutCoulMSMDielectric()
|
||||
|
||||
void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,ii,j,jj,inum,jnum,itype,jtype,itable;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul;
|
||||
double fpair_i,fpair_j;
|
||||
double fraction,table;
|
||||
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double egamma,fgamma,prefactor,prefactorE,efield_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, ii, j, jj, inum, jnum, itype, jtype, itable;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair, fcoul;
|
||||
double fpair_i, fpair_j;
|
||||
double fraction, table;
|
||||
double r, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj;
|
||||
double egamma, fgamma, prefactor, prefactorE, efield_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
double rsq;
|
||||
int eflag_old = eflag;
|
||||
|
||||
if (force->kspace->scalar_pressure_flag && vflag) {
|
||||
if (vflag > 2)
|
||||
error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' "
|
||||
"to obtain per-atom virial with kspace_style MSM");
|
||||
error->all(FLERR,
|
||||
"Must use 'kspace_modify pressure/scalar no' "
|
||||
"to obtain per-atom virial with kspace_style MSM");
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
if (ftmp) memory->destroy(ftmp);
|
||||
nmax = atom->nmax;
|
||||
memory->create(ftmp,nmax,3,"pair:ftmp");
|
||||
memory->create(ftmp, nmax, 3, "pair:ftmp");
|
||||
}
|
||||
memset(&ftmp[0][0],0,nmax*3*sizeof(double));
|
||||
memset(&ftmp[0][0], 0, nmax * 3 * sizeof(double));
|
||||
|
||||
// must switch on global energy computation if not already on
|
||||
|
||||
if (eflag == 0 || eflag == 2) {
|
||||
eflag++;
|
||||
}
|
||||
if (eflag == 0 || eflag == 2) { eflag++; }
|
||||
}
|
||||
|
||||
if (!efield || atom->nmax > nmax) {
|
||||
if (efield) memory->destroy(efield);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
}
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
double *eps = atom->epsilon;
|
||||
double** norm = atom->mu;
|
||||
double* curvature = atom->curvature;
|
||||
double* area = atom->area;
|
||||
double **norm = atom->mu;
|
||||
double *curvature = atom->curvature;
|
||||
double *area = atom->area;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
@ -139,10 +135,10 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -156,24 +152,24 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
|
||||
if (rsq < cut_coulsq && rsq > EPSILON) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
prefactor = qqrd2e * qtmp * q[j] / r;
|
||||
egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul);
|
||||
fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j]/r;
|
||||
prefactorE = q[j] / r;
|
||||
efield_i = prefactorE * fgamma;
|
||||
if (factor_coul < 1.0) efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
@ -181,95 +177,98 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
efield_i = q[j] * table/qqrd2e;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = qtmp * q[j] * table;
|
||||
efield_i = q[j] * table / qqrd2e;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = qtmp * q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j] * table/qqrd2e;
|
||||
efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] * table / qqrd2e;
|
||||
efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
} else
|
||||
forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
if (!(force->kspace->scalar_pressure_flag && vflag)) {
|
||||
|
||||
fpair_i = (forcecoul*etmp + factor_lj*forcelj) * r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
fpair_i = (forcecoul * etmp + factor_lj * forcelj) * r2inv;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
fpair_j = (forcecoul*eps[j] + factor_lj*forcelj) * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
fpair_j = (forcecoul * eps[j] + factor_lj * forcelj) * r2inv;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
} else {
|
||||
|
||||
// separate LJ and Coulombic forces
|
||||
|
||||
fpair = (factor_lj*forcelj) * r2inv;
|
||||
fpair = (factor_lj * forcelj) * r2inv;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
f[i][0] += delx * fpair;
|
||||
f[i][1] += dely * fpair;
|
||||
f[i][2] += delz * fpair;
|
||||
if (newton_pair) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
|
||||
fpair_i = (forcecoul*etmp) * r2inv;
|
||||
ftmp[i][0] += delx*fpair_i;
|
||||
ftmp[i][1] += dely*fpair_i;
|
||||
ftmp[i][2] += delz*fpair_i;
|
||||
fpair_i = (forcecoul * etmp) * r2inv;
|
||||
ftmp[i][0] += delx * fpair_i;
|
||||
ftmp[i][1] += dely * fpair_i;
|
||||
ftmp[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
fpair_j = (forcecoul*eps[j]) * r2inv;
|
||||
ftmp[j][0] -= delx*fpair_j;
|
||||
ftmp[j][1] -= dely*fpair_j;
|
||||
ftmp[j][2] -= delz*fpair_j;
|
||||
fpair_j = (forcecoul * eps[j]) * r2inv;
|
||||
ftmp[j][0] -= delx * fpair_j;
|
||||
ftmp[j][1] -= dely * fpair_j;
|
||||
ftmp[j][2] -= delz * fpair_j;
|
||||
}
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(etmp+eps[j])*egamma;
|
||||
ecoul = prefactor * (etmp + eps[j]) * egamma;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j]*(etmp+eps[j]) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
ecoul = qtmp * q[j] * (etmp + eps[j]) * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
|
||||
if (eflag_old && rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,evdwl,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -277,7 +276,7 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
if (force->kspace->scalar_pressure_flag && vflag) {
|
||||
for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0;
|
||||
for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul / 3.0;
|
||||
for (int i = 0; i < nmax; i++) {
|
||||
f[i][0] += ftmp[i][0];
|
||||
f[i][1] += ftmp[i][1];
|
||||
@ -288,63 +287,62 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype,
|
||||
double rsq,
|
||||
double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj, double &fforce)
|
||||
{
|
||||
double r2inv,r6inv,r,egamma,fgamma,prefactor;
|
||||
double fraction,table,forcecoul,forcelj,phicoul,philj;
|
||||
double r2inv, r6inv, r, egamma, fgamma, prefactor;
|
||||
double fraction, table, forcecoul, forcelj, phicoul, philj;
|
||||
int itable;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
|
||||
fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
|
||||
prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
|
||||
egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul);
|
||||
fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul);
|
||||
forcecoul = prefactor * fgamma;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup_single;
|
||||
rsq_lookup_single.f = rsq;
|
||||
itable = rsq_lookup_single.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = atom->q[i] * atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = atom->q[i]*atom->q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = atom->q[i] * atom->q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
}
|
||||
}
|
||||
} else forcecoul = 0.0;
|
||||
} else
|
||||
forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fforce = (forcecoul + factor_lj*forcelj) * r2inv;
|
||||
fforce = (forcecoul + factor_lj * forcelj) * r2inv;
|
||||
|
||||
double eng = 0.0;
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
phicoul = prefactor*egamma;
|
||||
phicoul = prefactor * egamma;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
phicoul = atom->q[i]*atom->q[j] * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
phicoul = atom->q[i] * atom->q[j] * table;
|
||||
}
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
|
||||
eng += phicoul;
|
||||
}
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
eng += factor_lj*philj;
|
||||
philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
eng += factor_lj * philj;
|
||||
}
|
||||
|
||||
return eng;
|
||||
@ -357,9 +355,9 @@ double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype,
|
||||
void PairLJCutCoulMSMDielectric::init_style()
|
||||
{
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair lj/cut/coul/msm/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair lj/cut/coul/msm/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
|
||||
@ -367,13 +365,12 @@ void PairLJCutCoulMSMDielectric::init_style()
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
if (force->kspace == NULL) error->all(FLERR, "Pair style requires a KSpace style");
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
// setup force tables
|
||||
|
||||
if (ncoultablebits) init_tables(cut_coul,cut_respa);
|
||||
if (ncoultablebits) init_tables(cut_coul, cut_respa);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -381,9 +378,9 @@ void PairLJCutCoulMSMDielectric::init_style()
|
||||
void *PairLJCutCoulMSMDielectric::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 0;
|
||||
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
|
||||
if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul;
|
||||
dim = 2;
|
||||
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str,"sigma") == 0) return (void *) sigma;
|
||||
if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
|
||||
if (strcmp(str, "sigma") == 0) return (void *) sigma;
|
||||
return NULL;
|
||||
}
|
||||
|
||||
@ -36,12 +36,12 @@ class PairLJCutCoulMSMDielectric : public PairLJCutCoulLong {
|
||||
double **efield;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
double **ftmp;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/ Sandia National Laboratories
|
||||
@ -16,42 +15,42 @@
|
||||
Contributing author: Trung Nguyen (Northwestern)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
#include "math_vector.h"
|
||||
#include "pair_lj_long_coul_long_dielectric.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "atom_vec_dielectric.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "integrate.h"
|
||||
#include "kspace.h"
|
||||
#include "math_const.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "neighbor.h"
|
||||
#include "respa.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using namespace MathExtra;
|
||||
|
||||
#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_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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp)
|
||||
: PairLJLongCoulLong(lmp)
|
||||
PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp) : PairLJLongCoulLong(lmp)
|
||||
{
|
||||
respa_enable = 0;
|
||||
cut_respa = NULL;
|
||||
@ -79,9 +78,9 @@ void PairLJLongCoulLongDielectric::init_style()
|
||||
PairLJLongCoulLongDielectric::init_style();
|
||||
|
||||
avec = (AtomVecDielectric *) atom->style_match("dielectric");
|
||||
if (!avec) error->all(FLERR,"Pair lj/long/coul/long/dielectric requires atom style dielectric");
|
||||
if (!avec) error->all(FLERR, "Pair lj/long/coul/long/dielectric requires atom style dielectric");
|
||||
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
int irequest = neighbor->request(this, instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
}
|
||||
@ -92,26 +91,25 @@ void PairLJLongCoulLongDielectric::init_style()
|
||||
|
||||
void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
|
||||
{
|
||||
double evdwl,ecoul,fpair;
|
||||
double evdwl, ecoul, fpair;
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
double **x = atom->x, *x0 = x[0];
|
||||
double **f = atom->f, *f0 = f[0], *fi = f0;
|
||||
double *q = atom->q;
|
||||
double *eps = avec->epsilon;
|
||||
double** norm = avec->mu;
|
||||
double* curvature = avec->curvature;
|
||||
double* area = avec->area;
|
||||
double **norm = avec->mu;
|
||||
double *curvature = avec->curvature;
|
||||
double *area = avec->area;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
@ -119,46 +117,51 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
int i,ii,j,jj,inum,jnum,itype,jtype,itable;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz;
|
||||
int order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6);
|
||||
int i, ii, j, jj, inum, jnum, itype, jtype, itable;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz;
|
||||
int order1 = ewald_order & (1 << 1), order6 = ewald_order & (1 << 6);
|
||||
int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni;
|
||||
double qi = 0.0, qri = 0.0;
|
||||
double fpair_i,fpair_j;
|
||||
double fraction,table;
|
||||
double fpair_i, fpair_j;
|
||||
double fraction, table;
|
||||
double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti;
|
||||
double grij,expm2,prefactor,t,erfc,prefactorE,efield_i,epot_i;
|
||||
double r,rsq,r2inv,force_coul,force_lj,factor_coul,factor_lj;
|
||||
double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
|
||||
vector xi, d;
|
||||
double grij, expm2, prefactor, t, erfc, prefactorE, efield_i, epot_i;
|
||||
double r, rsq, r2inv, force_coul, force_lj, factor_coul, factor_lj;
|
||||
double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2;
|
||||
double xi[3], d[3];
|
||||
|
||||
ineighn = (ineigh = list->ilist)+list->inum;
|
||||
ineighn = (ineigh = list->ilist) + list->inum;
|
||||
|
||||
for (; ineigh<ineighn; ++ineigh) { // loop over my atoms
|
||||
i = *ineigh; fi = f0+3*i;
|
||||
for (; ineigh < ineighn; ++ineigh) { // loop over my atoms
|
||||
i = *ineigh;
|
||||
fi = f0 + 3 * i;
|
||||
qtmp = q[i];
|
||||
|
||||
if (order1) qri = (qi = q[i])*qqrd2e; // initialize constants
|
||||
if (order1) qri = (qi = q[i]) * qqrd2e; // initialize constants
|
||||
offseti = offset[typei = type[i]];
|
||||
lj1i = lj1[typei]; lj2i = lj2[typei]; lj3i = lj3[typei]; lj4i = lj4[typei];
|
||||
cutsqi = cutsq[typei]; cut_ljsqi = cut_ljsq[typei];
|
||||
memcpy(xi, x0+(i+(i<<1)), sizeof(vector));
|
||||
jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];
|
||||
lj1i = lj1[typei];
|
||||
lj2i = lj2[typei];
|
||||
lj3i = lj3[typei];
|
||||
lj4i = lj4[typei];
|
||||
cutsqi = cutsq[typei];
|
||||
cut_ljsqi = cut_ljsq[typei];
|
||||
memcpy(xi, x0 + (i + (i << 1)), 3 * sizeof(double));
|
||||
jneighn = (jneigh = list->firstneigh[i]) + list->numneigh[i];
|
||||
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i][0];
|
||||
efield[i][1] = sf*norm[i][1];
|
||||
efield[i][2] = sf*norm[i][2];
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i][0];
|
||||
efield[i][1] = sf * norm[i][1];
|
||||
efield[i][2] = sf * norm[i][2];
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
|
||||
epot[i] = 0;
|
||||
|
||||
for (; jneigh<jneighn; ++jneigh) { // loop over neighbors
|
||||
for (; jneigh < jneighn; ++jneigh) { // loop over neighbors
|
||||
j = *jneigh;
|
||||
ni = sbmask(j);
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
@ -168,28 +171,28 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq >= cutsq[typei][typej]) continue;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
r = sqrt(rsq);
|
||||
|
||||
if (order1 && (rsq < cut_coulsq)) { // coulombic
|
||||
if (order1 && (rsq < cut_coulsq)) { // coulombic
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
|
||||
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;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
force_coul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) force_coul -= (1.0-factor_coul)*prefactor;
|
||||
expm2 = exp(-grij * grij);
|
||||
t = 1.0 / (1.0 + EWALD_P * grij);
|
||||
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
|
||||
prefactor = qqrd2e * qtmp * q[j] / r;
|
||||
force_coul = prefactor * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) force_coul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j]/r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] / r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
epot_i = efield_i;
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
@ -197,111 +200,113 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
force_coul = qtmp*q[j] * table;
|
||||
efield_i = q[j] * table/qqrd2e;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
force_coul = qtmp * q[j] * table;
|
||||
efield_i = q[j] * table / qqrd2e;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
force_coul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = qtmp * q[j] * table;
|
||||
force_coul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j] * table/qqrd2e;
|
||||
efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] * table / qqrd2e;
|
||||
efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
}
|
||||
epot_i = efield_i;
|
||||
}
|
||||
}
|
||||
else epot_i = efield_i = force_coul = ecoul = 0.0;
|
||||
} else
|
||||
epot_i = efield_i = force_coul = ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsqi[typej]) { // lj
|
||||
if (rsq < cut_ljsqi[typej]) { // lj
|
||||
if (order6) { // long-range lj
|
||||
if(!ndisptablebits || rsq <= tabinnerdispsq) { // series real space
|
||||
double rn = r2inv*r2inv*r2inv;
|
||||
double x2 = g2*rsq, a2 = 1.0/x2;
|
||||
x2 = a2*exp(-x2)*lj4i[typej];
|
||||
if (!ndisptablebits || rsq <= tabinnerdispsq) { // series real space
|
||||
double rn = r2inv * r2inv * r2inv;
|
||||
double x2 = g2 * rsq, a2 = 1.0 / x2;
|
||||
x2 = a2 * exp(-x2) * lj4i[typej];
|
||||
if (ni == 0) {
|
||||
force_lj =
|
||||
(rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
|
||||
force_lj = (rn *= rn) * lj1i[typej] -
|
||||
g8 * (((6.0 * a2 + 6.0) * a2 + 3.0) * a2 + 1.0) * x2 * rsq;
|
||||
if (eflag) evdwl = rn * lj3i[typej] - g6 * ((a2 + 1.0) * a2 + 0.5) * x2;
|
||||
} else { // special case
|
||||
double f = special_lj[ni], t = rn * (1.0 - f);
|
||||
force_lj = f * (rn *= rn) * lj1i[typej] -
|
||||
g8 * (((6.0 * a2 + 6.0) * a2 + 3.0) * a2 + 1.0) * x2 * rsq + t * lj2i[typej];
|
||||
if (eflag)
|
||||
evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
|
||||
evdwl = f * rn * lj3i[typej] - g6 * ((a2 + 1.0) * a2 + 0.5) * x2 + t * lj4i[typej];
|
||||
}
|
||||
else { // special case
|
||||
double f = special_lj[ni], t = rn*(1.0-f);
|
||||
force_lj = f*(rn *= rn)*lj1i[typej]-
|
||||
g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
|
||||
if (eflag)
|
||||
evdwl = f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej];
|
||||
}
|
||||
}
|
||||
else { // table real space
|
||||
} else { // table real space
|
||||
union_int_float_t disp_t;
|
||||
disp_t.f = rsq;
|
||||
const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
|
||||
double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
|
||||
double rn = r2inv*r2inv*r2inv;
|
||||
const int disp_k = (disp_t.i & ndispmask) >> ndispshiftbits;
|
||||
double f_disp = (rsq - rdisptable[disp_k]) * drdisptable[disp_k];
|
||||
double rn = r2inv * r2inv * r2inv;
|
||||
if (ni == 0) {
|
||||
force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
|
||||
if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
|
||||
}
|
||||
else { // special case
|
||||
double f = special_lj[ni], t = rn*(1.0-f);
|
||||
force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
|
||||
if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
|
||||
force_lj = (rn *= rn) * lj1i[typej] -
|
||||
(fdisptable[disp_k] + f_disp * dfdisptable[disp_k]) * lj4i[typej];
|
||||
if (eflag)
|
||||
evdwl = rn * lj3i[typej] -
|
||||
(edisptable[disp_k] + f_disp * dedisptable[disp_k]) * lj4i[typej];
|
||||
} else { // special case
|
||||
double f = special_lj[ni], t = rn * (1.0 - f);
|
||||
force_lj = f * (rn *= rn) * lj1i[typej] -
|
||||
(fdisptable[disp_k] + f_disp * dfdisptable[disp_k]) * lj4i[typej] +
|
||||
t * lj2i[typej];
|
||||
if (eflag)
|
||||
evdwl = f * rn * lj3i[typej] -
|
||||
(edisptable[disp_k] + f_disp * dedisptable[disp_k]) * lj4i[typej] +
|
||||
t * lj4i[typej];
|
||||
}
|
||||
}
|
||||
}
|
||||
else { // cut lj
|
||||
double rn = r2inv*r2inv*r2inv;
|
||||
} else { // cut lj
|
||||
double rn = r2inv * r2inv * r2inv;
|
||||
if (ni == 0) {
|
||||
force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
|
||||
if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
|
||||
}
|
||||
else { // special case
|
||||
force_lj = rn * (rn * lj1i[typej] - lj2i[typej]);
|
||||
if (eflag) evdwl = rn * (rn * lj3i[typej] - lj4i[typej]) - offseti[typej];
|
||||
} else { // special case
|
||||
double f = special_lj[ni];
|
||||
force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
|
||||
if (eflag)
|
||||
evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
|
||||
force_lj = f * rn * (rn * lj1i[typej] - lj2i[typej]);
|
||||
if (eflag) evdwl = f * (rn * (rn * lj3i[typej] - lj4i[typej]) - offseti[typej]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
else force_lj = evdwl = 0.0;
|
||||
else
|
||||
force_lj = evdwl = 0.0;
|
||||
|
||||
fpair = (force_coul*etmp+force_lj)*r2inv;
|
||||
f[i][0] += delx*fpair_i;
|
||||
f[i][1] += dely*fpair_i;
|
||||
f[i][2] += delz*fpair_i;
|
||||
fpair = (force_coul * etmp + force_lj) * r2inv;
|
||||
f[i][0] += delx * fpair_i;
|
||||
f[i][1] += dely * fpair_i;
|
||||
f[i][2] += delz * fpair_i;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
efield[i][0] += delx*efield_i;
|
||||
efield[i][1] += dely*efield_i;
|
||||
efield[i][2] += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
efield[i][0] += delx * efield_i;
|
||||
efield[i][1] += dely * efield_i;
|
||||
efield[i][2] += delz * efield_i;
|
||||
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (newton_pair && j >= nlocal) {
|
||||
|
||||
fpair_j = (force_coul*eps[j] + factor_lj*force_lj) * r2inv;
|
||||
f[j][0] -= delx*fpair_j;
|
||||
f[j][1] -= dely*fpair_j;
|
||||
f[j][2] -= delz*fpair_j;
|
||||
fpair_j = (force_coul * eps[j] + factor_lj * force_lj) * r2inv;
|
||||
f[j][0] -= delx * fpair_j;
|
||||
f[j][1] -= dely * fpair_j;
|
||||
f[j][2] -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(etmp+eps[j])*erfc;
|
||||
ecoul = prefactor * (etmp + eps[j]) * erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j]*(etmp+eps[j]) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
ecoul = qtmp * q[j] * (etmp + eps[j]) * table;
|
||||
}
|
||||
ecoul *= 0.5;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally_full(i,evdwl,ecoul,fpair_i,delx,dely,delz);
|
||||
if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
|
||||
@ -310,57 +315,61 @@ void PairLJLongCoulLongDielectric::compute(int eflag, int vflag)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJLongCoulLongDielectric::single(int i, int j, int itype, int jtype,
|
||||
double rsq, double factor_coul, double factor_lj,
|
||||
double &fforce)
|
||||
double PairLJLongCoulLongDielectric::single(int i, int j, int itype, int jtype, double rsq,
|
||||
double factor_coul, double factor_lj, double &fforce)
|
||||
{
|
||||
double r2inv, r6inv, force_coul, force_lj, ei, ej;
|
||||
double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q;
|
||||
double g2 = g_ewald_6 * g_ewald_6, g6 = g2 * g2 * g2, g8 = g6 * g2, *q = atom->q;
|
||||
double *eps = avec->epsilon;
|
||||
|
||||
double eng = 0.0;
|
||||
if (eps[i] == 1) ei = 0;
|
||||
else ei = eps[i];
|
||||
if (eps[j] == 1) ej = 0;
|
||||
else ej = eps[j];
|
||||
if (eps[i] == 1)
|
||||
ei = 0;
|
||||
else
|
||||
ei = eps[i];
|
||||
if (eps[j] == 1)
|
||||
ej = 0;
|
||||
else
|
||||
ej = eps[j];
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
if ((ewald_order&2) && (rsq < cut_coulsq)) { // coulombic
|
||||
if (!ncoultablebits || rsq <= tabinnersq) { // series real space
|
||||
double r = sqrt(rsq), x = g_ewald*r;
|
||||
double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
|
||||
r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
|
||||
force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r;
|
||||
eng += (t-r)*(ei+ej)*0.5;
|
||||
}
|
||||
else { // table real space
|
||||
r2inv = 1.0 / rsq;
|
||||
if ((ewald_order & 2) && (rsq < cut_coulsq)) { // coulombic
|
||||
if (!ncoultablebits || rsq <= tabinnersq) { // series real space
|
||||
double r = sqrt(rsq), x = g_ewald * r;
|
||||
double s = force->qqrd2e * q[i] * q[j], t = 1.0 / (1.0 + EWALD_P * x);
|
||||
r = s * (1.0 - factor_coul) / r;
|
||||
s *= g_ewald * exp(-x * x);
|
||||
force_coul = (t *= ((((t * A5 + A4) * t + A3) * t + A2) * t + A1) * s / x) + EWALD_F * s - r;
|
||||
eng += (t - r) * (ei + ej) * 0.5;
|
||||
} else { // table real space
|
||||
union_int_float_t t;
|
||||
t.f = rsq;
|
||||
const int k = (t.i & ncoulmask) >> ncoulshiftbits;
|
||||
double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
|
||||
t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
|
||||
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
|
||||
eng += qiqj*(etable[k]+f*detable[k]-t.f)*(ei+ej)*0.5;
|
||||
double f = (rsq - rtable[k]) * drtable[k], qiqj = q[i] * q[j];
|
||||
t.f = (1.0 - factor_coul) * (ctable[k] + f * dctable[k]);
|
||||
force_coul = qiqj * (ftable[k] + f * dftable[k] - t.f);
|
||||
eng += qiqj * (etable[k] + f * detable[k] - t.f) * (ei + ej) * 0.5;
|
||||
}
|
||||
} else force_coul = 0.0;
|
||||
} else
|
||||
force_coul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) { // lennard-jones
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
if (ewald_order&64) { // long-range
|
||||
double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
|
||||
x2 = a2*exp(-x2)*lj4[itype][jtype];
|
||||
force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
|
||||
g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
|
||||
eng += factor_lj*r6inv*lj3[itype][jtype]-
|
||||
g6*((a2+1.0)*a2+0.5)*x2+t*lj4[itype][jtype];
|
||||
if (rsq < cut_ljsq[itype][jtype]) { // lennard-jones
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
if (ewald_order & 64) { // long-range
|
||||
double x2 = g2 * rsq, a2 = 1.0 / x2, t = r6inv * (1.0 - factor_lj);
|
||||
x2 = a2 * exp(-x2) * lj4[itype][jtype];
|
||||
force_lj = factor_lj * (r6inv *= r6inv) * lj1[itype][jtype] -
|
||||
g8 * (((6.0 * a2 + 6.0) * a2 + 3.0) * a2 + a2) * x2 * rsq + t * lj2[itype][jtype];
|
||||
eng += factor_lj * r6inv * lj3[itype][jtype] - g6 * ((a2 + 1.0) * a2 + 0.5) * x2 +
|
||||
t * lj4[itype][jtype];
|
||||
} else { // cut
|
||||
force_lj = factor_lj * r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
eng += factor_lj *
|
||||
(r6inv * (r6inv * lj3[itype][jtype] - lj4[itype][jtype]) - offset[itype][jtype]);
|
||||
}
|
||||
else { // cut
|
||||
force_lj = factor_lj*r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
|
||||
eng += factor_lj*(r6inv*(r6inv*lj3[itype][jtype]-
|
||||
lj4[itype][jtype])-offset[itype][jtype]);
|
||||
}
|
||||
} else force_lj = 0.0;
|
||||
} else
|
||||
force_lj = 0.0;
|
||||
|
||||
fforce = (force_coul*eps[i]+force_lj)*r2inv;
|
||||
fforce = (force_coul * eps[i] + force_lj) * r2inv;
|
||||
return eng;
|
||||
}
|
||||
|
||||
@ -32,15 +32,15 @@ class PairLJLongCoulLongDielectric : public PairLJLongCoulLong {
|
||||
void init_style();
|
||||
double single(int, int, int, int, double, double, double, double &);
|
||||
|
||||
double** efield;
|
||||
double* epot;
|
||||
double **efield;
|
||||
double *epot;
|
||||
|
||||
protected:
|
||||
class AtomVecDielectric* avec;
|
||||
class AtomVecDielectric *avec;
|
||||
int nmax;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -34,7 +34,6 @@
|
||||
#include "remap_wrap.h"
|
||||
|
||||
#include <cmath>
|
||||
//#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
@ -30,9 +30,9 @@ class PPPMDielectric : public PPPM {
|
||||
virtual ~PPPMDielectric();
|
||||
virtual void compute(int, int);
|
||||
|
||||
double** efield;
|
||||
double* phi;
|
||||
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
|
||||
double **efield;
|
||||
double *phi;
|
||||
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
|
||||
|
||||
void qsum_qsq();
|
||||
|
||||
@ -42,11 +42,10 @@ class PPPMDielectric : public PPPM {
|
||||
void fieldforce_ik();
|
||||
void fieldforce_ad();
|
||||
|
||||
class AtomVecDielectric* avec;
|
||||
|
||||
class AtomVecDielectric *avec;
|
||||
};
|
||||
|
||||
}
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@ -19,16 +19,17 @@
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include "omp_compat.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
@ -37,31 +38,26 @@ using namespace MathConst;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulCutDielectricOMP::PairLJCutCoulCutDielectricOMP(LAMMPS *lmp) :
|
||||
PairLJCutCoulCutDielectric(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
PairLJCutCoulCutDielectric(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulCutDielectricOMP::~PairLJCutCoulCutDielectricOMP()
|
||||
{
|
||||
}
|
||||
PairLJCutCoulCutDielectricOMP::~PairLJCutCoulCutDielectricOMP() {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulCutDielectricOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
if (eflag || vflag) {
|
||||
ev_setup(eflag,vflag);
|
||||
} else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
@ -69,7 +65,7 @@ void PairLJCutCoulCutDielectricOMP::compute(int eflag, int vflag)
|
||||
const int inum = list->inum;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel default(none) shared(eflag,vflag)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
@ -81,49 +77,55 @@ void PairLJCutCoulCutDielectricOMP::compute(int eflag, int vflag)
|
||||
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
|
||||
else eval<1,1,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<1, 1, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 1, 0>(ifrom, ito, thr);
|
||||
} else {
|
||||
if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
|
||||
else eval<1,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<1, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
} else {
|
||||
if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
|
||||
else eval<0,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<0, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<0, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
|
||||
thr->timer(Timer::PAIR);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData *const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype,itable;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fpair_i,fpair_j;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc,prefactorE,efield_i,epot_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, jnum, itype, jtype, itable;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
|
||||
double fpair_i, fpair_j;
|
||||
double r, rsq, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj;
|
||||
double grij, expm2, prefactor, t, erfc, prefactorE, efield_i, epot_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const double * _noalias const q = atom->q;
|
||||
const double * _noalias const eps = atom->epsilon;
|
||||
const dbl3_t * _noalias const norm = (dbl3_t *) atom->mu[0];
|
||||
const double * _noalias const curvature = atom->curvature;
|
||||
const double * _noalias const area = atom->area;
|
||||
const int * _noalias const type = atom->type;
|
||||
const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const double *_noalias const q = atom->q;
|
||||
const double *_noalias const eps = atom->epsilon;
|
||||
const dbl3_t *_noalias const norm = (dbl3_t *) atom->mu[0];
|
||||
const double *_noalias const curvature = atom->curvature;
|
||||
const double *_noalias const area = atom->area;
|
||||
const int *_noalias const type = atom->type;
|
||||
const int nlocal = atom->nlocal;
|
||||
const double * _noalias const special_coul = force->special_coul;
|
||||
const double * _noalias const special_lj = force->special_lj;
|
||||
const double *_noalias const special_coul = force->special_coul;
|
||||
const double *_noalias const special_lj = force->special_lj;
|
||||
const double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp,extmp,eytmp,eztmp;
|
||||
double fxtmp, fytmp, fztmp, extmp, eytmp, eztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
@ -142,16 +144,16 @@ void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData * const t
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp=fytmp=fztmp=0.0;
|
||||
extmp=eytmp=eztmp=0.0;
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
extmp = eytmp = eztmp = 0.0;
|
||||
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i].x;
|
||||
efield[i][1] = sf*norm[i].y;
|
||||
efield[i][2] = sf*norm[i].z;
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i].x;
|
||||
efield[i][1] = sf * norm[i].y;
|
||||
efield[i][2] = sf * norm[i].z;
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -167,56 +169,59 @@ void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData * const t
|
||||
delx = xtmp - x[j].x;
|
||||
dely = ytmp - x[j].y;
|
||||
delz = ztmp - x[j].z;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
|
||||
if (rsq < cut_coulsq[itype][jtype] && rsq > EPSILON) {
|
||||
efield_i = q[j]*sqrt(r2inv);
|
||||
forcecoul = qqrd2e * qtmp*efield_i;
|
||||
efield_i = q[j] * sqrt(r2inv);
|
||||
forcecoul = qqrd2e * qtmp * efield_i;
|
||||
epot_i = efield_i;
|
||||
} else epot_i = efield_i = forcecoul = 0.0;
|
||||
} else
|
||||
epot_i = efield_i = forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
} else forcelj = 0.0;
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fpair_i = (factor_coul*etmp*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
fpair_i = (factor_coul * etmp * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
|
||||
fxtmp += delx*fpair_i;
|
||||
fytmp += dely*fpair_i;
|
||||
fztmp += delz*fpair_i;
|
||||
fxtmp += delx * fpair_i;
|
||||
fytmp += dely * fpair_i;
|
||||
fztmp += delz * fpair_i;
|
||||
|
||||
efield_i *= (factor_coul*etmp*r2inv);
|
||||
extmp += delx*efield_i;
|
||||
eytmp += dely*efield_i;
|
||||
eztmp += delz*efield_i;
|
||||
efield_i *= (factor_coul * etmp * r2inv);
|
||||
extmp += delx * efield_i;
|
||||
eytmp += dely * efield_i;
|
||||
eztmp += delz * efield_i;
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (NEWTON_PAIR || j >= nlocal) {
|
||||
fpair_j = (factor_coul*eps[j]*forcecoul + factor_lj*forcelj) * r2inv;
|
||||
f[j].x -= delx*fpair_j;
|
||||
f[j].y -= dely*fpair_j;
|
||||
f[j].z -= delz*fpair_j;
|
||||
fpair_j = (factor_coul * eps[j] * forcecoul + factor_lj * forcelj) * r2inv;
|
||||
f[j].x -= delx * fpair_j;
|
||||
f[j].y -= dely * fpair_j;
|
||||
f[j].z -= delz * fpair_j;
|
||||
}
|
||||
|
||||
if (EFLAG) {
|
||||
if (rsq < cut_coulsq[itype][jtype]) {
|
||||
ecoul = factor_coul * qqrd2e * qtmp*q[j]*(etmp+eps[j])*sqrt(r2inv);
|
||||
} else ecoul = 0.0;
|
||||
ecoul = factor_coul * qqrd2e * qtmp * q[j] * (etmp + eps[j]) * sqrt(r2inv);
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
ecoul *= 0.5;
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, ecoul, fpair, delx, dely, delz, thr);
|
||||
}
|
||||
}
|
||||
f[i].x += fxtmp;
|
||||
@ -227,4 +232,3 @@ void PairLJCutCoulCutDielectricOMP::eval(int iifrom, int iito, ThrData * const t
|
||||
efield[i][2] += eztmp;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -20,9 +20,9 @@
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neigh_list.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
@ -30,41 +30,37 @@
|
||||
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_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
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLongDielectricOMP::PairLJCutCoulLongDielectricOMP(LAMMPS *lmp) :
|
||||
PairLJCutCoulLongDielectric(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
PairLJCutCoulLongDielectric(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJCutCoulLongDielectricOMP::~PairLJCutCoulLongDielectricOMP()
|
||||
{
|
||||
}
|
||||
PairLJCutCoulLongDielectricOMP::~PairLJCutCoulLongDielectricOMP() {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJCutCoulLongDielectricOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
if (eflag || vflag) {
|
||||
ev_setup(eflag,vflag);
|
||||
} else evflag = vflag_fdotr = 0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
if (atom->nmax > nmax) {
|
||||
memory->destroy(efield);
|
||||
memory->destroy(epot);
|
||||
nmax = atom->nmax;
|
||||
memory->create(efield,nmax,3,"pair:efield");
|
||||
memory->create(epot,nmax,"pair:epot");
|
||||
memory->create(efield, nmax, 3, "pair:efield");
|
||||
memory->create(epot, nmax, "pair:epot");
|
||||
}
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
@ -72,7 +68,7 @@ void PairLJCutCoulLongDielectricOMP::compute(int eflag, int vflag)
|
||||
const int inum = list->inum;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel default(none) shared(eflag,vflag)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
@ -84,50 +80,55 @@ void PairLJCutCoulLongDielectricOMP::compute(int eflag, int vflag)
|
||||
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
|
||||
else eval<1,1,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<1, 1, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 1, 0>(ifrom, ito, thr);
|
||||
} else {
|
||||
if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
|
||||
else eval<1,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<1, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
} else {
|
||||
if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
|
||||
else eval<0,0,0>(ifrom, ito, thr);
|
||||
if (force->newton_pair)
|
||||
eval<0, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<0, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
|
||||
thr->timer(Timer::PAIR);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData * const thr)
|
||||
void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData *const thr)
|
||||
{
|
||||
int i,j,ii,jj,jnum,itype,jtype,itable;
|
||||
double qtmp,etmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
|
||||
double fraction,table;
|
||||
double r,rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
|
||||
double grij,expm2,prefactor,t,erfc,prefactorE,efield_i,epot_i;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
int i, j, ii, jj, jnum, itype, jtype, itable;
|
||||
double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
|
||||
double fraction, table;
|
||||
double r, rsq, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj;
|
||||
double grij, expm2, prefactor, t, erfc, prefactorE, efield_i, epot_i;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
|
||||
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const double * _noalias const q = atom->q;
|
||||
const double * _noalias const eps = atom->epsilon;
|
||||
const dbl3_t * _noalias const norm = (dbl3_t *) atom->mu[0];
|
||||
const double * _noalias const curvature = atom->curvature;
|
||||
const double * _noalias const area = atom->area;
|
||||
const int * _noalias const type = atom->type;
|
||||
const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0];
|
||||
dbl3_t *_noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const double *_noalias const q = atom->q;
|
||||
const double *_noalias const eps = atom->epsilon;
|
||||
const dbl3_t *_noalias const norm = (dbl3_t *) atom->mu[0];
|
||||
const double *_noalias const curvature = atom->curvature;
|
||||
const double *_noalias const area = atom->area;
|
||||
const int *_noalias const type = atom->type;
|
||||
const int nlocal = atom->nlocal;
|
||||
const double * _noalias const special_coul = force->special_coul;
|
||||
const double * _noalias const special_lj = force->special_lj;
|
||||
const double *_noalias const special_coul = force->special_coul;
|
||||
const double *_noalias const special_lj = force->special_lj;
|
||||
const double qqrd2e = force->qqrd2e;
|
||||
double fxtmp,fytmp,fztmp,extmp,eytmp,eztmp;
|
||||
|
||||
double fxtmp, fytmp, fztmp, extmp, eytmp, eztmp;
|
||||
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
@ -146,17 +147,17 @@ void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData * const
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp=fytmp=fztmp=0.0;
|
||||
extmp=eytmp=eztmp=0.0;
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
extmp = eytmp = eztmp = 0.0;
|
||||
|
||||
// self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al.
|
||||
|
||||
double curvature_threshold = sqrt(area[i]);
|
||||
if (curvature[i] < curvature_threshold) {
|
||||
double sf = curvature[i]/(4.0*MY_PIS*curvature_threshold) * area[i]*q[i];
|
||||
efield[i][0] = sf*norm[i].x;
|
||||
efield[i][1] = sf*norm[i].y;
|
||||
efield[i][2] = sf*norm[i].z;
|
||||
double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
|
||||
efield[i][0] = sf * norm[i].x;
|
||||
efield[i][1] = sf * norm[i].y;
|
||||
efield[i][2] = sf * norm[i].z;
|
||||
} else {
|
||||
efield[i][0] = efield[i][1] = efield[i][2] = 0;
|
||||
}
|
||||
@ -170,26 +171,26 @@ void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData * const
|
||||
delx = xtmp - x[j].x;
|
||||
dely = ytmp - x[j].y;
|
||||
delz = ztmp - x[j].z;
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r2inv = 1.0 / rsq;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
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;
|
||||
prefactor = qqrd2e * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
expm2 = exp(-grij * grij);
|
||||
t = 1.0 / (1.0 + EWALD_P * grij);
|
||||
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
|
||||
prefactor = qqrd2e * qtmp * q[j] / r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j]/r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F*grij*expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] / r;
|
||||
efield_i = prefactorE * (erfc + EWALD_F * grij * expm2);
|
||||
if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
epot_i = efield_i;
|
||||
|
||||
} else {
|
||||
@ -198,65 +199,68 @@ void PairLJCutCoulLongDielectricOMP::eval(int iifrom, int iito, ThrData * const
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = qtmp*q[j] * table;
|
||||
efield_i = q[j] * table/qqrd2e;
|
||||
table = ftable[itable] + fraction * dftable[itable];
|
||||
forcecoul = qtmp * q[j] * table;
|
||||
efield_i = q[j] * table / qqrd2e;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
table = ctable[itable] + fraction * dctable[itable];
|
||||
prefactor = qtmp * q[j] * table;
|
||||
forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
|
||||
prefactorE = q[j] * table/qqrd2e;
|
||||
efield_i -= (1.0-factor_coul)*prefactorE;
|
||||
prefactorE = q[j] * table / qqrd2e;
|
||||
efield_i -= (1.0 - factor_coul) * prefactorE;
|
||||
}
|
||||
epot_i = efield_i;
|
||||
}
|
||||
} else epot_i = efield_i = forcecoul = 0.0;
|
||||
} else
|
||||
epot_i = efield_i = forcecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
r6inv = r2inv*r2inv*r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
|
||||
r6inv = r2inv * r2inv * r2inv;
|
||||
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
|
||||
forcelj *= factor_lj;
|
||||
} else forcelj = 0.0;
|
||||
} else
|
||||
forcelj = 0.0;
|
||||
|
||||
fpair = (forcecoul + forcelj) * r2inv;
|
||||
|
||||
fxtmp += delx*fpair;
|
||||
fytmp += dely*fpair;
|
||||
fztmp += delz*fpair;
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
|
||||
efield_i *= (etmp*r2inv);
|
||||
extmp += delx*efield_i;
|
||||
eytmp += dely*efield_i;
|
||||
eztmp += delz*efield_i;
|
||||
efield_i *= (etmp * r2inv);
|
||||
extmp += delx * efield_i;
|
||||
eytmp += dely * efield_i;
|
||||
eztmp += delz * efield_i;
|
||||
epot[i] += epot_i;
|
||||
|
||||
if (NEWTON_PAIR || j < nlocal) {
|
||||
f[j].x -= delx*fpair;
|
||||
f[j].y -= dely*fpair;
|
||||
f[j].z -= delz*fpair;
|
||||
f[j].x -= delx * fpair;
|
||||
f[j].y -= dely * fpair;
|
||||
f[j].z -= delz * fpair;
|
||||
}
|
||||
|
||||
if (EFLAG) {
|
||||
if (rsq < cut_coulsq) {
|
||||
if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(etmp+eps[j])*erfc;
|
||||
ecoul = prefactor * (etmp + eps[j]) * erfc;
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = qtmp*q[j]*(etmp+eps[j]) * table;
|
||||
table = etable[itable] + fraction * detable[itable];
|
||||
ecoul = qtmp * q[j] * (etmp + eps[j]) * table;
|
||||
}
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
} else ecoul = 0.0;
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
|
||||
if (rsq < cut_ljsq[itype][jtype]) {
|
||||
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
|
||||
offset[itype][jtype];
|
||||
evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
} else evdwl = 0.0;
|
||||
} else
|
||||
evdwl = 0.0;
|
||||
}
|
||||
|
||||
if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
|
||||
evdwl,ecoul,fpair,delx,dely,delz,thr);
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, ecoul, fpair, delx, dely, delz, thr);
|
||||
}
|
||||
}
|
||||
f[i].x += fxtmp;
|
||||
|
||||
Reference in New Issue
Block a user