diff --git a/gmso/formats/mcf.py b/gmso/formats/mcf.py
index 4558a3038..1b3601703 100644
--- a/gmso/formats/mcf.py
+++ b/gmso/formats/mcf.py
@@ -3,16 +3,21 @@
import warnings
import networkx as nx
-import sympy
+import numpy as np
+import symengine
import unyt as u
from gmso import __version__
from gmso.core.topology import Topology
+from gmso.core.views import PotentialFilters
from gmso.exceptions import GMSOError
from gmso.formats.formats_registry import saves_as
from gmso.lib.potential_templates import PotentialTemplateLibrary
from gmso.utils.compatibility import check_compatibility
-from gmso.utils.conversions import convert_ryckaert_to_opls
+from gmso.utils.conversions import (
+ convert_opls_to_ryckaert,
+ convert_ryckaert_to_fourier,
+)
__all__ = ["write_mcf"]
@@ -24,7 +29,7 @@ def write_mcf(top, filename):
"""Generate a Cassandra MCF from a gmso.core.Topology object.
The MCF file stores the topology information for a single
- species (i.e., compound) in the Cassandra Monte Carlo software
+ species (i.e., molecule) in the Cassandra Monte Carlo software
(https://cassandra.nd.edu). The gmso.Topology object provided to
this function should therefore be for a single molecule with the
relevant forcefield parameters. One MCF file will be required
@@ -34,56 +39,67 @@ def write_mcf(top, filename):
----------
top : gmso.core.Topology
Topology object
- filename : str
- Path of the output file
+ filename : str or list
+ Path of the output file(s). If a string is provided
+ and the Topology object has more than one subtopology,
+ an integer suffix will be appended to the string.
+ If a list is provided and there is more than one subtopology,
+ the names of the output files will be
+ each element in the list. The number of element in the list
+ should match the
+ number of unique subtopologies.
Notes
-----
Atom indexing begins at 1. See
https://cassandra.nd.edu/index.php/documentation for a complete
description of the MCF format.
-
"""
- _check_compatibility(top)
-
- # Identify atoms in rings and Cassandra 'fragments'
- in_ring, frag_list, frag_conn = _id_rings_fragments(top)
-
- # TODO: What oh what to do about subtops?
- # For now refuse topologies with subtops as MCF writer is for
- # single molecules
- if len(top.unique_site_labels("molecule")) > 1:
- raise GMSOError(
- "MCF writer does not support multiple molecules. "
- "Please provide a single molecule as an gmso.Topology "
- "object to the MCF writer."
- )
+ subtops = []
+ for molecule in top.unique_site_labels(name_only=True):
+ subtops.append(top.create_subtop("molecule", (molecule, 1)))
+
+ if len(subtops) > 1:
+ if len(filename) != len(subtops):
+ raise ValueError(
+ "write_mcf: Number of filenames must match number of unique species in the Topology object"
+ )
- # Now we write the MCF file
- with open(filename, "w") as mcf:
- header = (
- "!***************************************"
- "****************************************\n"
- "!Molecular connectivity file\n"
- "!***************************************"
- "****************************************\n"
- "!File {} written by gmso {} at {}\n\n".format(
- filename, __version__, str(datetime.datetime.now())
+ for idx, subtop in enumerate(subtops):
+ _check_compatibility(subtop)
+
+ # Identify atoms in rings and Cassandra 'fragments'
+ in_ring, frag_list, frag_conn = _id_rings_fragments(subtop)
+
+ if len(subtops) > 1:
+ if isinstance(filename, str):
+ filename = filename[:-4] + f"_{idx}.mcf"
+ elif isinstance(filename, list):
+ filename = filename[idx]
+
+ # Now we write the MCF file
+ with open(filename, "w") as mcf:
+ header = (
+ "!***************************************"
+ "****************************************\n"
+ "!Molecular connectivity file\n"
+ "!***************************************"
+ "****************************************\n"
+ f"!File {filename} written by gmso {__version__} "
+ f"at {str(datetime.datetime.now())}\n\n"
)
- )
- mcf.write(header)
- _write_atom_information(mcf, top, in_ring)
- _write_bond_information(mcf, top)
- _write_angle_information(mcf, top)
- _write_dihedral_information(mcf, top)
- # TODO: Add improper information
- # _write_improper_information(mcf, top)
- _write_fragment_information(mcf, top, frag_list, frag_conn)
- _write_intrascaling_information(mcf, top)
+ mcf.write(header)
+ _write_atom_information(mcf, subtop, in_ring)
+ _write_bond_information(mcf, subtop)
+ _write_angle_information(mcf, subtop)
+ _write_dihedral_information(mcf, subtop)
+ _write_improper_information(mcf, subtop)
+ _write_fragment_information(mcf, subtop, frag_list, frag_conn)
+ _write_intrascaling_information(mcf, subtop)
- # That's all, folks!
- mcf.write("\n\nEND\n")
+ # That's all, folks!
+ mcf.write("\n\nEND\n")
def _id_rings_fragments(top):
@@ -102,16 +118,21 @@ def _id_rings_fragments(top):
Atom ids belonging to each fragment
frag_conn : list
Fragment ids of connected fragments
-
"""
# Identify atoms in rings
bond_graph = nx.Graph()
bond_graph.add_edges_from(
- [[bond.atom1.idx, bond.atom2.idx] for bond in top.bonds]
+ [
+ [
+ top.get_index(bond.connection_members[0]),
+ top.get_index(bond.connection_members[1]),
+ ]
+ for bond in top.bonds
+ ]
)
if len(top.bonds) == 0:
warnings.warn(
- "No bonds found. Cassandra will interpet " "this as a rigid species"
+ "No bonds found. Cassandra will interpet this as a rigid species"
)
in_ring = [False] * len(top.sites)
frag_list = []
@@ -119,7 +140,7 @@ def _id_rings_fragments(top):
return in_ring, frag_list, frag_conn
# Check if entire molecule is connected. Warn if not.
- if nx.is_connected(bond_graph) == False:
+ if not nx.is_connected(bond_graph):
raise ValueError(
"Not all components of the molecule are connected. "
"MCF files are for a single molecule and thus "
@@ -143,21 +164,24 @@ def _id_rings_fragments(top):
i: list(bond_graph.neighbors(i))
for i in range(bond_graph.number_of_nodes())
}
- # First ID fused rings
- fused_rings = []
- rings_to_remove = []
- for i in range(len(all_rings)):
- ring1 = all_rings[i]
- for j in range(i + 1, len(all_rings)):
- ring2 = all_rings[j]
- shared_atoms = list(set(ring1) & set(ring2))
- if len(shared_atoms) == 2:
- fused_rings.append(list(set(ring1 + ring2)))
- rings_to_remove.append(ring1)
- rings_to_remove.append(ring2)
- for ring in rings_to_remove:
- all_rings.remove(ring)
- all_rings = all_rings + fused_rings
+
+ # Handle fused/adjoining rings
+ rings_changed = True
+ while rings_changed:
+ rings_changed = False
+ for ring1 in all_rings:
+ if rings_changed:
+ break
+ for ring2 in all_rings:
+ if ring1 == ring2:
+ continue
+ if len(set(ring1) & set(ring2)) > 0:
+ all_rings.remove(ring1)
+ all_rings.remove(ring2)
+ all_rings.append(list(set(ring1 + ring2)))
+ rings_changed = True
+ break
+
# ID fragments which contain a ring
for ring in all_rings:
adjacent_atoms = []
@@ -209,53 +233,66 @@ def _write_atom_information(mcf, top, in_ring):
Topology object
in_ring : list
Boolean for each atom idx True if atom belongs to a ring
-
"""
- names = [site.name for site in top.sites]
- types = [site.atom_type.name for site in top.sites]
+ # Based upon Cassandra; updated following 1.2.2 release
+ max_element_length = 6
+ max_atomtype_length = 20
+
+ sites, names, atypes_list = zip(
+ *[(site, site.name, site.atom_type.name) for site in top.sites]
+ )
# Check constraints on atom type length and element name length
- # TODO: Update these following Cassandra release
- # to be more reasonable values
n_unique_names = len(set(names))
for name in names:
- if len(name) > 2:
+ if len(name) > max_element_length:
warnings.warn(
- "Warning, name {} will be shortened "
- "to two characters. Please confirm your final "
- "MCF.".format(name)
+ f"Name: {name} will be shortened to {max_element_length}"
+ "characters. Please confirm your final MCF."
)
# Confirm that shortening names to two characters does not
# cause two previously unique atom names to become identical.
- names = [name[:2] for name in names]
+ names = [name[:max_element_length] for name in names]
if len(set(names)) < n_unique_names:
warnings.warn(
- "Warning, the number of unique names has been "
- "reduced due to shortening the name to two "
- "characters."
+ "The number of unique names has been reduced due"
+ f"to shortening the name to {max_element_length} characters."
)
- n_unique_types = len(set(types))
- for itype in types:
- if len(itype) > 6:
+ pfilter = PotentialFilters.UNIQUE_SORTED_NAMES
+ n_unique_types = top.atom_types(filter_by=pfilter)
+
+ for type_ in n_unique_types:
+ if len(type_.name) > max_atomtype_length:
warnings.warn(
- "Warning, type name {} will be shortened to six "
- "characters as {}. Please confirm your final "
- "MCF.".format(itype, itype[-6:])
+ f"Type name: {type_.name} will be shortened to "
+ f"{max_atomtype_length} characters as "
+ f"{type[-max_atomtype_length:]}. Please confirm your final MCF."
)
- types = [itype[-6:] for itype in types]
- if len(set(types)) < n_unique_types:
+ atypes_list = [itype[-max_atomtype_length:] for itype in atypes_list]
+
+ if len(set(atypes_list)) < len(n_unique_types):
warnings.warn(
- "Warning, the number of unique atomtypes has been "
- "reduced due to shortening the atomtype name to six "
- "characters."
+ "The number of unique atomtypes has been reduced due to "
+ f"shortening the atomtype name to {max_atomtype_length} characters."
+ )
+
+ # Check charge neutrality
+ net_q = 0.0
+ for idx, site in enumerate(range(top.n_sites)):
+ net_q += top.sites[idx].charge.in_units(u.elementary_charge).value
+
+ if not np.isclose(net_q, 0.0):
+ raise ValueError(
+ "Net charge of the system is not zero. "
+ "Cassandra MFC requires a neutral system."
)
# Detect VDW style
vdw_styles = set()
- for site in top.sites:
- vdw_styles.add(_get_vdw_style(site.atom_type))
+ for site in sites:
+ vdw_styles.add(_get_vdw_style(site))
if len(vdw_styles) > 1:
raise GMSOError(
"More than one vdw_style detected. "
@@ -283,12 +320,12 @@ def _write_atom_information(mcf, top, in_ring):
"{:<4d} "
"{:<6s} "
"{:<2s} "
- "{:7.3f} "
- "{:7.3f} ".format(
+ "{:8.4f} "
+ "{:12.8f} ".format(
idx + 1,
- types[idx],
+ atypes_list[idx],
names[idx],
- site.mass.in_units(u.amu).value,
+ site.atom_type.mass.in_units(u.amu).value,
site.charge.in_units(u.elementary_charge).value,
)
)
@@ -354,8 +391,8 @@ def _write_bond_information(mcf, top):
"{:s} "
"{:10.5f}\n".format(
idx + 1,
- bond.connection_members[0].idx + 1, # TODO: Confirm the +1 here
- bond.connection_members[1].idx + 1,
+ top.get_index(bond.connection_members[0]) + 1,
+ top.get_index(bond.connection_members[1]) + 1,
"fixed",
bond.connection_type.parameters["r_eq"]
.in_units(u.Angstrom)
@@ -374,11 +411,10 @@ def _write_angle_information(mcf, top):
top : Topology
Topology object
"""
- # TODO: Add support for fixed angles
- angle_style = "harmonic"
header = (
"\n!Angle Format\n"
"!index i j k type parameters\n"
+ '!type="fixed", parms=equilibrium_angle\n'
'!type="harmonic", parms=force_constant equilibrium_angle\n'
"\n# Angle_Info\n"
)
@@ -387,27 +423,38 @@ def _write_angle_information(mcf, top):
mcf.write("{:d}\n".format(len(top.angles)))
for idx, angle in enumerate(top.angles):
mcf.write(
- "{:<4d} "
- "{:<4d} "
- "{:<4d} "
- "{:<4d} "
- "{:s} "
- "{:10.5f} "
- "{:10.5f}\n".format(
- idx + 1,
- angle.connection_members[0].idx + 1,
- angle.connection_members[1].idx
- + 1, # TODO: Confirm order for angles i-j-k
- angle.connection_members[2].idx + 1,
- angle_style,
- (0.5 * angle.connection_type.parameters["k"] / u.kb)
- .in_units("K/rad**2")
- .value, # TODO: k vs. k/2. conversion
- angle.connection_type.parameters["theta_eq"]
- .in_units(u.degree)
- .value,
- )
+ f"{idx + 1:<4d} "
+ f"{top.get_index(angle.connection_members[0]) + 1:<4d} "
+ f"{top.get_index(angle.connection_members[1]) + 1:<4d} "
+ f"{top.get_index(angle.connection_members[2]) + 1:<4d} "
)
+ angle_style = _get_angle_style(angle)
+ if angle_style == "fixed":
+ mcf.write(
+ "{:s} "
+ "{:10.5f}\n".format(
+ angle_style,
+ angle.connection_type.parameters["theta_eq"]
+ .in_units(u.degree)
+ .value,
+ )
+ )
+ elif angle_style == "harmonic":
+ mcf.write(
+ "{:s} "
+ "{:10.5f} "
+ "{:10.5f}\n".format(
+ angle_style,
+ (0.5 * angle.connection_type.parameters["k"] / u.kb)
+ .in_units("K/rad**2")
+ .value, # TODO: k vs. k/2. conversion
+ angle.connection_type.parameters["theta_eq"]
+ .in_units(u.degree)
+ .value,
+ )
+ )
+ else:
+ raise GMSOError("Unsupported angle style for Cassandra MCF writer")
def _write_dihedral_information(mcf, top):
@@ -419,8 +466,6 @@ def _write_dihedral_information(mcf, top):
The file object of the Cassandra mcf being written
top : Topology
Topology object
- dihedral_style : string
- Dihedral style for Cassandra to use
"""
# Dihedral info
header = (
@@ -435,7 +480,6 @@ def _write_dihedral_information(mcf, top):
mcf.write(header)
- # TODO: Are impropers buried in dihedrals?
mcf.write("{:d}\n".format(len(top.dihedrals)))
for idx, dihedral in enumerate(top.dihedrals):
mcf.write(
@@ -445,20 +489,39 @@ def _write_dihedral_information(mcf, top):
"{:<4d} "
"{:<4d} ".format(
idx + 1,
- dihedral.connection_members[0].idx + 1,
- dihedral.connection_members[1].idx + 1,
- dihedral.connection_members[2].idx + 1,
- dihedral.connection_members[3].idx + 1,
+ top.get_index(dihedral.connection_members[0]) + 1,
+ top.get_index(dihedral.connection_members[1]) + 1,
+ top.get_index(dihedral.connection_members[2]) + 1,
+ top.get_index(dihedral.connection_members[3]) + 1,
)
)
dihedral_style = _get_dihedral_style(dihedral)
+ dihedral_style = dihedral_style.upper()
# If ryckaert, convert to opls
- if dihedral_style == "ryckaert":
- dihedral.connection_type = convert_ryckaert_to_opls(
+ if dihedral_style == "RYCKAERT":
+ dihedral.connection_type = convert_ryckaert_to_fourier(
dihedral.connection_type
)
- dihedral_style = "opls"
- if dihedral_style == "opls":
+ # The GMSO Fourier style is named OPLS in Cassandra. See
+ # https://cassandra-mc.readthedocs.io/en/latest/guides/forcefield.html?highlight=opls#dihedrals
+
+ dihedral_style = "FOURIER"
+
+ if dihedral_style == "OPLS":
+ dihedral.connection_type = convert_opls_to_ryckaert(
+ dihedral.connection_type
+ )
+ dihedral.connection_type = convert_ryckaert_to_fourier(
+ dihedral.connection_type
+ )
+ dihedral_style = "FOURIER"
+
+ if dihedral_style == "FOURIER":
+ # Cassandra supports the OPLS (GMSO's Fourier) functional form up 4 terms, namely
+ # E = a0 + a1 * (1 + cos(phi)) + a2 * (1 - cos(2*phi)) + a3 * (1 + cos(3*phi))
+ # The GMSO Fourier potential has terms up to 0.5 * k4 * ( 1 - cos ( 4 * phi))
+ # So we need to exclude the last term in the GMSO topology.
+ dihedral_style = "OPLS"
mcf.write(
"{:s} "
"{:10.5f} "
@@ -484,7 +547,7 @@ def _write_dihedral_information(mcf, top):
.value,
)
)
- elif dihedral_style == "charmm":
+ elif dihedral_style == "CHARMM":
mcf.write(
"{:s} "
"{:10.5f} "
@@ -500,12 +563,12 @@ def _write_dihedral_information(mcf, top):
.value,
)
)
- elif dihedral_style == "harmonic":
+ elif dihedral_style == "HARMONIC":
mcf.write(
"{:s} "
"{:10.5f} "
"{:10.5f}\n".format(
- dihedral_style,
+ dihedral_style.lower(),
0.5
* dihedral.connection_type.parameters["k"]
.in_units("kJ/mol")
@@ -542,19 +605,19 @@ def _write_improper_information(mcf, top):
mcf.write(header)
mcf.write("{:d}\n".format(len(top.impropers)))
- improper_type = "harmonic"
+ improper_style = "harmonic"
for i, improper in enumerate(top.impropers):
mcf.write(
"{:<4d} {:<4d} {:<4d} {:<4d} {:<4d}"
- " {:s} {:8.3f} {:8.3f}\n".format(
+ " {:s} {:10.5f} {:10.5f}\n".format(
i + 1,
- improper.atom1.idx + 1,
- improper.atom2.idx + 1,
- improper.atom3.idx + 1,
- improper.atom4.idx + 1,
- improper_type,
- improper.type.psi_k * KCAL_TO_KJ,
- improper.type.psi_eq,
+ top.get_index(improper.connection_members[0]) + 1,
+ top.get_index(improper.connection_members[1]) + 1,
+ top.get_index(improper.connection_members[2]) + 1,
+ top.get_index(improper.connection_members[3]) + 1,
+ improper_style,
+ 0.5 * improper.connection_type.parameters["k"],
+ improper.connection_type.parameters["phi_eq"],
)
)
@@ -642,47 +705,62 @@ def _check_compatibility(top):
"""Check Topology object for compatibility with Cassandra MCF format."""
if not isinstance(top, Topology):
raise GMSOError("MCF writer requires a Topology object.")
- if not all([site.atom_type.name for site in top.sites]):
+ if not all([site.atom_type for site in top.sites]):
raise GMSOError(
"MCF writing not supported without parameterized forcefield."
)
accepted_potentials = (
potential_templates["LennardJonesPotential"],
potential_templates["MiePotential"],
+ potential_templates["FixedBondPotential"],
+ potential_templates["HarmonicBondPotential"],
potential_templates["HarmonicAnglePotential"],
+ potential_templates["FixedAnglePotential"],
potential_templates["PeriodicTorsionPotential"],
potential_templates["OPLSTorsionPotential"],
+ potential_templates["FourierTorsionPotential"],
potential_templates["RyckaertBellemansTorsionPotential"],
)
check_compatibility(top, accepted_potentials)
-def _get_vdw_style(atom_type):
+def _get_vdw_style(site):
"""Return the vdw style."""
vdw_styles = {
"LJ": potential_templates["LennardJonesPotential"],
"Mie": potential_templates["MiePotential"],
}
- return _get_potential_style(vdw_styles, atom_type)
+ return _get_potential_style(vdw_styles, site.atom_type)
+
+
+def _get_angle_style(angle):
+ """Return the angle style."""
+ angle_styles = {
+ "harmonic": potential_templates["HarmonicAnglePotential"],
+ "fixed": potential_templates["FixedAnglePotential"],
+ }
+
+ return _get_potential_style(angle_styles, angle.connection_type)
-def _get_dihedral_style(dihedral_type):
+def _get_dihedral_style(dihedral):
"""Return the dihedral style."""
dihedral_styles = {
"charmm": potential_templates["PeriodicTorsionPotential"],
"harmonic": potential_templates["HarmonicTorsionPotential"],
"opls": potential_templates["OPLSTorsionPotential"],
+ "fourier": potential_templates["FourierTorsionPotential"],
"ryckaert": potential_templates["RyckaertBellemansTorsionPotential"],
}
- return _get_potential_style(dihedral_styles, dihedral_type)
+ return _get_potential_style(dihedral_styles, dihedral.connection_type)
def _get_potential_style(styles, potential):
"""Return the potential style."""
for style, ref in styles.items():
if ref.independent_variables == potential.independent_variables:
- if sympy.simplify(ref.expression - potential.expression) == 0:
+ if symengine.expand(ref.expression - potential.expression) == 0:
return style
return False
diff --git a/gmso/lib/jsons/FixedAnglePotential.json b/gmso/lib/jsons/FixedAnglePotential.json
new file mode 100644
index 000000000..c32eef5ff
--- /dev/null
+++ b/gmso/lib/jsons/FixedAnglePotential.json
@@ -0,0 +1,8 @@
+{
+ "name": "FixedAnglePotential",
+ "expression": "DiracDelta(theta-theta_eq)",
+ "independent_variables": "theta",
+ "expected_parameters_dimensions": {
+ "theta_eq": "angle"
+ }
+}
diff --git a/gmso/lib/jsons/FixedBondPotential.json b/gmso/lib/jsons/FixedBondPotential.json
new file mode 100644
index 000000000..f4845f402
--- /dev/null
+++ b/gmso/lib/jsons/FixedBondPotential.json
@@ -0,0 +1,8 @@
+{
+ "name": "FixedBondPotential",
+ "expression": "DiracDelta(r-r_eq)",
+ "independent_variables": "r",
+ "expected_parameters_dimensions": {
+ "r_eq": "length"
+ }
+}
diff --git a/gmso/tests/base_test.py b/gmso/tests/base_test.py
index b5dc15395..a91d410e4 100644
--- a/gmso/tests/base_test.py
+++ b/gmso/tests/base_test.py
@@ -263,6 +263,14 @@ def typed_water_system(self, water_system):
top = apply(top, ff)
return top
+ @pytest.fixture
+ def typed_tip3p_rigid_system(self, water_system):
+ top = water_system
+ top.identify_connections()
+ ff = ForceField(get_path("tip3p-rigid.xml"))
+ top = apply(top, ff)
+ return top
+
@pytest.fixture
def foyer_fullerene(self):
from foyer.tests.utils import get_fn
diff --git a/gmso/tests/files/ethane-rigid.xml b/gmso/tests/files/ethane-rigid.xml
new file mode 100644
index 000000000..cb857f7ef
--- /dev/null
+++ b/gmso/tests/files/ethane-rigid.xml
@@ -0,0 +1,23 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/gmso/tests/files/ethylene.xml b/gmso/tests/files/ethylene.xml
index 6345773dc..7059a1d48 100644
--- a/gmso/tests/files/ethylene.xml
+++ b/gmso/tests/files/ethylene.xml
@@ -51,21 +51,21 @@
-
-
-
-
-
-
-
+
+
+
+
+
+
+
-
-
-
-
-
-
+
+
+
+
+
+
diff --git a/gmso/tests/files/tip3p-rigid.xml b/gmso/tests/files/tip3p-rigid.xml
new file mode 100644
index 000000000..956424a50
--- /dev/null
+++ b/gmso/tests/files/tip3p-rigid.xml
@@ -0,0 +1,38 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/gmso/tests/test_mcf.py b/gmso/tests/test_mcf.py
index bfe1279b4..37b255917 100644
--- a/gmso/tests/test_mcf.py
+++ b/gmso/tests/test_mcf.py
@@ -1,130 +1,383 @@
+import re
+import subprocess
+
+import mbuild as mb
import numpy as np
import pytest
import unyt as u
+from gmso.core.forcefield import ForceField
from gmso.exceptions import EngineIncompatibilityError
+from gmso.external import from_mbuild
from gmso.formats.mcf import write_mcf
+from gmso.parameterization import apply
from gmso.tests.base_test import BaseTest
+from gmso.tests.utils import get_path
+from gmso.utils.conversions import convert_ryckaert_to_fourier
+from gmso.utils.io import get_fn, has_cassandra, has_parmed, import_
+if has_cassandra:
+ mc = import_("mosdef_cassandra")
-class TestMCF(BaseTest):
- def test_write_lj_simple(self, n_typed_ar_system):
- top = n_typed_ar_system(n_sites=1)
- top.save("ar.mcf")
+if has_parmed:
+ pmd = import_("parmed")
- def test_write_mie_simple(self, n_typed_xe_mie):
- top = n_typed_xe_mie()
- top.save("xe.mcf")
+def parse_mcf(filename):
+ mcf_data = []
+ mcf_idx = {}
+ with open(filename) as f:
+ for line in f:
+ mcf_data.append(line.strip().split())
+
+ for idx, line in enumerate(mcf_data):
+ if len(line) > 1:
+ if line[1] == "Atom_Info":
+ mcf_idx["Atom_Info"] = idx
+ if len(line) > 1:
+ if line[1] == "Bond_Info":
+ mcf_idx["Bond_Info"] = idx
+ if len(line) > 1:
+ if line[1] == "Angle_Info":
+ mcf_idx["Angle_Info"] = idx
+ if len(line) > 1:
+ if line[1] == "Dihedral_Info":
+ mcf_idx["Dihedral_Info"] = idx
+ if len(line) > 1:
+ if line[1] == "Fragment_Info":
+ mcf_idx["Fragment_Info"] = idx
+ if len(line) > 1:
+ if line[1] == "Fragment_Connectivity":
+ mcf_idx["Fragment_Connectivity"] = idx
+ if len(line) > 1:
+ if line[1] == "Intra_Scaling":
+ mcf_idx["Intra_Scaling"] = idx
+
+ return mcf_data, mcf_idx
+
+
+def run_cassandra(cassandra, inp_file):
+ """Calls Cassandra. Taken from mosdef_cassandra v0.3.2"""
+ cassandra_cmd = "{cassandra} {inp_file}".format(
+ cassandra=cassandra, inp_file=inp_file
+ )
+ p = subprocess.Popen(
+ cassandra_cmd,
+ shell=True,
+ stdout=subprocess.PIPE,
+ stderr=subprocess.PIPE,
+ universal_newlines=True,
+ )
+ out, err = p.communicate()
+
+ if p.returncode != 0 or "error" in err.lower() or "error" in out.lower():
+ return 1, out, err
+ return 0, out, err
+
+
+def is_charge_neutral(mcf_data, mcf_idx):
+ n_sites = int(mcf_data[mcf_idx["Atom_Info"] + 1][0])
+ net_q = 0.0
+ for idx, site in enumerate(range(n_sites)):
+ net_q += float(mcf_data[mcf_idx["Atom_Info"] + 2 + idx][4])
+ return np.isclose(
+ net_q,
+ 0.0,
+ )
+
+
+class TestMCF(BaseTest):
def test_write_lj_full(self, n_typed_ar_system):
top = n_typed_ar_system(n_sites=1)
top.save("ar.mcf")
- mcf_data = []
- with open("ar.mcf") as f:
- for line in f:
- mcf_data.append(line.strip().split())
+ mcf_data, mcf_idx = parse_mcf("ar.mcf")
- for idx, line in enumerate(mcf_data):
- if len(line) > 1:
- if line[1] == "Atom_Info":
- atom_section_start = idx
+ assert is_charge_neutral(mcf_data, mcf_idx)
- assert mcf_data[atom_section_start + 1][0] == "1"
- assert mcf_data[atom_section_start + 2][1] == "Ar"
- assert mcf_data[atom_section_start + 2][2] == "Ar"
- assert mcf_data[atom_section_start + 2][5] == "LJ"
+ assert mcf_data[mcf_idx["Atom_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][1] == "Ar"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][2] == "Ar"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][5] == "LJ"
assert np.isclose(
- float(mcf_data[atom_section_start + 2][3]),
- top.sites[0].mass.in_units(u.amu).value,
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][3]),
+ top.sites[0].atom_type.mass.in_units(u.amu).value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][4]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][4]),
top.sites[0].charge.in_units(u.elementary_charge).value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][6]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][6]),
(top.sites[0].atom_type.parameters["epsilon"] / u.kb)
.in_units(u.K)
.value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][7]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][7]),
top.sites[0]
.atom_type.parameters["sigma"]
.in_units(u.Angstrom)
.value,
)
+ assert mcf_data[mcf_idx["Bond_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Dihedral_Info"] + 1][0] == "0"
+
+ assert mcf_data[mcf_idx["Fragment_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2] == ["1", "1", "1"]
+
+ assert mcf_data[mcf_idx["Fragment_Connectivity"] + 1][0] == "0"
+
+ assert np.allclose(float(mcf_data[-5][0]), 0.0)
+ assert np.allclose(float(mcf_data[-5][1]), 0.0)
+ assert np.allclose(float(mcf_data[-5][2]), 0.5)
+ assert np.allclose(float(mcf_data[-5][3]), 1.0)
+ assert np.allclose(float(mcf_data[-4][0]), 0.0)
+ assert np.allclose(float(mcf_data[-4][1]), 0.0)
+ assert np.allclose(float(mcf_data[-4][2]), 0.5)
+ assert np.allclose(float(mcf_data[-4][3]), 1.0)
+
+ def test_write_not_neutral(self, n_typed_ar_system):
+ top = n_typed_ar_system(n_sites=1)
+ top.sites[0].charge = 1.0 * u.elementary_charge
+ with pytest.raises(ValueError):
+ top.save("ar.mcf")
def test_write_mie_full(self, n_typed_xe_mie):
top = n_typed_xe_mie()
top.save("xe.mcf")
- mcf_data = []
- with open("xe.mcf") as f:
- for line in f:
- mcf_data.append(line.strip().split())
-
- for idx, line in enumerate(mcf_data):
- if len(line) > 1:
- if line[1] == "Atom_Info":
- atom_section_start = idx
+ mcf_data, mcf_idx = parse_mcf("xe.mcf")
+ assert is_charge_neutral(mcf_data, mcf_idx)
# Check a some atom info
- assert mcf_data[atom_section_start + 1][0] == "1"
- assert mcf_data[atom_section_start + 2][1] == "Xe"
- assert mcf_data[atom_section_start + 2][2] == "Xe"
- assert mcf_data[atom_section_start + 2][5] == "Mie"
+ assert mcf_data[mcf_idx["Atom_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][1] == "Xe"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][2] == "Xe"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][5] == "Mie"
assert np.isclose(
- float(mcf_data[atom_section_start + 2][3]),
- top.sites[0].mass.in_units(u.amu).value,
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][3]),
+ top.sites[0].atom_type.mass.in_units(u.amu).value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][4]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][4]),
top.sites[0].charge.in_units(u.elementary_charge).value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][6]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][6]),
(top.sites[0].atom_type.parameters["epsilon"] / u.kb)
.in_units(u.K)
.value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][7]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][7]),
top.sites[0]
.atom_type.parameters["sigma"]
.in_units(u.Angstrom)
.value,
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][8]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][8]),
top.sites[0].atom_type.parameters["n"],
)
assert np.isclose(
- float(mcf_data[atom_section_start + 2][9]),
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][9]),
top.sites[0].atom_type.parameters["m"],
)
- def test_modified_potentials(self, n_typed_ar_system):
- top = n_typed_ar_system(n_sites=1)
+ assert mcf_data[mcf_idx["Bond_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Dihedral_Info"] + 1][0] == "0"
+
+ assert mcf_data[mcf_idx["Fragment_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2] == ["1", "1", "1"]
+
+ assert mcf_data[mcf_idx["Fragment_Connectivity"] + 1][0] == "0"
+
+ assert np.allclose(float(mcf_data[-5][0]), 0.0)
+ assert np.allclose(float(mcf_data[-5][1]), 0.0)
+ assert np.allclose(float(mcf_data[-5][2]), 0.5)
+ assert np.allclose(float(mcf_data[-5][3]), 1.0)
+ assert np.allclose(float(mcf_data[-4][0]), 0.0)
+ assert np.allclose(float(mcf_data[-4][1]), 0.0)
+ assert np.allclose(float(mcf_data[-4][2]), 0.5)
+ assert np.allclose(float(mcf_data[-4][3]), 1.0)
+
+ def test_write_single_fragment_two_atoms(self):
+ """
+ The main purpose of his test is to check for valid
+ fragment information ouput for molecules that have
+ no angles.
+ """
+ ethane = mb.load(get_fn("ethane_ua.mol2"))
+ top = from_mbuild(ethane)
+ ff = ForceField(get_path("ethane-rigid.xml"))
+ top.identify_connections()
+ apply(top, ff, remove_untyped=True)
+ write_mcf(top, "ethane-rigid.mcf")
+
+ mcf_data, mcf_idx = parse_mcf("ethane-rigid.mcf")
+ assert is_charge_neutral(mcf_data, mcf_idx)
+
+ # Assert number of fragments
+ assert mcf_data[mcf_idx["Fragment_Info"] + 1][0] == "1"
+ # Assert number of atoms in the first fragment
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2][1] == "2"
+ # Assert atom IDs in the first fragment
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2][2] == "1"
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2][3] == "2"
+ assert mcf_data[mcf_idx["Fragment_Connectivity"] + 1][0] == "0"
+
+ def test_modified_incompatible_expressions(self, typed_ethane):
+ top = typed_ethane
+
+ # Test that we can't write a MCF with a modified potential
next(iter(top.atom_types)).set_expression("sigma + epsilon*r")
+ with pytest.raises(EngineIncompatibilityError):
+ top.save("lj.mcf")
+ next(iter(top.bond_types)).set_expression("k*(r-r_eq)**3")
with pytest.raises(EngineIncompatibilityError):
- top.save("out.mcf")
+ top.save("bond.mcf")
- alternate_lj = "4*epsilon*sigma**12/r**12 - 4*epsilon*sigma**6/r**6"
- next(iter(top.atom_types)).set_expression(alternate_lj)
+ # Modified angles
+ next(iter(top.angle_types)).set_expression(
+ "0.5 * k*(theta-theta_eq)**3"
+ )
+ with pytest.raises(EngineIncompatibilityError):
+ top.save("angle.mcf")
- top.save("ar.mcf")
+ # Modified dihedrals
+ # next(iter(top.dihedral_types)).set_expression("c0 * c1 * c2 * c3 * c4 * c5")
+ # with pytest.raises(EngineIncompatibilityError):
+ # top.save("dihedral.mcf")
- def test_scaling_factors(self, n_typed_ar_system):
- top = n_typed_ar_system(n_sites=1)
- top.save("ar.mcf")
- mcf_data = []
- with open("ar.mcf") as f:
- for line in f:
- mcf_data.append(line.strip().split())
+ def test_typed_ethylene(self):
+ ethylene = mb.load("C=C", smiles=True)
+ top = from_mbuild(ethylene)
+ ff = ForceField(get_path("ethylene.xml"))
+ top.identify_connections()
+ apply(top, ff, remove_untyped=True)
+ write_mcf(top, "ethylene.mcf")
+
+ mcf_data, mcf_idx = parse_mcf("ethylene.mcf")
+
+ assert is_charge_neutral(mcf_data, mcf_idx)
+
+ # Check atom info
+ assert mcf_data[mcf_idx["Atom_Info"] + 1][0] == "6"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][1] == "opls_143"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][2] == "C"
+ assert mcf_data[mcf_idx["Atom_Info"] + 2][5] == "LJ"
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][3]),
+ top.sites[0].atom_type.mass.in_units(u.amu).value,
+ )
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][4]),
+ top.sites[0].charge.in_units(u.elementary_charge).value,
+ )
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][6]),
+ (top.sites[0].atom_type.parameters["epsilon"] / u.kb)
+ .in_units(u.K)
+ .value,
+ )
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Atom_Info"] + 2][7]),
+ top.sites[0]
+ .atom_type.parameters["sigma"]
+ .in_units(u.Angstrom)
+ .value,
+ )
+
+ # Check bond section
+ assert mcf_data[mcf_idx["Bond_Info"] + 1][0] == "5"
+ assert mcf_data[mcf_idx["Bond_Info"] + 2][3] == "fixed"
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Bond_Info"] + 2][4]),
+ top.bonds[0]
+ .bond_type.parameters["r_eq"]
+ .in_units(u.Angstrom)
+ .value,
+ )
+
+ # Check angle section
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "6"
+ assert mcf_data[mcf_idx["Angle_Info"] + 2][4] == "harmonic"
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Angle_Info"] + 2][6]),
+ top.angles[0]
+ .angle_type.parameters["theta_eq"]
+ .in_units(u.degree)
+ .value,
+ )
+
+ assert np.isclose(
+ 2.0 * float(mcf_data[mcf_idx["Angle_Info"] + 2][5]),
+ (top.angles[0].angle_type.parameters["k"] / u.kb)
+ .in_units(u.K / u.radian**2)
+ .value,
+ )
+
+ # Check dihedral section
+ assert mcf_data[mcf_idx["Dihedral_Info"] + 1][0] == "4"
+ dihedral_style = mcf_data[mcf_idx["Dihedral_Info"] + 2][5].lower()
+ assert dihedral_style.lower() == "opls"
+
+ k = list(ff.dihedral_types.keys())
+ dihedral_type = ff.dihedral_types[k[0]]
+ ff_coeffs = convert_ryckaert_to_fourier(
+ dihedral_type
+ ).parameters.values()
+
+ # We need to drop the last coefficient in the fourier GMSO dihedral type
+ # (the term that has cos(4*phi) since it is not supported in the equivalent OPLS
+ # dihedral type in Cassandra)
+
+ ff_coeffs = np.array([float(x) for x in ff_coeffs])[:-1]
+ mcf_coeffs = np.array(
+ [float(x) for x in mcf_data[mcf_idx["Dihedral_Info"] + 2][6:]]
+ )
+
+ assert np.allclose(ff_coeffs, 2.0 * mcf_coeffs)
+
+ def test_fixed_angles(self, typed_tip3p_rigid_system):
+ top = typed_tip3p_rigid_system
+ write_mcf(top, "tip3p-rigid.mcf")
+
+ mcf_data, mcf_idx = parse_mcf("tip3p-rigid.mcf")
+
+ assert is_charge_neutral(mcf_data, mcf_idx)
+
+ # Check angle section
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Angle_Info"] + 2][4] == "fixed"
+ assert np.isclose(
+ float(mcf_data[mcf_idx["Angle_Info"] + 2][5]),
+ top.angles[0]
+ .angle_type.parameters["theta_eq"]
+ .in_units(u.degree)
+ .value,
+ )
+
+ def test_top_with_n_ar_system(self, n_typed_ar_system):
+ top = n_typed_ar_system(n_sites=10)
+ write_mcf(top, "top-10ar.mcf")
+
+ mcf_data, mcf_idx = parse_mcf("top-10ar.mcf")
+
+ assert is_charge_neutral(mcf_data, mcf_idx)
+
+ assert mcf_data[mcf_idx["Atom_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Bond_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Dihedral_Info"] + 1][0] == "0"
+ assert mcf_data[mcf_idx["Fragment_Info"] + 1][0] == "1"
+ assert mcf_data[mcf_idx["Fragment_Info"] + 2] == ["1", "1", "1"]
+ assert mcf_data[mcf_idx["Fragment_Connectivity"] + 1][0] == "0"
assert np.allclose(float(mcf_data[-5][0]), 0.0)
assert np.allclose(float(mcf_data[-5][1]), 0.0)
assert np.allclose(float(mcf_data[-5][2]), 0.5)
@@ -133,19 +386,181 @@ def test_scaling_factors(self, n_typed_ar_system):
assert np.allclose(float(mcf_data[-4][1]), 0.0)
assert np.allclose(float(mcf_data[-4][2]), 0.5)
assert np.allclose(float(mcf_data[-4][3]), 1.0)
- top.set_lj_scale([0.1, 0.2, 0.5])
- top.set_electrostatics_scale([0.2, 0.4, 0.6])
-
- top.save("ar.mcf", overwrite=True)
- mcf_data = []
- with open("ar.mcf") as f:
- for line in f:
- mcf_data.append(line.strip().split())
- assert np.allclose(float(mcf_data[-5][0]), 0.1)
- assert np.allclose(float(mcf_data[-5][1]), 0.2)
- assert np.allclose(float(mcf_data[-5][2]), 0.5)
- assert np.allclose(float(mcf_data[-5][3]), 1.0)
- assert np.allclose(float(mcf_data[-4][0]), 0.2)
- assert np.allclose(float(mcf_data[-4][1]), 0.4)
- assert np.allclose(float(mcf_data[-4][2]), 0.6)
- assert np.allclose(float(mcf_data[-4][3]), 1.0)
+
+ def test_top_with_mixture(self):
+ pass
+
+ @pytest.mark.skipif(not has_cassandra, reason="cassandra is not installed")
+ def test_in_cassandra(self, typed_ethane):
+ """
+ This test runs a single point energy calculation in Cassandra using an MCF
+ generated by gmso and compare the total energy
+ to the energy of a simulation run with a MCF file generated using mosdef_cassandra
+ (which involves using a parmed.Structure)
+ """
+ from mosdef_cassandra.utils.detect import detect_cassandra_binaries
+ from mosdef_cassandra.writers.writers import write_input
+
+ from gmso.external.convert_parmed import to_parmed
+
+ # First run the mosdef_cassandra simulation. Mosdef_cassandra generates an input file
+ # as well as an MCF. Later, we will use a mosdef_cassandra
+ # input file as a template and replace the MCF file with the GMSO MCF file
+
+ box = mb.Box([3.0, 3.0, 3.0])
+ species = to_parmed(typed_ethane)
+ system = mc.System([box], [species], mols_to_add=[[5]])
+ ensemble = "nvt"
+ moveset = mc.MoveSet(ensemble, [species])
+ mc.run(
+ system=system,
+ moveset=moveset,
+ run_type="equilibration",
+ run_length=0,
+ temperature=300.0 * u.K,
+ run_name="nvt_mbuild",
+ seeds=[12356, 64321],
+ )
+
+ py, fraglib_setup, cassandra = detect_cassandra_binaries()
+
+ # TODO: not sure why the cassandra MCF writer of mBuild
+ # outputs a different intramolecular exclusions relative
+ # to the GMSO writer. This is a temporary fix to make the
+ # test pass. We should investigate this further.
+ # Also, try to use the function top.set_lj_scale() and see how
+ # to update subtopologies.
+
+ write_mcf(typed_ethane, "gmso.mcf")
+ mcf_data, mcf_idx = parse_mcf("gmso.mcf")
+ mcf_data[mcf_idx["Intra_Scaling"] + 1][0:4] = [
+ "0.0",
+ "0.0",
+ "0.0",
+ "1.0",
+ ]
+ mcf_data[mcf_idx["Intra_Scaling"] + 2][0:4] = [
+ "0.0",
+ "0.0",
+ "0.0",
+ "1.0",
+ ]
+ with open("gmso.mcf", mode="w") as f:
+ for line in mcf_data:
+ f.write(" ".join(line) + "\n")
+
+ inp_file = write_input(
+ system=system,
+ moveset=moveset,
+ run_type="equilibration",
+ run_length=0,
+ temperature=300.0 * u.K,
+ run_name="nvt_gmso",
+ seeds=[12356, 64321],
+ )
+
+ with open(inp_file, mode="r") as f:
+ lines = f.read()
+ lines = lines.replace("species1.mcf", "gmso.mcf")
+ # The fragment files section is empty unless
+ # restart option is used in mosdef_cassandra.
+ # See the mosdef_cassandra.writers.inp_functions.generate_input
+ lines = lines.replace(
+ "# Fragment_Files",
+ "# Fragment_Files\nspecies1/frag1/frag1.dat 1\nspecies1/frag2/frag2.dat 2\n",
+ )
+
+ with open(inp_file, mode="w") as f:
+ f.writelines(lines)
+
+ # Run the simulation with the GMSO MCF file
+ code, out, err = run_cassandra(cassandra, inp_file)
+
+ assert code == 0
+ assert "complete" in out
+
+ # Parse log files
+ with open("nvt_gmso.out.log", mode="r") as f:
+ lines = f.readlines()
+ energy = None
+ for line in lines:
+ if "Total system energy" in line:
+ energy = float(line.split()[-1])
+ break
+
+ with open("nvt_mbuild.out.log", mode="r") as f:
+ lines = f.readlines()
+ energy_ref = 0.0
+ for line in lines:
+ if "Total system energy" in line:
+ energy_ref = float(line.split()[-1])
+ break
+
+ assert np.isclose(energy, energy_ref, rtol=1e-3)
+
+ @pytest.mark.skipif(not has_parmed, reason="ParmEd is not installed")
+ def test_parmed_vs_gmso(self, parmed_ethane):
+ """
+ This test compares the output of a MCF file generated
+ by gmso to the output of a MCF file generated by parmed.
+ The Parmed MCF file is generated through mbuild.
+ """
+ from gmso.external.convert_parmed import from_parmed
+
+ mb.formats.cassandramcf.write_mcf(
+ parmed_ethane,
+ "parmed-ethane.mcf",
+ dihedral_style="opls",
+ angle_style="harmonic",
+ )
+
+ top = from_parmed(parmed_ethane)
+ write_mcf(top, "gmso-ethane.mcf")
+
+ mcf_data_pmd, mcf_idx_pmd = parse_mcf("parmed-ethane.mcf")
+ mcf_data_gmso, mcf_idx_gmso = parse_mcf("gmso-ethane.mcf")
+ skip_lines = [3]
+ float_pattern = r"[+-]?[0-9]*[.][0-9]*"
+ for i, (line_pmd, line_gmso) in enumerate(
+ zip(mcf_data_pmd, mcf_data_gmso)
+ ):
+ if i in skip_lines:
+ continue
+
+ pmd_parms = np.array(
+ re.findall(float_pattern, " ".join(line_pmd)), dtype=np.float64
+ )
+ gmso_parms = np.array(
+ re.findall(float_pattern, " ".join(line_gmso)), dtype=np.float64
+ )
+
+ assert np.allclose(pmd_parms, gmso_parms, rtol=1e-3)
+
+ # TODO: can we read top.mcf into the simulation somehow?
+ # TODO: or use it to validate that the simulation was done
+ # correctly?
+
+ def test_untyped_top(self):
+ pass
+
+ def test_top_with_ring(self, typed_benzene_ua_system):
+ top = typed_benzene_ua_system
+ write_mcf(top, "benzene-ua.mcf")
+
+ mcf_data, mcf_idx = parse_mcf("benzene-ua.mcf")
+
+ assert is_charge_neutral(mcf_data, mcf_idx)
+
+ assert mcf_data[mcf_idx["Atom_Info"] + 1][0] == "6"
+
+ for idx in range(0, 6):
+ last_label = mcf_data[mcf_idx["Atom_Info"] + 2 + idx][-1]
+ assert last_label == "ring"
+
+ assert mcf_data[mcf_idx["Bond_Info"] + 1][0] == "6"
+ assert mcf_data[mcf_idx["Angle_Info"] + 1][0] == "6"
+ assert mcf_data[mcf_idx["Dihedral_Info"] + 1][0] == "6"
+
+ assert mcf_data[mcf_idx["Fragment_Info"] + 1][0] == "1"
+ frag_atoms = mcf_data[mcf_idx["Fragment_Info"] + 2][1:]
+ assert set(frag_atoms) == set([str(i) for i in range(1, 7)])
diff --git a/gmso/tests/test_potential_templates.py b/gmso/tests/test_potential_templates.py
index e8ac358ab..4487188af 100644
--- a/gmso/tests/test_potential_templates.py
+++ b/gmso/tests/test_potential_templates.py
@@ -180,6 +180,12 @@ def test_periodic_improper_potential(self, templates):
"phi_eq": ud.angle,
}
+ def test_fixed_bond_potential(self, templates):
+ potential = templates["FixedBondPotential"]
+ assert potential.name == "FixedBondPotential"
+ assert potential.expression == sympy.sympify("DiracDelta(r-r_eq)")
+ assert potential.independent_variables == {sympy.sympify("r")}
+
def test_harmonic_bond_potential(self, templates):
harmonic_bond_potential = templates["HarmonicBondPotential"]
assert harmonic_bond_potential.name == "HarmonicBondPotential"
@@ -195,6 +201,14 @@ def test_harmonic_bond_potential(self, templates):
"r_eq": ud.length,
}
+ def test_fixed_angle_potential(self, templates):
+ potential = templates["FixedAnglePotential"]
+ assert potential.name == "FixedAnglePotential"
+ assert potential.expression == sympy.sympify(
+ "DiracDelta(theta-theta_eq)"
+ )
+ assert potential.independent_variables == {sympy.sympify("theta")}
+
def test_harmonic_angle_potential(self, templates):
harmonic_angle_potential = templates["HarmonicAnglePotential"]
assert harmonic_angle_potential.name == "HarmonicAnglePotential"
diff --git a/gmso/tests/test_reference_xmls.py b/gmso/tests/test_reference_xmls.py
index 43b769b64..c8272d3ed 100644
--- a/gmso/tests/test_reference_xmls.py
+++ b/gmso/tests/test_reference_xmls.py
@@ -452,7 +452,7 @@ def test_ethylene_forcefield(self):
"4*epsilon*((sigma/r)**12 - (sigma/r)**6)",
"0.5 * k * (r-r_eq)**2",
"0.5 * k * (theta-theta_eq)**2",
- "c_0 + c_1 * cos(psi) + c_2 * cos(psi)**2 + c_3 * cos(psi)**3 + c_4 * cos(psi)**4 + c_5 * cos(psi)**5",
+ "c0 + c1 * cos(phi) + c2 * cos(phi)**2 + c3 * cos(phi)**3 + c4 * cos(phi)**4 + c5 * cos(phi)**5",
]
]
@@ -559,7 +559,7 @@ def test_ethylene_forcefield(self):
assert_allclose_units(
ethylene.dihedral_types[
"opls_144~opls_143~opls_143~opls_144"
- ].parameters["c_0"],
+ ].parameters["c0"],
58.576 * u.Unit("kJ/mol"),
rtol=1e-5,
atol=1e-8,
@@ -567,7 +567,7 @@ def test_ethylene_forcefield(self):
assert_allclose_units(
ethylene.dihedral_types[
"opls_144~opls_143~opls_143~opls_144"
- ].parameters["c_2"],
+ ].parameters["c2"],
-58.576 * u.Unit("kJ/mol"),
rtol=1e-5,
atol=1e-8,
@@ -575,7 +575,7 @@ def test_ethylene_forcefield(self):
assert_allclose_units(
ethylene.dihedral_types[
"opls_144~opls_143~opls_143~opls_144"
- ].parameters["c_5"],
+ ].parameters["c5"],
0.0 * u.Unit("kJ/mol"),
rtol=1e-5,
atol=1e-8,
diff --git a/gmso/utils/io.py b/gmso/utils/io.py
index 8dc3b632d..9405108e8 100644
--- a/gmso/utils/io.py
+++ b/gmso/utils/io.py
@@ -205,6 +205,14 @@ def import_(module):
except ImportError:
has_pandas = False
+try:
+ import mosdef_cassandra
+
+ has_cassandra = True
+ del mosdef_cassandra
+except ImportError:
+ has_cassandra = False
+
def run_from_ipython():
"""Verify that the code is running in an ipython kernel."""