diff --git a/doc/src/molecule.rst b/doc/src/molecule.rst index b930a9fc65..05821e6076 100644 --- a/doc/src/molecule.rst +++ b/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 -(not identical) to the data file read by the :doc:`read_data ` -commands, and is as follows. +The format of an individual molecule file looks similar but is +different than that of a data file read by the :doc:`read_data ` +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; -it typically contains a description of the file or a comment from the software -that created the file. +first line of the header and thus of the molecule file is *always* +skipped; it typically contains a description of the file or a comment +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 -* Nb *bonds* = # of bonds Nb in molecule, default = 0 -* Na *angles* = # of angles Na in molecule, default = 0 -* 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 + .. list-table:: + :header-rows: 1 + :widths: auto -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. + * - Number(s) + - Keyword + - Meaning + - Default Value + * - N + - atoms + - # 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 ` 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 diff --git a/src/atom.cpp b/src/atom.cpp index b604c54e6b..c08df16614 100644 --- a/src/atom.cpp +++ b/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::vectorAtom::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) diff --git a/src/molecule.cpp b/src/molecule.cpp index 6e2d3891d3..f83d8658df 100644 --- a/src/molecule.cpp +++ b/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), - num_angle(nullptr), angle_type(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), - angle_atom3(nullptr), num_dihedral(nullptr), dihedral_type(nullptr), dihedral_atom1(nullptr), - dihedral_atom2(nullptr), dihedral_atom3(nullptr), dihedral_atom4(nullptr), - num_improper(nullptr), improper_type(nullptr), improper_atom1(nullptr), improper_atom2(nullptr), - improper_atom3(nullptr), improper_atom4(nullptr), nspecial(nullptr), special(nullptr), - shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), avec_body(nullptr), - ibodyparams(nullptr), dbodyparams(nullptr), fragmentmask(nullptr), dx(nullptr), dxcom(nullptr), - dxbody(nullptr), quat_external(nullptr), fp(nullptr), count(nullptr) + radius(nullptr), rmass(nullptr), mu(nullptr), num_bond(nullptr), bond_type(nullptr), + bond_atom(nullptr), num_angle(nullptr), angle_type(nullptr), angle_atom1(nullptr), + angle_atom2(nullptr), angle_atom3(nullptr), num_dihedral(nullptr), dihedral_type(nullptr), + dihedral_atom1(nullptr), dihedral_atom2(nullptr), dihedral_atom3(nullptr), + dihedral_atom4(nullptr), num_improper(nullptr), improper_type(nullptr), improper_atom1(nullptr), + improper_atom2(nullptr), improper_atom3(nullptr), improper_atom4(nullptr), nspecial(nullptr), + special(nullptr), shake_flag(nullptr), shake_atom(nullptr), shake_type(nullptr), + avec_body(nullptr), ibodyparams(nullptr), dbodyparams(nullptr), fragmentmask(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++) diff --git a/src/molecule.h b/src/molecule.h index 06a1211ea3..faba036aab 100644 --- a/src/molecule.h +++ b/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 *); diff --git a/unittest/formats/test_molecule_file.cpp b/unittest/formats/test_molecule_file.cpp index 8fe1fc1eab..c798d2f4c2 100644 --- a/unittest/formats/test_molecule_file.cpp +++ b/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);