Updated AtomVecDielectric to support atom properties like full (molecule, bonds, angles, dihedrals, special)

src/USER-DIELECTRIC/atom_vec_dielectric.cpp

@@ -11,16 +11,8 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
-#include <cmath>
-#include <cstdlib>
 #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_grow = (char *)
-  fields_copy = (char *) "q mu3 area ed em epsilon curvature";
+    "q molecule num_bond bond_type bond_atom "
-  fields_comm = (char *) "q mu3 area ed em epsilon curvature";
+    "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 "
-  fields_comm_vel = (char *) "q mu3 area ed em epsilon curvature";
+    "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 = (char *) "q molecule mu3 area ed em epsilon curvature";
-  fields_border_vel = (char *) "q mu3 area ed em epsilon curvature";
+  fields_border_vel = (char *) "q molecule mu3 area ed em epsilon curvature";
-  fields_exchange = (char *) "q mu3 area ed em epsilon curvature";
+  fields_exchange = (char *)
-  fields_restart = (char * ) "q mu3 area ed em epsilon curvature";
+    "q molecule num_bond bond_type bond_atom "
-  fields_create = (char *) "q mu3 area ed em epsilon curvature";
+    "num_angle angle_type angle_atom1 angle_atom2 angle_atom3 "
-  fields_data_atom = (char *) "id molecule type q x mu3 area ed em epsilon curvature";
+    "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;
+}
 

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:
-  
 };
 
 }

src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp

@@ -15,21 +15,21 @@
    Contributing author: Trung Nguyen (Northwestern)
 ------------------------------------------------------------------------- */
 
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
 #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 <cmath>
+#include <cstring>
 
 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 *q_real = atom->q_unscaled;
-  double* eps = avec->epsilon;
+  double* eps = atom->epsilon;
-  double** norm = avec->mu;
+  double** norm = atom->mu;
-  double* curvature = avec->curvature;
+  double* curvature = atom->curvature;
-  double* area = avec->area;
+  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];

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])

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) {

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;

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
   // --------------------------------------------------------------------