Merge pull request #4046 from akohlmey/molecule-dipole

Add support for Dipoles section to molecule files

doc/src/molecule.rst

@@ -126,14 +126,50 @@ molecule (header keyword = inertia).
 Format of a molecule file
 """""""""""""""""""""""""
 
-The format of an individual molecule file is similar but
+The format of an individual molecule file looks similar but is
-(not identical) to the data file read by the :doc:`read_data <read_data>`
+different than that of a data file read by the :doc:`read_data <read_data>`
-commands, and is as follows.
+commands.  Here is a simple example for a TIP3P water molecule:
+
+.. code-block::
+
+   # Water molecule. TIP3P geometry
+   # header section:
+   3 atoms
+   2 bonds
+   1 angles
+
+   # body section:
+   Coords
+
+   1    0.00000  -0.06556   0.00000
+   2    0.75695   0.52032   0.00000
+   3   -0.75695   0.52032   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+
+   Charges
+
+   1       -0.834
+   2        0.417
+   3        0.417
+
+   Bonds
+
+   1   1      1      2
+   2   1      1      3
+
+   Angles
+
+   1   1      2      1      3
 
 A molecule file has a header and a body.  The header appears first.  The
-first line of the header and thus of the molecule file is *always* skipped;
+first line of the header and thus of the molecule file is *always*
-it typically contains a description of the file or a comment from the software
+skipped; it typically contains a description of the file or a comment
-that created the file.
+from the software that created the file.
 
 Then lines are read one line at a time.  Lines can have a trailing
 comment starting with '#' that is ignored.  There *must* be at least one
@@ -158,25 +194,62 @@ appear if the value(s) are different than the default, except when
 defining a *body* particle, which requires setting the number of
 *atoms* to 1, and setting the *inertia* in a specific section (see below).
 
-* N *atoms* = # of atoms N in molecule, default = 0
+   .. list-table::
-* Nb *bonds* = # of bonds Nb in molecule, default = 0
+      :header-rows: 1
-* Na *angles* = # of angles Na in molecule, default = 0
+      :widths: auto
-* Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0
-* Ni *impropers* = # of impropers Ni in molecule, default = 0
-* Nf *fragments* = # of fragments Nf in molecule, default = 0
-* Ninteger Ndouble *body* = # of integer and floating-point values
-  in body particle, default = 0
-* Mtotal *mass* = total mass of molecule
-* Xc Yc Zc *com* = coordinates of center-of-mass of molecule
-* Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule
 
-For *mass*, *com*, and *inertia*, the default is for LAMMPS to
+      * - Number(s)
-calculate this quantity itself if needed, assuming the molecules
+        - Keyword
-consist of a set of point particles or finite-size particles (with a
+        - Meaning
-non-zero diameter) that do not overlap.  If finite-size particles in
+        - Default Value
-the molecule do overlap, LAMMPS will not account for the overlap
+      * - N
-effects when calculating any of these 3 quantities, so you should
+        - atoms
-pre-compute them yourself and list the values in the file.
+        - # of atoms N in molecule
+        - 0
+      * - Nb
+        - bonds
+        - # of bonds Nb in molecule
+        - 0
+      * - Na
+        - angles
+        - # of angles Na in molecule
+        - 0
+      * - Nd
+        - dihedrals
+        - # of dihedrals Nd in molecule
+        - 0
+      * - Ni
+        - impropers
+        - # of impropers Ni in molecule
+        - 0
+      * - Nf
+        - fragments
+        - # of fragments Nf in molecule
+        - 0
+      * - Ninteger Ndouble
+        - body
+        - # of integer and floating-point values in body particle
+        - 0
+      * - Mtotal
+        - mass
+        - total mass of molecule
+        - computed
+      * - Xc Yc Zc
+        - com
+        - coordinates of center-of-mass of molecule
+        - computed
+      * - Ixx Iyy Izz Ixy Ixz Iyz
+        - inertia
+        - 6 components of inertia tensor of molecule
+        - computed
+
+For *mass*, *com*, and *inertia*, the default is for LAMMPS to calculate
+this quantity itself if needed, assuming the molecules consist of a set
+of point particles or finite-size particles (with a non-zero diameter)
+that do **not** overlap.  If finite-size particles in the molecule
+**do** overlap, LAMMPS will not account for the overlap effects when
+calculating any of these 3 quantities, so you should pre-compute them
+yourself and list the values in the file.
 
 The mass and center-of-mass coordinates (Xc,Yc,Zc) are
 self-explanatory.  The 6 moments of inertia (ixx,iyy,izz,ixy,ixz,iyz)
@@ -188,7 +261,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 +376,21 @@ not listed, the default diameter of each atom in the molecule is 1.0.
 
 ----------
 
+..versionadded:: TBD
+
+*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,28 +142,30 @@ 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
 
+  if (title.empty()) title = "(no title)";
   if (me == 0)
     utils::logmesg(lmp,
-                   "Read molecule template {}:\n  {} molecules\n"
+                   "Read molecule template {}:\n{}\n"
+                   "  {} molecules\n"
                    "  {} fragments\n"
                    "  {} atoms with max type {}\n"
                    "  {} bonds with max type {}\n"
                    "  {} angles with max type {}\n"
                    "  {} dihedrals with max type {}\n"
                    "  {} impropers with max type {}\n",
-                   id, nmolecules, nfragments, natoms, ntypes, nbonds, nbondtypes, nangles,
+                   id, title, nmolecules, nfragments, natoms, ntypes, nbonds, nbondtypes, nangles,
                    nangletypes, ndihedrals, ndihedraltypes, nimpropers, nimpropertypes);
 }
 
@@ -423,6 +426,8 @@ void Molecule::read(int flag)
     if (eof == nullptr) error->one(FLERR, "Unexpected end of molecule file");
   }
 
+  if (flag == 0) title = utils::trim(line);
+
   // read header lines
   // skip blank lines or lines that start with "#"
   // stop when read an unrecognized line
@@ -572,6 +577,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 +959,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 +1873,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;
 
@@ -1869,6 +1915,7 @@ void Molecule::check_attributes()
 
 void Molecule::initialize()
 {
+  title.clear();
   natoms = 0;
   nbonds = nangles = ndihedrals = nimpropers = 0;
   ntypes = 0;
@@ -1880,7 +1927,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 +1990,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 +2215,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

@@ -25,6 +25,8 @@ class Molecule : protected Pointers {
                // else 0 if not first in set
   int last;    // 1 if last molecule in set, else 0
 
+  std::string title;    // title string of the molecule file
+
   // number of atoms,bonds,etc in molecule
   // nibody,ndbody = # of integer/double fields in body
 
@@ -41,7 +43,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 +65,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 +145,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 *);

unittest/formats/test_molecule_file.cpp

@@ -32,6 +32,8 @@ using testing::StrEq;
 
 using utils::split_words;
 
+const double EPSILON = 5.0e-14;
+
 #define test_name test_info_->name()
 
 static void create_molecule_files(const std::string &h2o_filename, const std::string &co2_filename)
@@ -145,7 +147,7 @@ protected:
         fclose(fp);
 
         command(fmt::format("molecule {} {} {}", name, file, args));
-        remove(file.c_str());
+        platform::unlink(file.c_str());
     }
 };
 
@@ -184,7 +186,7 @@ TEST_F(MoleculeFileTest, badargs)
     TEST_FAILURE(
         ".*Illegal molecule command.*",
         run_mol_cmd(test_name, "scale", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
-    remove("badargs.mol");
+    platform::unlink("moltest_badargs.mol");
 }
 
 TEST_F(MoleculeFileTest, noatom)
@@ -193,14 +195,14 @@ TEST_F(MoleculeFileTest, noatom)
                  run_mol_cmd(test_name, "",
                              "Comment\n0 atoms\n1 bonds\n\n"
                              " Coords\n\nBonds\n\n 1 1 2\n"););
-    remove("noatom.mol");
+    platform::unlink("moltest_noatom.mol");
 }
 
 TEST_F(MoleculeFileTest, empty)
 {
     TEST_FAILURE(".*ERROR: Unexpected end of molecule file.*",
                  run_mol_cmd(test_name, "", "Comment\n\n"););
-    remove("empty.mol");
+    platform::unlink("moltest_empty.mol");
 }
 
 TEST_F(MoleculeFileTest, nospecial)
@@ -210,7 +212,7 @@ TEST_F(MoleculeFileTest, nospecial)
                              "Comment\n3 atoms\n\n2 bonds\n\n"
                              " Coords\n\n 1 1.0 1.0 1.0\n 2 1.0 1.0 0.0\n 3 1.0 0.0 1.0\n"
                              " Bonds\n\n 1 1 1 2\n 2 1 1 3\n"););
-    remove("nospecial.mol");
+    platform::unlink("moltest_nospecial.mol");
 }
 
 TEST_F(MoleculeFileTest, minimal)
@@ -218,7 +220,7 @@ TEST_F(MoleculeFileTest, minimal)
     BEGIN_CAPTURE_OUTPUT();
     run_mol_cmd(test_name, "", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n");
     auto output = END_CAPTURE_OUTPUT();
-    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*1 molecules.*\n"
+    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*Comment.*\n.*1 molecules.*\n"
                                       ".*0 fragments.*\n.*1 atoms.*\n.*0 bonds.*"));
 }
 
@@ -230,7 +232,7 @@ TEST_F(MoleculeFileTest, notype)
     command("create_box 1 box");
     run_mol_cmd(test_name, "", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n");
     auto output = END_CAPTURE_OUTPUT();
-    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*1 molecules.*\n"
+    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*Comment.*\n.*1 molecules.*\n"
                                       ".*0 fragments.*\n.*1 atoms.*\n.*0 bonds.*"));
     TEST_FAILURE(".*ERROR: Create_atoms molecule must have atom types.*",
                  command("create_atoms 0 single 0.0 0.0 0.0 mol notype 542465"););
@@ -259,7 +261,7 @@ TEST_F(MoleculeFileTest, twomols)
                 " Coords\n\n 1 0.0 0.0 0.0\n 2 0.0 0.0 1.0\n"
                 " Molecules\n\n 1 1\n 2 2\n\n Types\n\n 1 1\n 2 2\n\n");
     auto output = END_CAPTURE_OUTPUT();
-    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*2 molecules.*\n"
+    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*Comment.*\n.*2 molecules.*\n"
                                       ".*0 fragments.*\n.*2 atoms with max type 2.*\n.*0 bonds.*"));
     ASSERT_EQ(lmp->atom->nmolecule, 1);
     auto mols = lmp->atom->get_molecule_by_id(test_name);
@@ -273,10 +275,10 @@ TEST_F(MoleculeFileTest, twofiles)
     auto output = END_CAPTURE_OUTPUT();
     ASSERT_THAT(
         output,
-        ContainsRegex(".*Read molecule template twomols:.*\n.*1 molecules.*\n"
+        ContainsRegex(".*Read molecule template twomols:.*\n.*Water.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 2.*\n.*2 bonds with max type 1.*\n"
                       ".*1 angles with max type 1.*\n.*0 dihedrals.*\n.*0 impropers.*\n"
-                      ".*Read molecule template twomols:.*\n.*1 molecules.*\n"
+                      ".*Read molecule template twomols:.*\n.*CO2.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 4.*\n.*2 bonds with max type 2.*\n"
                       ".*1 angles with max type 2.*\n.*0 dihedrals.*"));
     BEGIN_CAPTURE_OUTPUT();
@@ -306,7 +308,7 @@ TEST_F(MoleculeFileTest, labelmap)
     auto output = END_CAPTURE_OUTPUT();
     ASSERT_THAT(
         output,
-        ContainsRegex(".*Read molecule template h2olabel:.*\n.*1 molecules.*\n"
+        ContainsRegex(".*Read molecule template h2olabel:.*\n.*Water.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 2.*\n.*2 bonds with max type 1.*\n"
                       ".*1 angles with max type 1.*\n.*0 dihedrals.*\n.*0 impropers.*"));
     BEGIN_CAPTURE_OUTPUT();
@@ -314,7 +316,7 @@ TEST_F(MoleculeFileTest, labelmap)
     output = END_CAPTURE_OUTPUT();
     ASSERT_THAT(
         output,
-        ContainsRegex(".*Read molecule template co2label:.*\n.*1 molecules.*\n"
+        ContainsRegex(".*Read molecule template co2label:.*\n.*CO2.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 4.*\n.*2 bonds with max type 2.*\n"
                       ".*1 angles with max type 2.*\n.*0 dihedrals.*"));
     BEGIN_CAPTURE_OUTPUT();
@@ -328,12 +330,12 @@ TEST_F(MoleculeFileTest, labelmap)
     auto second = output.substr(mark);
     ASSERT_THAT(
         first,
-        ContainsRegex(".*Read molecule template h2onum:.*\n.*1 molecules.*\n"
+        ContainsRegex(".*Read molecule template h2onum:.*\n.*Water.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 2.*\n.*2 bonds with max type 1.*\n"
                       ".*1 angles with max type 1.*\n.*0 dihedrals.*\n.*0 impropers.*\n"));
     ASSERT_THAT(
         second,
-        ContainsRegex(".*Read molecule template co2num:.*\n.*1 molecules.*\n"
+        ContainsRegex(".*Read molecule template co2num:.*\n.*CO2.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 4.*\n.*2 bonds with max type 2.*\n"
                       ".*1 angles with max type 2.*\n.*0 dihedrals.*"));
     ASSERT_EQ(lmp->atom->nmolecule, 4);
@@ -379,7 +381,7 @@ TEST_F(MoleculeFileTest, bonds)
                 " 1 1 1 2\n"
                 " 2 2 1 3\n\n");
     auto output = END_CAPTURE_OUTPUT();
-    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*1 molecules.*\n"
+    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*Comment.*\n.*1 molecules.*\n"
                                       ".*0 fragments.*\n.*4 atoms.*type.*2.*\n"
                                       ".*2 bonds.*type.*2.*\n.*0 angles.*"));
 
@@ -404,6 +406,60 @@ TEST_F(MoleculeFileTest, bonds)
     END_HIDE_OUTPUT();
 }
 
+TEST_F(MoleculeFileTest, dipoles)
+{
+    if (!LAMMPS::is_installed_pkg("DIPOLE")) GTEST_SKIP();
+    BEGIN_CAPTURE_OUTPUT();
+    command("atom_style dipole");
+    command("region box block 0 1 0 1 0 1");
+    command("create_box 2 box");
+    run_mol_cmd(test_name, "",
+                "# Dumbbell with dipole molecule file.\n\n"
+                "2 atoms\n\n"
+                "Coords\n\n1 -1.0 0.0 0.0\n2  1.0 0.0 0.0\n\n"
+                "Types\n\n1 1\n2 2\n\n"
+                "Dipoles\n\n1 1.0 0.0 0.0\n2 1.0 1.0 0.0\n\n");
+    auto output = END_CAPTURE_OUTPUT();
+    ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*Dumbbell.*\n.*1 molecules.*\n"
+                                      ".*0 fragments.*\n.*2 atoms.*type.*2.*\n"));
+
+    BEGIN_CAPTURE_OUTPUT();
+    command("mass * 1.0");
+    command("create_atoms 0 single 0.5 0.5 0.5 mol dipoles 67235 rotate 90.0 0.0 0.0 1.0");
+    output = END_CAPTURE_OUTPUT();
+    ASSERT_THAT(output, ContainsRegex(".*Created 2 atoms.*"));
+
+    Molecule *mol = lmp->atom->molecules[0];
+    ASSERT_EQ(mol->natoms, 2);
+    ASSERT_EQ(lmp->atom->natoms, 2);
+    mol->compute_mass();
+    mol->compute_com();
+    EXPECT_NEAR(mol->masstotal, 2.0, EPSILON);
+    EXPECT_NEAR(mol->com[0], 0.0, EPSILON);
+    EXPECT_NEAR(mol->com[1], 0.0, EPSILON);
+    EXPECT_NEAR(mol->com[2], 0.0, EPSILON);
+    EXPECT_EQ(mol->comatom, 1);
+    ASSERT_NE(mol->mu, nullptr);
+    EXPECT_NEAR(mol->mu[0][0], 1.0, EPSILON);
+    EXPECT_NEAR(mol->mu[0][1], 0.0, EPSILON);
+    EXPECT_NEAR(mol->mu[0][2], 0.0, EPSILON);
+    EXPECT_NEAR(mol->mu[1][0], 1.0, EPSILON);
+    EXPECT_NEAR(mol->mu[1][1], 1.0, EPSILON);
+    EXPECT_NEAR(mol->mu[1][2], 0.0, EPSILON);
+    EXPECT_NEAR(mol->maxextent, 2.0, EPSILON);
+    // dipoles should be rotated by 90 degrees clockwise around the z axis
+    double **mu = lmp->atom->mu;
+    ASSERT_NE(mu, nullptr);
+    EXPECT_NEAR(mu[0][0], 0.0, EPSILON);
+    EXPECT_NEAR(mu[0][1], 1.0, EPSILON);
+    EXPECT_NEAR(mu[0][2], 0.0, EPSILON);
+    EXPECT_NEAR(mu[0][3], 1.0, EPSILON);
+    EXPECT_NEAR(mu[1][0], -1.0, EPSILON);
+    EXPECT_NEAR(mu[1][1], 1.0, EPSILON);
+    EXPECT_NEAR(mu[1][2], 0.0, EPSILON);
+    EXPECT_NEAR(mu[1][3], sqrt(2.0), EPSILON);
+}
+
 int main(int argc, char **argv)
 {
     MPI_Init(&argc, &argv);