Add support for "Dipoles" section in molecule file

doc/src/molecule.rst

@@ -188,7 +188,7 @@ internally.
 
 These are the allowed section keywords for the body of the file.
 
-* *Coords, Types, Molecules, Fragments, Charges, Diameters, Masses* = atom-property sections
+* *Coords, Types, Molecules, Fragments, Charges, Diameters, Dipoles, Masses* = atom-property sections
 * *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections
 * *Special Bond Counts, Special Bonds* = special neighbor info
 * *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info
@@ -303,6 +303,19 @@ not listed, the default diameter of each atom in the molecule is 1.0.
 
 ----------
 
+*Dipoles* section:
+
+* one line per atom
+* line syntax: ID mux muy muz
+* mux,muy,muz = x-, y-, and z-component of point dipole vector of atom
+
+This section is only allowed for :doc:`atom styles <atom_style>` that
+support particles with point dipoles, e.g. atom_style dipole.  If not
+listed, the default dipole component of each atom in the molecule is set
+to 0.0.
+
+----------
+
 *Masses* section:
 
 * one line per atom

src/atom.cpp

@@ -26,6 +26,7 @@
 #include "input.h"
 #include "label_map.h"
 #include "math_const.h"
+#include "math_extra.h"
 #include "memory.h"
 #include "modify.h"
 #include "molecule.h"
@@ -2112,6 +2113,15 @@ std::vector<Molecule *>Atom::get_molecule_by_id(const std::string &id)
 void Atom::add_molecule_atom(Molecule *onemol, int iatom, int ilocal, tagint offset)
 {
   if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
+  if (onemol->muflag && mu_flag) {
+    double r[3], rotmat[3][3];
+    MathExtra::quat_to_mat(onemol->quat_external, rotmat);
+    MathExtra::matvec(rotmat, onemol->mu[iatom], r);
+    mu[ilocal][0] = r[0];
+    mu[ilocal][1] = r[1];
+    mu[ilocal][2] = r[2];
+    mu[ilocal][3] = sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
+  }
   if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
   if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
   else if (rmass_flag)

src/molecule.cpp

@@ -41,15 +41,16 @@ using namespace LAMMPS_NS;
 
 Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
     Pointers(lmp), id(nullptr), x(nullptr), type(nullptr), molecule(nullptr), q(nullptr),
-    radius(nullptr), rmass(nullptr), num_bond(nullptr), bond_type(nullptr), bond_atom(nullptr),
+    radius(nullptr), rmass(nullptr), mu(nullptr), num_bond(nullptr), bond_type(nullptr),
-    num_angle(nullptr), angle_type(nullptr), angle_atom1(nullptr), angle_atom2(nullptr),
+    bond_atom(nullptr), num_angle(nullptr), angle_type(nullptr), angle_atom1(nullptr),
-    angle_atom3(nullptr), num_dihedral(nullptr), dihedral_type(nullptr), dihedral_atom1(nullptr),
+    angle_atom2(nullptr), angle_atom3(nullptr), num_dihedral(nullptr), dihedral_type(nullptr),
-    dihedral_atom2(nullptr), dihedral_atom3(nullptr), dihedral_atom4(nullptr),
+    dihedral_atom1(nullptr), dihedral_atom2(nullptr), dihedral_atom3(nullptr),
-    num_improper(nullptr), improper_type(nullptr), improper_atom1(nullptr), improper_atom2(nullptr),
+    dihedral_atom4(nullptr), num_improper(nullptr), improper_type(nullptr), improper_atom1(nullptr),
-    improper_atom3(nullptr), improper_atom4(nullptr), nspecial(nullptr), special(nullptr),
+    improper_atom2(nullptr), improper_atom3(nullptr), improper_atom4(nullptr), nspecial(nullptr),
-    shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), avec_body(nullptr),
+    special(nullptr), shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr),
-    ibodyparams(nullptr), dbodyparams(nullptr), fragmentmask(nullptr), dx(nullptr), dxcom(nullptr),
+    avec_body(nullptr), ibodyparams(nullptr), dbodyparams(nullptr), fragmentmask(nullptr),
-    dxbody(nullptr), quat_external(nullptr), fp(nullptr), count(nullptr)
+    dx(nullptr), dxcom(nullptr), dxbody(nullptr), quat_external(nullptr), fp(nullptr),
+    count(nullptr)
 {
   me = comm->me;
 
@@ -132,7 +133,7 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
 
   // initialize all fields to empty
 
-  initialize();
+  Molecule::initialize();
 
   // scan file for sizes of all fields and allocate storage for them
 
@@ -141,14 +142,14 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
     if (fp == nullptr)
       error->one(FLERR, "Cannot open molecule file {}: {}", arg[ifile], utils::getsyserror());
   }
-  read(0);
+  Molecule::read(0);
   if (me == 0) fclose(fp);
-  allocate();
+  Molecule::allocate();
 
   // read file again to populate all fields
 
   if (me == 0) fp = fopen(arg[ifile], "r");
-  read(1);
+  Molecule::read(1);
   if (me == 0) fclose(fp);
 
   // stats
@@ -572,6 +573,12 @@ void Molecule::read(int flag)
         diameters(line);
       else
         skip_lines(natoms, line, keyword);
+    } else if (keyword == "Dipoles") {
+      muflag = 1;
+      if (flag)
+        dipoles(line);
+      else
+        skip_lines(natoms, line, keyword);
     } else if (keyword == "Masses") {
       rmassflag = 1;
       if (flag)
@@ -948,6 +955,40 @@ void Molecule::diameters(char *line)
   }
 }
 
+/* ----------------------------------------------------------------------
+   read charges from file
+------------------------------------------------------------------------- */
+
+void Molecule::dipoles(char *line)
+{
+  for (int i = 0; i < natoms; i++) count[i] = 0;
+  try {
+    for (int i = 0; i < natoms; i++) {
+      readline(line);
+
+      ValueTokenizer values(utils::trim_comment(line));
+      if ((int) values.count() != 4)
+        error->all(FLERR, "Invalid line in Dipoles section of molecule file: {}", line);
+
+      int iatom = values.next_int() - 1;
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR, "Invalid atom index in Dipoles section of molecule file");
+
+      count[iatom]++;
+      mu[iatom][0] = values.next_double();
+      mu[iatom][1] = values.next_double();
+      mu[iatom][2] = values.next_double();
+    }
+  } catch (TokenizerException &e) {
+    error->all(FLERR, "Invalid line in Dipoles section of molecule file: {}\n{}", e.what(), line);
+  }
+
+  for (int i = 0; i < natoms; i++) {
+    if (count[i] == 0)
+      error->all(FLERR, "Atom {} missing in Dipoles section of molecule file", i + 1);
+  }
+}
+
 /* ----------------------------------------------------------------------
    read masses from file
 ------------------------------------------------------------------------- */
@@ -1828,6 +1869,7 @@ void Molecule::check_attributes()
 
   int mismatch = 0;
   if (qflag && !atom->q_flag) mismatch = 1;
+  if (muflag && !atom->mu_flag) mismatch = 1;
   if (radiusflag && !atom->radius_flag) mismatch = 1;
   if (rmassflag && !atom->rmass_flag) mismatch = 1;
 
@@ -1880,7 +1922,7 @@ void Molecule::initialize()
   bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
   maxspecial = 0;
 
-  xflag = typeflag = moleculeflag = fragmentflag = qflag = radiusflag = rmassflag = 0;
+  xflag = typeflag = moleculeflag = fragmentflag = qflag = radiusflag = muflag = rmassflag = 0;
   bondflag = angleflag = dihedralflag = improperflag = 0;
   nspecialflag = specialflag = 0;
   shakeflag = shakeflagflag = shakeatomflag = shaketypeflag = 0;
@@ -1943,6 +1985,7 @@ void Molecule::allocate()
       for (int j = 0; j < natoms; j++) fragmentmask[i][j] = 0;
   }
   if (qflag) memory->create(q, natoms, "molecule:q");
+  if (muflag) memory->create(mu, natoms, 3, "molecule:mu");
   if (radiusflag) memory->create(radius, natoms, "molecule:radius");
   if (rmassflag) memory->create(rmass, natoms, "molecule:rmass");
 
@@ -2167,6 +2210,11 @@ void Molecule::print()
     for (int i = 0; i < natoms; i++)
       printf("    %d %g\n",i+1,radius[i]);
   }
+  if (muflag) {
+    printf(  "Dipoles:\n");
+    for (int i = 0; i < natoms; i++)
+      printf("    %d %g %g %g\n",i+1,mu[i][0],mu[i][1],mu[i][2]);
+  }
   if (rmassflag) {
     printf(  "Masses:\n");
     for (int i = 0; i < natoms; i++)

src/molecule.h

@@ -41,7 +41,7 @@ class Molecule : protected Pointers {
 
   // 1 if attribute defined in file, 0 if not
 
-  int xflag, typeflag, moleculeflag, fragmentflag, qflag, radiusflag, rmassflag;
+  int xflag, typeflag, moleculeflag, fragmentflag, qflag, radiusflag, muflag, rmassflag;
   int bondflag, angleflag, dihedralflag, improperflag;
   int nspecialflag, specialflag;
   int shakeflag, shakeflagflag, shakeatomflag, shaketypeflag;
@@ -63,6 +63,7 @@ class Molecule : protected Pointers {
   double *q;           // charge on each atom
   double *radius;      // radius of each atom
   double *rmass;       // mass of each atom
+  double **mu;         // dipole vector of each atom
 
   int *num_bond;    // bonds, angles, dihedrals, impropers for each atom
   int **bond_type;
@@ -142,6 +143,7 @@ class Molecule : protected Pointers {
   void fragments(char *);
   void charges(char *);
   void diameters(char *);
+  void dipoles(char *);
   void masses(char *);
   void bonds(int, char *);
   void angles(int, char *);