diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp index 1c2a99eb19..10876459fd 100644 --- a/src/USER-DIELECTRIC/atom_vec_dielectric.cpp +++ b/src/USER-DIELECTRIC/atom_vec_dielectric.cpp @@ -11,16 +11,8 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include -#include #include "atom_vec_dielectric.h" #include "atom.h" -#include "comm.h" -#include "domain.h" -#include "modify.h" -#include "fix.h" -#include "memory.h" -#include "error.h" using namespace LAMMPS_NS; @@ -28,30 +20,67 @@ using namespace LAMMPS_NS; AtomVecDielectric::AtomVecDielectric(LAMMPS *lmp) : AtomVec(lmp) { - molecular = Atom::ATOMIC; + molecular = Atom::MOLECULAR; + bonds_allow = angles_allow = dihedrals_allow = impropers_allow = 1; mass_type = PER_TYPE; - - atom->q_flag = atom->mu_flag = 1; + atom->molecule_flag = atom->q_flag = atom->mu_flag = 1; + // strings with peratom variables to include in each AtomVec method // strings cannot contain fields in corresponding AtomVec default strings // order of fields in a string does not matter // except: fields_data_atom & fields_data_vel must match data file - fields_grow = (char *) "q mu3 area ed em epsilon curvature"; - fields_copy = (char *) "q mu3 area ed em epsilon curvature"; - fields_comm = (char *) "q mu3 area ed em epsilon curvature"; - fields_comm_vel = (char *) "q mu3 area ed em epsilon curvature"; + fields_grow = (char *) + "q molecule num_bond bond_type bond_atom " + "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 " + "num_dihedral dihedral_type dihedral_atom1 dihedral_atom2 " + "dihedral_atom3 dihedral_atom4 " + "num_improper improper_type improper_atom1 improper_atom2 " + "improper_atom3 improper_atom4 " + "nspecial special " + "mu3 area ed em epsilon curvature q_unscaled"; + fields_copy = (char *) + "q molecule num_bond bond_type bond_atom " + "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 " + "num_dihedral dihedral_type dihedral_atom1 dihedral_atom2 " + "dihedral_atom3 dihedral_atom4 " + "num_improper improper_type improper_atom1 improper_atom2 " + "improper_atom3 improper_atom4 " + "nspecial special " + "mu3 area ed em epsilon curvature q_unscaled"; + fields_comm = (char *) ""; + fields_comm_vel = (char *) ""; fields_reverse = (char *) ""; - fields_border = (char *) "q mu3 area ed em epsilon curvature"; - fields_border_vel = (char *) "q mu3 area ed em epsilon curvature"; - fields_exchange = (char *) "q mu3 area ed em epsilon curvature"; - fields_restart = (char * ) "q mu3 area ed em epsilon curvature"; - fields_create = (char *) "q mu3 area ed em epsilon curvature"; - fields_data_atom = (char *) "id molecule type q x mu3 area ed em epsilon curvature"; + fields_border = (char *) "q molecule mu3 area ed em epsilon curvature"; + fields_border_vel = (char *) "q molecule mu3 area ed em epsilon curvature"; + fields_exchange = (char *) + "q molecule num_bond bond_type bond_atom " + "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 " + "num_dihedral dihedral_type dihedral_atom1 dihedral_atom2 " + "dihedral_atom3 dihedral_atom4 " + "num_improper improper_type improper_atom1 improper_atom2 " + "improper_atom3 improper_atom4 " + "nspecial special " + "mu3 area ed em epsilon curvature q_unscaled"; + fields_restart = (char *) + "q molecule num_bond bond_type bond_atom " + "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 " + "num_dihedral dihedral_type dihedral_atom1 dihedral_atom2 " + "dihedral_atom3 dihedral_atom4 " + "num_improper improper_type improper_atom1 improper_atom2 " + "improper_atom3 improper_atom4 " + "mu3 area ed em epsilon curvature q_unscaled"; + fields_create = (char *) + "q molecule num_bond num_angle num_dihedral num_improper nspecial " + "mu3 area ed em epsilon curvature q_unscaled"; + fields_data_atom = (char *) "id molecule type q x " + "mu3 area ed em epsilon curvature"; fields_data_vel = (char *) "id v"; setup_fields(); + + bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0; } /* ---------------------------------------------------------------------- @@ -61,6 +90,16 @@ AtomVecDielectric::AtomVecDielectric(LAMMPS *lmp) : AtomVec(lmp) void AtomVecDielectric::grow_pointers() { + num_bond = atom->num_bond; + bond_type = atom->bond_type; + num_angle = atom->num_angle; + angle_type = atom->angle_type; + num_dihedral = atom->num_dihedral; + dihedral_type = atom->dihedral_type; + num_improper = atom->num_improper; + improper_type = atom->improper_type; + nspecial = atom->nspecial; + mu = atom->mu; area = atom->area; ed = atom->ed; @@ -88,6 +127,14 @@ void AtomVecDielectric::create_atom_post(int ilocal) void AtomVecDielectric::data_atom_post(int ilocal) { + num_bond[ilocal] = 0; + num_angle[ilocal] = 0; + num_dihedral[ilocal] = 0; + num_improper[ilocal] = 0; + nspecial[ilocal][0] = 0; + nspecial[ilocal][1] = 0; + nspecial[ilocal][2] = 0; + double* q = atom->q; q_unscaled[ilocal] = q[ilocal]; q[ilocal] /= epsilon[ilocal]; @@ -103,6 +150,14 @@ void AtomVecDielectric::pack_data_post(int ilocal) q_unscaled[ilocal] = q[ilocal]/epsilon[ilocal]; } +/* ---------------------------------------------------------------------- + initialize other atom quantities after AtomVec::unpack_restart() +------------------------------------------------------------------------- */ - +void AtomVecDielectric::unpack_restart_init(int ilocal) +{ + nspecial[ilocal][0] = 0; + nspecial[ilocal][1] = 0; + nspecial[ilocal][2] = 0; +} diff --git a/src/USER-DIELECTRIC/atom_vec_dielectric.h b/src/USER-DIELECTRIC/atom_vec_dielectric.h index 797e41f112..6c933c5949 100644 --- a/src/USER-DIELECTRIC/atom_vec_dielectric.h +++ b/src/USER-DIELECTRIC/atom_vec_dielectric.h @@ -32,12 +32,17 @@ class AtomVecDielectric : public AtomVec { void create_atom_post(int); void data_atom_post(int); void pack_data_post(int); + void unpack_restart_init(int); + + private: + 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; double **mu; double *area,*ed,*em,*epsilon,*curvature,*q_unscaled; - - private: - }; } diff --git a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp index 2cb411e68f..724fac2fac 100644 --- a/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp @@ -15,21 +15,21 @@ Contributing author: Trung Nguyen (Northwestern) ------------------------------------------------------------------------- */ -#include -#include -#include -#include #include "pair_coul_cut_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 "memory.h" #include "math_const.h" -#include "error.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -74,11 +74,11 @@ void PairCoulCutDielectric::compute(int eflag, int vflag) double **x = atom->x; double **f = atom->f; double *q = atom->q; - double *q_real = avec->q_unscaled; - double* eps = avec->epsilon; - double** norm = avec->mu; - double* curvature = avec->curvature; - double* area = avec->area; + double *q_real = atom->q_unscaled; + 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; @@ -183,7 +183,7 @@ double PairCoulCutDielectric::single(int i, int j, int itype, int jtype, double &fforce) { double r2inv,forcecoul,phicoul,ei,ej; - double* eps = avec->epsilon; + double* eps = atom->epsilon; r2inv = 1.0/rsq; forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv)*eps[i]; 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 56bdd579d9..8a60ccd31b 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp @@ -208,7 +208,7 @@ double PairLJCutCoulCutDielectric::single(int i, int j, int itype, int jtype, double &fforce) { double r2inv,r6inv,forcecoul,forcelj,phicoul,ei,ej,philj; - double* eps = avec->epsilon; + double* eps = atom->epsilon; r2inv = 1.0/rsq; if (rsq < cut_coulsq[itype][jtype]) 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 613defa2f6..4967358aa4 100644 --- a/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp @@ -277,7 +277,7 @@ double PairLJCutCoulLongDielectric::single(int i, int j, int itype, int jtype, double r2inv,r6inv,r,grij,expm2,t,erfc,ei,ej,prefactor; double fraction,table,forcecoul,forcelj,phicoul,philj; int itable; - double *eps = avec->epsilon; + double *eps = atom->epsilon; r2inv = 1.0/rsq; if (rsq < cut_coulsq) { diff --git a/src/USER-DIELECTRIC/pppm_dielectric.cpp b/src/USER-DIELECTRIC/pppm_dielectric.cpp index 62779d79c4..d9c957aaa3 100644 --- a/src/USER-DIELECTRIC/pppm_dielectric.cpp +++ b/src/USER-DIELECTRIC/pppm_dielectric.cpp @@ -482,9 +482,6 @@ void PPPMDielectric::qsum_qsq() const int nlocal = atom->nlocal; double qsum_local(0.0), qsqsum_local(0.0); -#if defined(_OPENMP) -#pragma omp parallel for default(none) reduction(+:qsum_local,qsqsum_local) -#endif for (int i = 0; i < nlocal; i++) { double qtmp = eps[i]*q[i]; qsum_local += qtmp; diff --git a/src/atom.cpp b/src/atom.cpp index 63e90692e0..75181abbd7 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -520,7 +520,7 @@ void Atom::peratom_create() add_peratom("em",&em,DOUBLE,0); add_peratom("epsilon",&epsilon,DOUBLE,0); add_peratom("curvature",&curvature,DOUBLE,0); - add_peratom("q_unscaled",&curvature,DOUBLE,0); + add_peratom("q_unscaled",&q_unscaled,DOUBLE,0); // end of customization section // --------------------------------------------------------------------