diff --git a/gmso/tests/test_mcf.py b/gmso/tests/test_mcf.py index a0a3d592e..c6bbc7b83 100644 --- a/gmso/tests/test_mcf.py +++ b/gmso/tests/test_mcf.py @@ -1,4 +1,5 @@ import re +import subprocess import mbuild as mb import numpy as np @@ -48,10 +49,32 @@ def parse_mcf(filename): 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 @@ -369,9 +392,21 @@ def test_top_with_mixture(self): @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 - write_mcf(typed_ethane, "top.mcf") + # 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]]) @@ -381,11 +416,86 @@ def test_in_cassandra(self, typed_ethane): system=system, moveset=moveset, run_type="equilibration", - run_length=100, + 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="top", + 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() + 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() + 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): """