diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp index c804b05ed3..b0e3487807 100644 --- a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp @@ -13,6 +13,7 @@ ------------------------------------------------------------------------- */ #include "atom_vec_dielectric.h" + #include "atom.h" #include "citeme.h" diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.h b/src/USER-DIELECTRIC/atom_vec_dielectric.h index dbdd4ff53b..8e342eda08 100644 --- a/src/USER-DIELECTRIC/atom_vec_dielectric.h +++ b/src/USER-DIELECTRIC/atom_vec_dielectric.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 diff --git a/src/USER-DIELECTRIC/compute_efield_atom.h b/src/USER-DIELECTRIC/compute_efield_atom.h index 6c81150265..1961e15d69 100644 --- a/src/USER-DIELECTRIC/compute_efield_atom.h +++ b/src/USER-DIELECTRIC/compute_efield_atom.h @@ -44,7 +44,7 @@ class ComputeEfieldAtom : public Compute { double **efield; }; -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/USER-DIELECTRIC/fix_polarize_bem_gmres.h b/src/USER-DIELECTRIC/fix_polarize_bem_gmres.h index e680b63240..7671ddf1c2 100644 --- a/src/USER-DIELECTRIC/fix_polarize_bem_gmres.h +++ b/src/USER-DIELECTRIC/fix_polarize_bem_gmres.h @@ -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 diff --git a/src/USER-DIELECTRIC/fix_polarize_bem_icc.h b/src/USER-DIELECTRIC/fix_polarize_bem_icc.h index eca8ffcaee..6ba6acb1a0 100644 --- a/src/USER-DIELECTRIC/fix_polarize_bem_icc.h +++ b/src/USER-DIELECTRIC/fix_polarize_bem_icc.h @@ -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 diff --git a/src/USER-DIELECTRIC/fix_polarize_functional.h b/src/USER-DIELECTRIC/fix_polarize_functional.h index 704e976020..7a4b8c00da 100644 --- a/src/USER-DIELECTRIC/fix_polarize_functional.h +++ b/src/USER-DIELECTRIC/fix_polarize_functional.h @@ -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 diff --git a/src/USER-DIELECTRIC/msm_dielectric.h b/src/USER-DIELECTRIC/msm_dielectric.h index 1a36386f13..9874f5e0b2 100644 --- a/src/USER-DIELECTRIC/msm_dielectric.h +++ b/src/USER-DIELECTRIC/msm_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp index b6983586cf..1b540157e4 100644 --- a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -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; diff --git a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h index e10df64691..ef6fd3fb45 100644 --- a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h +++ b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp index a69ba9fd00..fd5110cdb6 100644 --- a/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp @@ -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 #include @@ -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); } - diff --git a/src/USER-DIELECTRIC/pair_coul_long_dielectric.h b/src/USER-DIELECTRIC/pair_coul_long_dielectric.h index 7f3bfebf5a..993555e452 100644 --- a/src/USER-DIELECTRIC/pair_coul_long_dielectric.h +++ b/src/USER-DIELECTRIC/pair_coul_long_dielectric.h @@ -34,12 +34,11 @@ class PairCoulLongDielectric : public PairCoulLong { double **efield; protected: - class AtomVecDielectric* avec; + class AtomVecDielectric *avec; int nmax; - }; -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp index e03002c3bd..82a5b0c99a 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -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 @@ -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; diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.h b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.h index 52e6ce1789..93bee549c8 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.h +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp index d1485bd32d..4f8b2d63c8 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.cpp @@ -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 -#include -#include -#include #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 +#include 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; diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h b/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h index 18b3c1bb59..0301ff4fbd 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_debye_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp index ecbb592e5c..7f5b67a22b 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -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 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; diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.h b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.h index 1d2b3c6123..06ad1c5e7f 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.h +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp index 4c5efe7c99..c273fc2132 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp @@ -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 @@ -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; } diff --git a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.h b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.h index e64409a2cb..3e9e06bbc9 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.h +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index 09c8de7e71..564a156a04 100644 --- a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -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 -#include -#include -#include -#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 +#include 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 (; ineighfirstneigh[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= 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; } diff --git a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.h b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.h index 0ec3321df4..83ba66a9a0 100644 --- a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.h +++ b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.h @@ -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 diff --git a/src/USER-DIELECTRIC/pppm_dielectric.cpp b/src/USER-DIELECTRIC/pppm_dielectric.cpp index 6d0f3f33b5..7fdb0b76c4 100644 --- a/src/USER-DIELECTRIC/pppm_dielectric.cpp +++ b/src/USER-DIELECTRIC/pppm_dielectric.cpp @@ -34,7 +34,6 @@ #include "remap_wrap.h" #include -//#include using namespace LAMMPS_NS; using namespace MathConst; diff --git a/src/USER-DIELECTRIC/pppm_dielectric.h b/src/USER-DIELECTRIC/pppm_dielectric.h index 24eb180c5f..7ea6a831eb 100644 --- a/src/USER-DIELECTRIC/pppm_dielectric.h +++ b/src/USER-DIELECTRIC/pppm_dielectric.h @@ -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 diff --git a/src/USER-OMP/pair_lj_cut_coul_cut_dielectric_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_cut_dielectric_omp.cpp index a056ae86c9..22ae388236 100644 --- a/src/USER-OMP/pair_lj_cut_coul_cut_dielectric_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_cut_dielectric_omp.cpp @@ -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 +#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 -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; } } - diff --git a/src/USER-OMP/pair_lj_cut_coul_long_dielectric_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_dielectric_omp.cpp index 39149e692f..6aa79e18a5 100644 --- a/src/USER-OMP/pair_lj_cut_coul_long_dielectric_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_long_dielectric_omp.cpp @@ -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 @@ -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 -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;