From 3327de217947493af7033168a8f6dda3330c985d Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 16 Sep 2024 12:40:38 -0500 Subject: [PATCH 1/9] initially create charmm parameter writer files and tests --- gmso/formats/prm_writer.py | 181 ++++++++++++++++++ gmso/lib/jsons/UreyBradleyAnglePotential.json | 11 ++ gmso/tests/test_par.py | 49 +++++ 3 files changed, 241 insertions(+) create mode 100644 gmso/formats/prm_writer.py create mode 100644 gmso/lib/jsons/UreyBradleyAnglePotential.json create mode 100644 gmso/tests/test_par.py diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py new file mode 100644 index 000000000..84667ee92 --- /dev/null +++ b/gmso/formats/prm_writer.py @@ -0,0 +1,181 @@ +"""CHARMM Par format.""" + +import warnings +from gmso.formats.formats_registry import saves_as + +__all__ = ["write_prm"] + +@saves_as("prm", "par") +def write_prm(topology, filename): + """Write CHARMM Parameter .prm file given a parametrized structure. + + Notes + ----- + Follows format according to + https://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/ + node25.html + Furthermore, ParmEd can support writing CHARMM par, rtf, str files + by converting the parmed.Structure into parmed.CharmmParameterSet + + Parmed stores rmin/2 in "rmin" + """ + # ATOMS + with open(filename, "w") as f: + f.write("ATOMS\n") + unique_atoms = set() + for site in topology.sites: + unique_atoms.add((site.atom_type.name, site.mass)) + for atom in unique_atoms: + f.write("MASS -1 {:8s} {:8.4f}\n".format(atom[0], atom[1])) + + f.write("\nBONDS\n") + unique_bonds = set() + for bond in topology.bonds: + atom1, atom2 = bond.connection_members + unique_bonds.add( + ( + atom1.atom_type.name, + atom2.atom_type.name, + bond.bond_type, + ) + ) + + for bond in unique_bonds: + # TODO: only works for harmonic bonds + f.write( + "{:8s} {:8s} {:.5f} {:.5f}\n".format( + bond[0], bond[1], bond[2].k, bond[2].req + ) + ) + + f.write("\nANGLES\n") + unique_harmonic_angles = set() + unique_ubs = set() + for angle in topology.angles: + atom1. atom2, atom3 = angle.connection_members + unique_ubs.add( + atom1.atom_type.name, + atom2.atom_type.name, + atom3.atom_type.name, + angle.angle_type, + angle.angle_type.name == "UreyBradleyAnglePotential" + ) + + for angle in unique_angles: + if angle[4]: # urey_bradley flag + f.write( + "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format( + angle[0], + angle[1], + angle[2], + angle[3].k, + angle[3].theteq, + angle[3].kub, + angle[3].r_eq, + ) + else: # assume harmonic angle + f.write( + "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format( + angle[0], + angle[1], + angle[2], + angle[3].k, + angle[3].theteq, + ) + + + # These dihedrals need to be PeriodicTorsion Style (Charmm style) + if len(topology.rb_torsions) > 0: + warnings.warn("RB Torsions detected, but unsupported in par writer") + f.write("\nDIHEDRALS\n") + unique_dihedrals = set() + for dihedral in topology.dihedrals: + # works for PeriodicTorsion + atom1, atom2, atom3, atom4 = dihedral.connection_members + unique_dihedrals.add( + ( + atom1.atom_type.name, + atom2.atom_type.name, + atom3.atom_type.name, + atom4.atom_type.name, + dihedral.dihedral_type, + ) + ) + + for dihedral in unique_dihedrals: + f.write( + "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( + dihedral[0], + dihedral[1], + dihedral[2], + dihedral[3], + dihedral[4].k, + dihedral[4].n, + dihedral[4].phi_eq, + ) + ) + + # TODO: No support for harmonic impropers + + f.write("\nIMPROPER\n") + unique_impropers = set() + for improper in topology.impropers: + atom1, atom2, atom3, atom4 = improper.connection_members + unique_impropers.add( + ( + atom1.atom_type.name, + atom2.atom_type.name, + atom3.atom_type.name, + atom4.atom_type.name, + improper.type, + ) + ) + for improper in unique_impropers: + # order of impropers goes + f.write( + "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( + improper[2], + improper[0], + improper[1], + improper[3], + improper[4].psi_k, + 0, + improper[4].psi_eq, + ) + ) + + sc_nb = [a for a in scnb] + if len(sc_nb) > 1: + warnings.warn( + "Multiple 1-4 LJ scalings were detected, " + "defaulting to first LJ scaling detected, {}".format(sc_nb[0]) + ) + sc_nb = sc_nb[0] + elif len(sc_nb) == 1: + sc_nb = sc_nb[0] + elif len(sc_nb) == 0: + warnings.warn("No 1-4 LJ scaling was detected, defaulting 1") + sc_nb = 1.0 + + f.write("\nNONBONDED\n") + unique_atypes = set() + for atom in topology.atoms: + unique_atypes.add(atom.atom_type) + for atype in unique_atypes: + # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) + f.write( + "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format( + atype.name, + 0.0, + -1 * atype.epsilon, + atype.rmin, + 0.0, + -1 * sc_nb * atype.epsilon, + atype.rmin, + ) + ) + + if topology.has_NBFIX(): + warnings.warn("NBFixes detected but unsupported in par writer") + + f.write("\nEND") diff --git a/gmso/lib/jsons/UreyBradleyAnglePotential.json b/gmso/lib/jsons/UreyBradleyAnglePotential.json new file mode 100644 index 000000000..69bd87a6d --- /dev/null +++ b/gmso/lib/jsons/UreyBradleyAnglePotential.json @@ -0,0 +1,11 @@ +{ + "name": "UreyBradleyPotential", + "expression": "0.5 * k * (theta-theta_eq)**2 + 0.5 * kub * (r-r_eq)**2", + "independent_variables": ["theta", "r"], + "expected_parameters_dimensions": { + "k":"energy/angle**2", + "theta_eq": "angle", + "kub": "energy/length**2", + "r_eq": "length" + } +} diff --git a/gmso/tests/test_par.py b/gmso/tests/test_par.py new file mode 100644 index 000000000..3b6675aba --- /dev/null +++ b/gmso/tests/test_par.py @@ -0,0 +1,49 @@ +import numpy as np +import pytest + +from gmso.formats.prm_writer import write_par +import mbuild as mb +from mbuild.tests.base_test import BaseTest +from mbuild.utils.io import get_fn, has_foyer + + +@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") +class TestPar(BaseTest): + @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") + def test_save_charmm(self): + cmpd = mb.load(get_fn("charmm_dihedral.mol2")) + for i in cmpd.particles(): + i.name = "_{}".format(i.name) + structure = cmpd.to_parmed( + box=cmpd.get_boundingbox(), + residues=set([p.parent.name for p in cmpd.particles()]), + ) + + from foyer import Forcefield + + ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")]) + structure = ff.apply(structure, assert_dihedral_params=False) + + write_par(structure, "charmm_dihedral.par") + + @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") + def test_save_forcefield(self, ethane): + ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") + + def test_par_parameters(self, ethane): + ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") + from parmed.charmm import CharmmParameterSet + + pset = CharmmParameterSet.load_set(pfile="ethane-opls.par") + assert len(pset.bond_types) == 3 + assert len(pset.angle_types) == 3 + assert len(pset.atom_types) == 2 + + def test_parameter_forms(self, ethane): + # test that the values are correct + # write mbuild file and compare to gmso file + pass + def test_raise_errors(self): + # only takes harmonic bonds + # only takes harmonic of ub angles + pass From 25f893f6a83fb60adc7c8714c8ef706107073070 Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 16 Sep 2024 17:42:10 -0500 Subject: [PATCH 2/9] function charmm file writers --- gmso/formats/__init__.py | 1 + gmso/formats/prm_writer.py | 179 +++++++++----------------- gmso/tests/test_par.py | 46 +++---- gmso/utils/files/charmm_dihedral.mol2 | 17 +++ gmso/utils/files/charmm_truncated.xml | 44 +++++++ 5 files changed, 144 insertions(+), 143 deletions(-) create mode 100644 gmso/utils/files/charmm_dihedral.mol2 create mode 100644 gmso/utils/files/charmm_truncated.xml diff --git a/gmso/formats/__init__.py b/gmso/formats/__init__.py index b961c299c..804f3e05b 100644 --- a/gmso/formats/__init__.py +++ b/gmso/formats/__init__.py @@ -10,6 +10,7 @@ from .lammpsdata import write_lammpsdata from .mcf import write_mcf from .mol2 import read_mol2 +from .prm_writer import write_prm from .top import write_top from .xyz import read_xyz, write_xyz diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index 84667ee92..88e97457f 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -1,11 +1,10 @@ """CHARMM Par format.""" -import warnings from gmso.formats.formats_registry import saves_as +from gmso.utils.units import LAMMPS_UnitSystems -__all__ = ["write_prm"] -@saves_as("prm", "par") +@saves_as(".prm", ".par") def write_prm(topology, filename): """Write CHARMM Parameter .prm file given a parametrized structure. @@ -19,163 +18,107 @@ def write_prm(topology, filename): Parmed stores rmin/2 in "rmin" """ + unit_system = LAMMPS_UnitSystems("real") # ang, kcal/mol, amu # ATOMS with open(filename, "w") as f: f.write("ATOMS\n") unique_atoms = set() for site in topology.sites: - unique_atoms.add((site.atom_type.name, site.mass)) - for atom in unique_atoms: - f.write("MASS -1 {:8s} {:8.4f}\n".format(atom[0], atom[1])) - - f.write("\nBONDS\n") - unique_bonds = set() - for bond in topology.bonds: - atom1, atom2 = bond.connection_members - unique_bonds.add( + unique_atoms.add( ( - atom1.atom_type.name, - atom2.atom_type.name, - bond.bond_type, + site.atom_type.name, + unit_system.convert_parameter( + site.atom_type.mass, n_decimals=6, name="mass" + ), ) ) + for atom in unique_atoms: + f.write("MASS -1 {:8s} {:8s}\n".format(atom[0], atom[1])) - for bond in unique_bonds: - # TODO: only works for harmonic bonds + # TODO: Works for harmonic bonds + f.write("\nBONDS\n") + for bond in topology.bond_types: + atom1, atom2 = bond.member_types f.write( "{:8s} {:8s} {:.5f} {:.5f}\n".format( - bond[0], bond[1], bond[2].k, bond[2].req + atom1, atom2, bond.parameters["k"], bond.parameters["r_eq"] ) ) f.write("\nANGLES\n") - unique_harmonic_angles = set() - unique_ubs = set() - for angle in topology.angles: - atom1. atom2, atom3 = angle.connection_members - unique_ubs.add( - atom1.atom_type.name, - atom2.atom_type.name, - atom3.atom_type.name, - angle.angle_type, - angle.angle_type.name == "UreyBradleyAnglePotential" - ) - - for angle in unique_angles: - if angle[4]: # urey_bradley flag + for angle in topology.angle_types: + atom1, atom2, atom3 = angle.member_types + if angle.name == "UreyBradleyAnglePotential": f.write( "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format( - angle[0], - angle[1], - angle[2], - angle[3].k, - angle[3].theteq, - angle[3].kub, - angle[3].r_eq, + atom1, + atom2, + atom3, + angle.parameters["k"], + angle.parameters["theta_eq"], + angle.parameters["kub"], + angle.parameters["r_eq"], ) - else: # assume harmonic angle + ) + else: # assume harmonic style: f.write( "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format( - angle[0], - angle[1], - angle[2], - angle[3].k, - angle[3].theteq, + atom1, + atom2, + atom3, + angle.parameters["k"], + angle.parameters["theta_eq"], ) - + ) # These dihedrals need to be PeriodicTorsion Style (Charmm style) - if len(topology.rb_torsions) > 0: - warnings.warn("RB Torsions detected, but unsupported in par writer") f.write("\nDIHEDRALS\n") - unique_dihedrals = set() - for dihedral in topology.dihedrals: + for dihedral in topology.dihedral_types: # works for PeriodicTorsion - atom1, atom2, atom3, atom4 = dihedral.connection_members - unique_dihedrals.add( - ( - atom1.atom_type.name, - atom2.atom_type.name, - atom3.atom_type.name, - atom4.atom_type.name, - dihedral.dihedral_type, - ) - ) - - for dihedral in unique_dihedrals: + atom1, atom2, atom3, atom4 = dihedral.member_types f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( - dihedral[0], - dihedral[1], - dihedral[2], - dihedral[3], - dihedral[4].k, - dihedral[4].n, - dihedral[4].phi_eq, + "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format( + atom1, + atom2, + atom3, + atom4, + dihedral.parameters["k"][0], + dihedral.parameters["n"][0], + dihedral.parameters["phi_eq"][0], ) ) # TODO: No support for harmonic impropers f.write("\nIMPROPER\n") - unique_impropers = set() - for improper in topology.impropers: - atom1, atom2, atom3, atom4 = improper.connection_members - unique_impropers.add( - ( - atom1.atom_type.name, - atom2.atom_type.name, - atom3.atom_type.name, - atom4.atom_type.name, - improper.type, - ) - ) - for improper in unique_impropers: - # order of impropers goes + for improper in topology.improper_types: + atom1, atom2, atom3, atom4 = improper.member_types f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5d} {:.5f}\n".format( - improper[2], - improper[0], - improper[1], - improper[3], - improper[4].psi_k, - 0, - improper[4].psi_eq, + "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format( + atom1, + atom2, + atom3, + atom4, + improper.parameters["phi_k"][0], + improper.parameters["n"][0], + improper.parameters["phi_eq"][0], ) ) - sc_nb = [a for a in scnb] - if len(sc_nb) > 1: - warnings.warn( - "Multiple 1-4 LJ scalings were detected, " - "defaulting to first LJ scaling detected, {}".format(sc_nb[0]) - ) - sc_nb = sc_nb[0] - elif len(sc_nb) == 1: - sc_nb = sc_nb[0] - elif len(sc_nb) == 0: - warnings.warn("No 1-4 LJ scaling was detected, defaulting 1") - sc_nb = 1.0 - f.write("\nNONBONDED\n") - unique_atypes = set() - for atom in topology.atoms: - unique_atypes.add(atom.atom_type) - for atype in unique_atypes: + nonbonded14 = topology.scaling_factors[0][2] + + for atype in topology.atom_types: # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) f.write( "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format( atype.name, + 0.0, # ignored, + atype.parameters["epsilon"], + atype.parameters["sigma"], 0.0, - -1 * atype.epsilon, - atype.rmin, - 0.0, - -1 * sc_nb * atype.epsilon, - atype.rmin, + atype.parameters["epsilon"] * nonbonded14, + atype.parameters["sigma"] * nonbonded14, ) ) - - if topology.has_NBFIX(): - warnings.warn("NBFixes detected but unsupported in par writer") - f.write("\nEND") diff --git a/gmso/tests/test_par.py b/gmso/tests/test_par.py index 3b6675aba..b999e3728 100644 --- a/gmso/tests/test_par.py +++ b/gmso/tests/test_par.py @@ -1,37 +1,32 @@ -import numpy as np -import pytest +from gmso import ForceField, Topology +from gmso.parameterization import apply +from gmso.tests.base_test import BaseTest +from gmso.utils.io import get_fn -from gmso.formats.prm_writer import write_par -import mbuild as mb -from mbuild.tests.base_test import BaseTest -from mbuild.utils.io import get_fn, has_foyer +# TODO: Sorting of atomtypes info in connection types +# TODO: Test of correct potential forms +# TODO: Handle iterating over unique types +# TODO: Handle multiple values in charmm dihedral k as np array -@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") class TestPar(BaseTest): - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") def test_save_charmm(self): - cmpd = mb.load(get_fn("charmm_dihedral.mol2")) - for i in cmpd.particles(): - i.name = "_{}".format(i.name) - structure = cmpd.to_parmed( - box=cmpd.get_boundingbox(), - residues=set([p.parent.name for p in cmpd.particles()]), - ) - - from foyer import Forcefield + top = Topology.load(get_fn("charmm_dihedral.mol2")) + for site in top.sites: + site.name = f"_{site.name}" - ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")]) - structure = ff.apply(structure, assert_dihedral_params=False) + ff = ForceField(get_fn("charmm_truncated.xml")) + ptop = apply( + top, ff, identify_connections=True, ignore_params=["dihedral", "improper"] + ) + ptop.save("charmm_dihedral.prm") - write_par(structure, "charmm_dihedral.par") + def test_par_parameters(self, ethane, oplsaa_forcefield): + from gmso.parameterization import apply - @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed") - def test_save_forcefield(self, ethane): - ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") + ptop = apply(ethane, oplsaa_forcefield) - def test_par_parameters(self, ethane): - ethane.save(filename="ethane-opls.par", forcefield_name="oplsaa") + ptop.save(filename="ethane-opls.par") from parmed.charmm import CharmmParameterSet pset = CharmmParameterSet.load_set(pfile="ethane-opls.par") @@ -43,6 +38,7 @@ def test_parameter_forms(self, ethane): # test that the values are correct # write mbuild file and compare to gmso file pass + def test_raise_errors(self): # only takes harmonic bonds # only takes harmonic of ub angles diff --git a/gmso/utils/files/charmm_dihedral.mol2 b/gmso/utils/files/charmm_dihedral.mol2 new file mode 100644 index 000000000..141e13a9b --- /dev/null +++ b/gmso/utils/files/charmm_dihedral.mol2 @@ -0,0 +1,17 @@ +@MOLECULE +DPPC +4 3 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 13.9528 17.1622 31.1323 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 CTL1 -0.1558 6.4170 4.4541 CTL1 1 DPPC 1.500000 + 2 OSL -0.7278 7.6038 3.8783 OSL 1 DPPC -0.780000 + 3 PL -0.5108 6.0856 5.8561 PL 1 DPPC -0.780000 + 4 OSLP -0.2550 5.1063 3.5611 OSLP 1 DPPC -0.570000 + +@BOND + 1 1 2 1 + 2 2 3 1 + 3 1 4 1 diff --git a/gmso/utils/files/charmm_truncated.xml b/gmso/utils/files/charmm_truncated.xml new file mode 100644 index 000000000..8df778965 --- /dev/null +++ b/gmso/utils/files/charmm_truncated.xml @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 13e0d6e02bf425a16913c942fdf10d4b67f38ad0 Mon Sep 17 00:00:00 2001 From: CalCraven Date: Fri, 27 Sep 2024 07:56:37 -0500 Subject: [PATCH 3/9] Added support and testing files for urey bradley potentials --- gmso/formats/prm_writer.py | 140 ++++++++++--- .../LAMMPSHarmonicImproperPotential.json | 9 + gmso/lib/jsons/UreyBradleyAnglePotential.json | 4 +- gmso/tests/{test_par.py => test_prm.py} | 31 ++- gmso/utils/files/charmm_amoeba.xml | 192 ++++++++++++++++++ gmso/utils/files/charmm_improper.mol2 | 21 ++ 6 files changed, 358 insertions(+), 39 deletions(-) create mode 100644 gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json rename gmso/tests/{test_par.py => test_prm.py} (55%) create mode 100644 gmso/utils/files/charmm_amoeba.xml create mode 100644 gmso/utils/files/charmm_improper.mol2 diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index 88e97457f..d47696048 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -2,6 +2,39 @@ from gmso.formats.formats_registry import saves_as from gmso.utils.units import LAMMPS_UnitSystems +from gmso.utils.compatibility import check_compatibility +from gmso.core.views import PotentialFilters +from gmso.lib.potential_templates import PotentialTemplateLibrary + + +def _validate_potential_compatibility(top): + """Check compatability of topology object potentials with LAMMPSDATA format.""" + pfilter = PotentialFilters.UNIQUE_EXPRESSION + pot_types = check_compatibility( + top, _accepted_potentials(), site_pfilter=pfilter, conn_pfilter=pfilter + ) + return pot_types + +def _accepted_potentials(): + """List of accepted potentials that LAMMPS can support.""" + templates = PotentialTemplateLibrary() + lennard_jones_potential = templates["LennardJonesPotential"] + harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] + harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] + ub_angle_potential = templates["UreyBradleyAnglePotential"] + periodic_torsion_potential = templates["PeriodicTorsionPotential"] + harmonic_improper_potential = templates["LAMMPSHarmonicImproperPotential"] + accepted_potentialsList = [ + lennard_jones_potential, + harmonic_bond_potential, + harmonic_angle_potential, + ub_angle_potential, + periodic_torsion_potential, + harmonic_improper_potential, + ] + return accepted_potentialsList + + @saves_as(".prm", ".par") @@ -18,6 +51,11 @@ def write_prm(topology, filename): Parmed stores rmin/2 in "rmin" """ + # Validation + potentialsMap = _validate_potential_compatibility(topology) + + + unit_system = LAMMPS_UnitSystems("real") # ang, kcal/mol, amu # ATOMS with open(filename, "w") as f: @@ -28,7 +66,7 @@ def write_prm(topology, filename): ( site.atom_type.name, unit_system.convert_parameter( - site.atom_type.mass, n_decimals=6, name="mass" + site.atom_type.mass, n_decimals=3, name="mass" ), ) ) @@ -40,8 +78,15 @@ def write_prm(topology, filename): for bond in topology.bond_types: atom1, atom2 = bond.member_types f.write( - "{:8s} {:8s} {:.5f} {:.5f}\n".format( - atom1, atom2, bond.parameters["k"], bond.parameters["r_eq"] + "{:8s} {:8s} {:7s} {:7s}\n".format( + + atom1, atom2, + unit_system.convert_parameter( + bond.parameters["k"], n_decimals=3 + ), + unit_system.convert_parameter( + bond.parameters["r_eq"], n_decimals=3 + ) ) ) @@ -50,24 +95,37 @@ def write_prm(topology, filename): atom1, atom2, atom3 = angle.member_types if angle.name == "UreyBradleyAnglePotential": f.write( - "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format( + "{:8s} {:8s} {:8s} {:7s} {:7s} {:7s} {:7s}\n".format( atom1, atom2, atom3, - angle.parameters["k"], - angle.parameters["theta_eq"], - angle.parameters["kub"], - angle.parameters["r_eq"], + unit_system.convert_parameter( + angle.parameters["k"], n_decimals=3 + ), + unit_system.convert_parameter( + angle.parameters["theta_eq"], n_decimals=3, name="theta_eq" # necessary for conversion to degrees not radians + ), + unit_system.convert_parameter( + angle.parameters["kub"], n_decimals=3 + ), + unit_system.convert_parameter( + angle.parameters["r_eq"], n_decimals=3 + ), ) ) else: # assume harmonic style: f.write( - "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format( + "{:8s} {:8s} {:8s} {:7s} {:7s}\n".format( atom1, atom2, atom3, - angle.parameters["k"], - angle.parameters["theta_eq"], + unit_system.convert_parameter( + angle.parameters["k"], n_decimals=3 + ), + unit_system.convert_parameter( + angle.parameters["theta_eq"], n_decimals=3, + name="theta_eq" + ), ) ) @@ -76,32 +134,44 @@ def write_prm(topology, filename): for dihedral in topology.dihedral_types: # works for PeriodicTorsion atom1, atom2, atom3, atom4 = dihedral.member_types - f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format( - atom1, - atom2, - atom3, - atom4, - dihedral.parameters["k"][0], - dihedral.parameters["n"][0], - dihedral.parameters["phi_eq"][0], + variable_dtypes = ["k", "n", "phi_eq"] + zipped_params = (dihedral.parameters[x] for x in variable_dtypes) + for k, n, phi_eq in zip(*zipped_params): + f.write( + "{:8s} {:8s} {:8s} {:8s} {:7s} {:7s} {:7s}\n".format( + atom1, + atom2, + atom3, + atom4, + unit_system.convert_parameter( + k, n_decimals=3 + ), + unit_system.convert_parameter( + n, n_decimals=3 + ), + unit_system.convert_parameter( + phi_eq, n_decimals=3, name="phi_eq" + ), + ) ) - ) - # TODO: No support for harmonic impropers f.write("\nIMPROPER\n") for improper in topology.improper_types: atom1, atom2, atom3, atom4 = improper.member_types f.write( - "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format( + "{:8s} {:8s} {:8s} {:8s} {:.5s} {:5.3f} {:.5s}\n".format( atom1, atom2, atom3, atom4, - improper.parameters["phi_k"][0], - improper.parameters["n"][0], - improper.parameters["phi_eq"][0], + unit_system.convert_parameter( + improper.parameters["k"], n_decimals=3 + ), + 0.0, + unit_system.convert_parameter( + improper.parameters["phi_eq"], n_decimals=3 + ), ) ) @@ -111,14 +181,22 @@ def write_prm(topology, filename): for atype in topology.atom_types: # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) f.write( - "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format( + "{:8s} {:8.3f} {:8s} {:8s} {:8.3f} {:8s} {:8s}\n".format( atype.name, 0.0, # ignored, - atype.parameters["epsilon"], - atype.parameters["sigma"], + unit_system.convert_parameter( + atype.parameters["epsilon"], n_decimals=3 + ), + unit_system.convert_parameter( + atype.parameters["sigma"], n_decimals=3 + ), 0.0, - atype.parameters["epsilon"] * nonbonded14, - atype.parameters["sigma"] * nonbonded14, + unit_system.convert_parameter( + atype.parameters["epsilon"], n_decimals=3 + ), + unit_system.convert_parameter( + atype.parameters["sigma"], n_decimals=3 + ), ) ) f.write("\nEND") diff --git a/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json b/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json new file mode 100644 index 000000000..f3a33261a --- /dev/null +++ b/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json @@ -0,0 +1,9 @@ +{ + "name": "HarmonicImproperPotential", + "expression": "k * (phi - phi_eq)**2", + "independent_variables": "phi", + "expected_parameters_dimensions": { + "k": "energy/angle**2", + "phi_eq": "angle" + } +} diff --git a/gmso/lib/jsons/UreyBradleyAnglePotential.json b/gmso/lib/jsons/UreyBradleyAnglePotential.json index 69bd87a6d..55c896913 100644 --- a/gmso/lib/jsons/UreyBradleyAnglePotential.json +++ b/gmso/lib/jsons/UreyBradleyAnglePotential.json @@ -1,7 +1,7 @@ { - "name": "UreyBradleyPotential", + "name": "UreyBradleyAnglePotential", "expression": "0.5 * k * (theta-theta_eq)**2 + 0.5 * kub * (r-r_eq)**2", - "independent_variables": ["theta", "r"], + "independent_variables": "theta,r", "expected_parameters_dimensions": { "k":"energy/angle**2", "theta_eq": "angle", diff --git a/gmso/tests/test_par.py b/gmso/tests/test_prm.py similarity index 55% rename from gmso/tests/test_par.py rename to gmso/tests/test_prm.py index b999e3728..e059e98e7 100644 --- a/gmso/tests/test_par.py +++ b/gmso/tests/test_prm.py @@ -1,7 +1,10 @@ +import pytest + from gmso import ForceField, Topology from gmso.parameterization import apply from gmso.tests.base_test import BaseTest from gmso.utils.io import get_fn +from gmso.exceptions import EngineIncompatibilityError # TODO: Sorting of atomtypes info in connection types # TODO: Test of correct potential forms @@ -17,22 +20,38 @@ def test_save_charmm(self): ff = ForceField(get_fn("charmm_truncated.xml")) ptop = apply( - top, ff, identify_connections=True, ignore_params=["dihedral", "improper"] + top, ff, identify_connections=True, ignore_params=["improper"] ) ptop.save("charmm_dihedral.prm") - def test_par_parameters(self, ethane, oplsaa_forcefield): + def test_par_parameters(self): from gmso.parameterization import apply + top = Topology.load(get_fn("charmm_improper.mol2")) + ff = ForceField(get_fn("charmm_amoeba.xml")) + for site in top.sites: + site.name = f"_{site.name}" - ptop = apply(ethane, oplsaa_forcefield) + ptop = apply(top, ff, identify_connections=True, + ignore_params=["improper", "dihedral"] + ) ptop.save(filename="ethane-opls.par") from parmed.charmm import CharmmParameterSet pset = CharmmParameterSet.load_set(pfile="ethane-opls.par") - assert len(pset.bond_types) == 3 - assert len(pset.angle_types) == 3 - assert len(pset.atom_types) == 2 + assert len(pset.bond_types) == 10 # each bond counted twice + assert len(pset.angle_types) == 10 # each angle counted twice + assert len(pset.dihedral_types) == 2 + assert len(pset.improper_types) == 1 + assert len(pset.atom_types) == 4 + + def test_prm_incompatibile_types(self, ethane, oplsaa_forcefield): + from gmso.parameterization import apply + + ptop = apply(ethane, oplsaa_forcefield, identify_connections=True) + with pytest.raises(EngineIncompatibilityError): + ptop.save(filename="ethane-opls.par") + def test_parameter_forms(self, ethane): # test that the values are correct diff --git a/gmso/utils/files/charmm_amoeba.xml b/gmso/utils/files/charmm_amoeba.xml new file mode 100644 index 000000000..1dda2feb1 --- /dev/null +++ b/gmso/utils/files/charmm_amoeba.xml @@ -0,0 +1,192 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 69.0 + 138.0 + 138.0 + + + 0.0 + 45.0 + 90.0 + + + 1.0 + 2.0 + 3.0 + + + + + + + + + + + + + + + + + + + + + diff --git a/gmso/utils/files/charmm_improper.mol2 b/gmso/utils/files/charmm_improper.mol2 new file mode 100644 index 000000000..7394123eb --- /dev/null +++ b/gmso/utils/files/charmm_improper.mol2 @@ -0,0 +1,21 @@ +@MOLECULE +DPPC +4 3 1 0 1 +SMALL +USER_CHARGES +@CRYSIN + 13.9528 17.1622 31.1323 90.0000 90.0000 90.0000 1 1 +@ATOM + 1 CTL1 -0.1558 6.4170 4.4541 CTL1 1 DPPC 1.500000 + 2 OSL -0.7278 7.6038 3.8783 OSL 1 DPPC -0.780000 + 3 PL -0.5108 6.0856 5.8561 PL 1 DPPC -0.780000 + 4 OSLP -0.2550 5.1063 3.5611 OSLP 1 DPPC -0.570000 + 5 OSLP -0.2550 5.1063 3.5611 OSLP 1 DPPC -0.570000 + 6 CTL1 -0.1558 6.4170 4.4541 CTL1 1 DPPC 1.500000 + +@BOND + 1 1 2 1 + 2 2 3 1 + 3 1 4 1 + 4 3 5 1 + 5 3 6 1 From 2c7df9c23ccdedaf1fdeeb52292b09826a9aedb2 Mon Sep 17 00:00:00 2001 From: CalCraven Date: Fri, 27 Sep 2024 07:59:18 -0500 Subject: [PATCH 4/9] Add extra checking for reading independent variables from forcefield as a list --- gmso/lib/potential_templates.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/gmso/lib/potential_templates.py b/gmso/lib/potential_templates.py index 9a3564eeb..5a075a8a4 100644 --- a/gmso/lib/potential_templates.py +++ b/gmso/lib/potential_templates.py @@ -69,7 +69,11 @@ def __init__( potential_expression=None, expected_parameters_dimensions=None, ): - if not isinstance(independent_variables, set): + if isinstance(independent_variables, set): + pass + elif isinstance(independent_variables, list): + independent_variables = set(independent_variables) + elif isinstance(independent_variables, str): independent_variables = set(independent_variables.split(",")) if potential_expression is None: From 7bfbb3f63d770bca4d73fa2a952fc26fc7be048a Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 30 Sep 2024 10:00:58 -0400 Subject: [PATCH 5/9] Add compatibility checks for charmm parm files --- gmso/formats/prm_writer.py | 66 ++++++++++++------- .../LAMMPSHarmonicImproperPotential.json | 2 +- gmso/utils/files/charmm_amoeba.xml | 4 +- 3 files changed, 44 insertions(+), 28 deletions(-) diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index d47696048..109639c2c 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -1,10 +1,10 @@ """CHARMM Par format.""" -from gmso.formats.formats_registry import saves_as -from gmso.utils.units import LAMMPS_UnitSystems -from gmso.utils.compatibility import check_compatibility from gmso.core.views import PotentialFilters +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.units import LAMMPS_UnitSystems def _validate_potential_compatibility(top): @@ -15,6 +15,7 @@ def _validate_potential_compatibility(top): ) return pot_types + def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() @@ -28,17 +29,15 @@ def _accepted_potentials(): lennard_jones_potential, harmonic_bond_potential, harmonic_angle_potential, - ub_angle_potential, + ub_angle_potential, periodic_torsion_potential, harmonic_improper_potential, ] return accepted_potentialsList - - @saves_as(".prm", ".par") -def write_prm(topology, filename): +def write_prm(topology, filename, strict_potentials=False): """Write CHARMM Parameter .prm file given a parametrized structure. Notes @@ -52,10 +51,20 @@ def write_prm(topology, filename): Parmed stores rmin/2 in "rmin" """ # Validation + # TODO: Use strict_x, (e.g. x=bonds) to validate what topology attrs to convert + if not strict_potentials: + default_parameterMaps = { # TODO: sites are not checked currently because gmso + # doesn't store pair potential eqn the same way as the connections. + "impropers": ["LAMMPSHarmonicImproperPotential"], + "dihedrals": ["PeriodicTorsionPotential"], + "angles": ["LAMMPSHarmonicAnglePotential", "UreyBradleyAnglePotential"], + "bonds": ["LAMMPSHarmonicBondPotential"], + # "sites":"LennardJonesPotential", + # "sites":"CoulombicPotential" + } + _try_default_potential_conversions(topology, default_parameterMaps) potentialsMap = _validate_potential_compatibility(topology) - - unit_system = LAMMPS_UnitSystems("real") # ang, kcal/mol, amu # ATOMS with open(filename, "w") as f: @@ -79,14 +88,12 @@ def write_prm(topology, filename): atom1, atom2 = bond.member_types f.write( "{:8s} {:8s} {:7s} {:7s}\n".format( - - atom1, atom2, - unit_system.convert_parameter( - bond.parameters["k"], n_decimals=3 - ), + atom1, + atom2, + unit_system.convert_parameter(bond.parameters["k"], n_decimals=3), unit_system.convert_parameter( bond.parameters["r_eq"], n_decimals=3 - ) + ), ) ) @@ -103,7 +110,9 @@ def write_prm(topology, filename): angle.parameters["k"], n_decimals=3 ), unit_system.convert_parameter( - angle.parameters["theta_eq"], n_decimals=3, name="theta_eq" # necessary for conversion to degrees not radians + angle.parameters["theta_eq"], + n_decimals=3, + name="theta_eq", # necessary for conversion to degrees not radians ), unit_system.convert_parameter( angle.parameters["kub"], n_decimals=3 @@ -123,8 +132,7 @@ def write_prm(topology, filename): angle.parameters["k"], n_decimals=3 ), unit_system.convert_parameter( - angle.parameters["theta_eq"], n_decimals=3, - name="theta_eq" + angle.parameters["theta_eq"], n_decimals=3, name="theta_eq" ), ) ) @@ -143,19 +151,14 @@ def write_prm(topology, filename): atom2, atom3, atom4, - unit_system.convert_parameter( - k, n_decimals=3 - ), - unit_system.convert_parameter( - n, n_decimals=3 - ), + unit_system.convert_parameter(k, n_decimals=3), + unit_system.convert_parameter(n, n_decimals=3), unit_system.convert_parameter( phi_eq, n_decimals=3, name="phi_eq" ), ) ) - f.write("\nIMPROPER\n") for improper in topology.improper_types: atom1, atom2, atom3, atom4 = improper.member_types @@ -200,3 +203,16 @@ def write_prm(topology, filename): ) ) f.write("\nEND") + + +def _try_default_potential_conversions(top, potentialsDict): + """Take a topology a convert all potentials to the style in potentialDict.""" "" + for pot_container in potentialsDict: + containerStr = pot_container[:-1] + "_types" + if getattr(top, containerStr): + for potential in potentialsDict[pot_container]: + top.convert_potential_styles({pot_container: potential}) + elif getattr(top, pot_container): + raise AttributeError( + f"Missing parameters in {pot_container} for {top.get_untyped(pot_container)}" + ) diff --git a/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json b/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json index f3a33261a..ab0246b7a 100644 --- a/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json +++ b/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json @@ -1,5 +1,5 @@ { - "name": "HarmonicImproperPotential", + "name": "LAMMPSHarmonicImproperPotential", "expression": "k * (phi - phi_eq)**2", "independent_variables": "phi", "expected_parameters_dimensions": { diff --git a/gmso/utils/files/charmm_amoeba.xml b/gmso/utils/files/charmm_amoeba.xml index 1dda2feb1..19d45cfcb 100644 --- a/gmso/utils/files/charmm_amoeba.xml +++ b/gmso/utils/files/charmm_amoeba.xml @@ -138,8 +138,8 @@ - - + + From ea5ec2894b84fe98f7f207447b6bf5600e5bf399 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 30 Sep 2024 14:02:23 +0000 Subject: [PATCH 6/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gmso/tests/test_prm.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/gmso/tests/test_prm.py b/gmso/tests/test_prm.py index e059e98e7..a45aaca72 100644 --- a/gmso/tests/test_prm.py +++ b/gmso/tests/test_prm.py @@ -1,10 +1,10 @@ import pytest from gmso import ForceField, Topology +from gmso.exceptions import EngineIncompatibilityError from gmso.parameterization import apply from gmso.tests.base_test import BaseTest from gmso.utils.io import get_fn -from gmso.exceptions import EngineIncompatibilityError # TODO: Sorting of atomtypes info in connection types # TODO: Test of correct potential forms @@ -19,28 +19,27 @@ def test_save_charmm(self): site.name = f"_{site.name}" ff = ForceField(get_fn("charmm_truncated.xml")) - ptop = apply( - top, ff, identify_connections=True, ignore_params=["improper"] - ) + ptop = apply(top, ff, identify_connections=True, ignore_params=["improper"]) ptop.save("charmm_dihedral.prm") def test_par_parameters(self): from gmso.parameterization import apply + top = Topology.load(get_fn("charmm_improper.mol2")) ff = ForceField(get_fn("charmm_amoeba.xml")) for site in top.sites: site.name = f"_{site.name}" - ptop = apply(top, ff, identify_connections=True, - ignore_params=["improper", "dihedral"] - ) + ptop = apply( + top, ff, identify_connections=True, ignore_params=["improper", "dihedral"] + ) ptop.save(filename="ethane-opls.par") from parmed.charmm import CharmmParameterSet pset = CharmmParameterSet.load_set(pfile="ethane-opls.par") - assert len(pset.bond_types) == 10 # each bond counted twice - assert len(pset.angle_types) == 10 # each angle counted twice + assert len(pset.bond_types) == 10 # each bond counted twice + assert len(pset.angle_types) == 10 # each angle counted twice assert len(pset.dihedral_types) == 2 assert len(pset.improper_types) == 1 assert len(pset.atom_types) == 4 @@ -52,7 +51,6 @@ def test_prm_incompatibile_types(self, ethane, oplsaa_forcefield): with pytest.raises(EngineIncompatibilityError): ptop.save(filename="ethane-opls.par") - def test_parameter_forms(self, ethane): # test that the values are correct # write mbuild file and compare to gmso file From 2f10015f23a71d48c7acde1bff012b5c076f69ff Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 30 Sep 2024 10:53:24 -0400 Subject: [PATCH 7/9] Fix unused variables in charmm writer and update epsilon in lj equation from 4*epsilon --- gmso/formats/prm_writer.py | 9 ++++++--- gmso/utils/files/charmm_amoeba.xml | 2 +- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index 109639c2c..6f595db76 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -20,6 +20,7 @@ def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() lennard_jones_potential = templates["LennardJonesPotential"] + lennard_jones_potential.expression /= 4 # no 4*epsilon term harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] ub_angle_potential = templates["UreyBradleyAnglePotential"] @@ -59,11 +60,11 @@ def write_prm(topology, filename, strict_potentials=False): "dihedrals": ["PeriodicTorsionPotential"], "angles": ["LAMMPSHarmonicAnglePotential", "UreyBradleyAnglePotential"], "bonds": ["LAMMPSHarmonicBondPotential"], - # "sites":"LennardJonesPotential", + "sites":["LennardJonesPotential"], # "sites":"CoulombicPotential" } _try_default_potential_conversions(topology, default_parameterMaps) - potentialsMap = _validate_potential_compatibility(topology) + _validate_potential_compatibility(topology) unit_system = LAMMPS_UnitSystems("real") # ang, kcal/mol, amu # ATOMS @@ -179,7 +180,7 @@ def write_prm(topology, filename, strict_potentials=False): ) f.write("\nNONBONDED\n") - nonbonded14 = topology.scaling_factors[0][2] + # nonbonded14 = topology.scaling_factors[0][2] for atype in topology.atom_types: # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) @@ -202,6 +203,8 @@ def write_prm(topology, filename, strict_potentials=False): ), ) ) + # TODO: Add NONBONDED section: NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch - + # !adm jr., 5/08/91, suggested cutoff scheme f.write("\nEND") diff --git a/gmso/utils/files/charmm_amoeba.xml b/gmso/utils/files/charmm_amoeba.xml index 19d45cfcb..c3d464e67 100644 --- a/gmso/utils/files/charmm_amoeba.xml +++ b/gmso/utils/files/charmm_amoeba.xml @@ -3,7 +3,7 @@ - + From 2f7657e5fd86de603957bd6215a91313e224a404 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 30 Sep 2024 14:54:07 +0000 Subject: [PATCH 8/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- gmso/formats/prm_writer.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index 6f595db76..d2ceb4378 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -20,7 +20,7 @@ def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() lennard_jones_potential = templates["LennardJonesPotential"] - lennard_jones_potential.expression /= 4 # no 4*epsilon term + lennard_jones_potential.expression /= 4 # no 4*epsilon term harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] ub_angle_potential = templates["UreyBradleyAnglePotential"] @@ -60,7 +60,7 @@ def write_prm(topology, filename, strict_potentials=False): "dihedrals": ["PeriodicTorsionPotential"], "angles": ["LAMMPSHarmonicAnglePotential", "UreyBradleyAnglePotential"], "bonds": ["LAMMPSHarmonicBondPotential"], - "sites":["LennardJonesPotential"], + "sites": ["LennardJonesPotential"], # "sites":"CoulombicPotential" } _try_default_potential_conversions(topology, default_parameterMaps) From 2fab3c0034d31a7de604ea408f393b34c591b652 Mon Sep 17 00:00:00 2001 From: CalCraven Date: Mon, 4 Nov 2024 21:13:10 -0600 Subject: [PATCH 9/9] fix failing tests with singularity templates --- gmso/formats/prm_writer.py | 41 +++++++++++++++++++---- gmso/lib/potential_templates.py | 3 +- gmso/tests/test_potential_templates.py | 4 --- gmso/tests/test_specific_ff_to_residue.py | 1 + gmso/utils/conversions.py | 7 ++-- gmso/utils/specific_ff_to_residue.py | 2 +- 6 files changed, 43 insertions(+), 15 deletions(-) diff --git a/gmso/formats/prm_writer.py b/gmso/formats/prm_writer.py index 6f595db76..42f08059a 100644 --- a/gmso/formats/prm_writer.py +++ b/gmso/formats/prm_writer.py @@ -20,7 +20,7 @@ def _accepted_potentials(): """List of accepted potentials that LAMMPS can support.""" templates = PotentialTemplateLibrary() lennard_jones_potential = templates["LennardJonesPotential"] - lennard_jones_potential.expression /= 4 # no 4*epsilon term + lennard_jones_potential.expression /= 4 harmonic_bond_potential = templates["LAMMPSHarmonicBondPotential"] harmonic_angle_potential = templates["LAMMPSHarmonicAnglePotential"] ub_angle_potential = templates["UreyBradleyAnglePotential"] @@ -54,24 +54,49 @@ def write_prm(topology, filename, strict_potentials=False): # Validation # TODO: Use strict_x, (e.g. x=bonds) to validate what topology attrs to convert if not strict_potentials: - default_parameterMaps = { # TODO: sites are not checked currently because gmso + templates = PotentialTemplateLibrary() + lennard_jones_potential = templates["LennardJonesPotential"] + lennard_jones_potential.expression /= 4 # no 4*epsilon term + default_parameterMaps = { # TODO: sites coulomb potential not checked here # doesn't store pair potential eqn the same way as the connections. "impropers": ["LAMMPSHarmonicImproperPotential"], "dihedrals": ["PeriodicTorsionPotential"], "angles": ["LAMMPSHarmonicAnglePotential", "UreyBradleyAnglePotential"], "bonds": ["LAMMPSHarmonicBondPotential"], - "sites":["LennardJonesPotential"], + "sites": [lennard_jones_potential], # "sites":"CoulombicPotential" } _try_default_potential_conversions(topology, default_parameterMaps) _validate_potential_compatibility(topology) unit_system = LAMMPS_UnitSystems("real") # ang, kcal/mol, amu + + # need a unique string to pair each Atom to the NonBonded force sections + # this uniqueness is akin to PotentialFilters.UNIQUE_PARAMETERS + # uniqueTypesDict = dict() # key is tuple of (sigma, epsilon) + # pfilter = PotentialFilters.UNIQUE_PARAMETERS + # [atype for atype in topology.atom_types(pfilter)] + # ATOMS + # columns: unique name, mass with open(filename, "w") as f: f.write("ATOMS\n") unique_atoms = set() for site in topology.sites: + # key = + # if key not in uniqueTypesDict: + # sigma + # epsilon + # name + # mass + # sigma 1-4 + # epsilon 1-4 + # else: + # vals = uniqueTypesDict[key] + + # unique_atoms.add(vals) + ######## + unique_atoms.add( ( site.atom_type.name, @@ -83,7 +108,7 @@ def write_prm(topology, filename, strict_potentials=False): for atom in unique_atoms: f.write("MASS -1 {:8s} {:8s}\n".format(atom[0], atom[1])) - # TODO: Works for harmonic bonds + # TODO: Works for harmonic bonds only f.write("\nBONDS\n") for bond in topology.bond_types: atom1, atom2 = bond.member_types @@ -180,10 +205,10 @@ def write_prm(topology, filename, strict_potentials=False): ) f.write("\nNONBONDED\n") + # columns: unique name, 0.0, epsilon, rmin/2, 0.0, epsilon 1-4, rmin/2 (1-4) # nonbonded14 = topology.scaling_factors[0][2] for atype in topology.atom_types: - # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4) f.write( "{:8s} {:8.3f} {:8s} {:8s} {:8.3f} {:8s} {:8s}\n".format( atype.name, @@ -211,7 +236,11 @@ def write_prm(topology, filename, strict_potentials=False): def _try_default_potential_conversions(top, potentialsDict): """Take a topology a convert all potentials to the style in potentialDict.""" "" for pot_container in potentialsDict: - containerStr = pot_container[:-1] + "_types" + if pot_container == "sites": # sites contain atom_types + containerStr = "atom_types" + else: + # connections contain "{connection}_types" + containerStr = pot_container[:-1] + "_types" if getattr(top, containerStr): for potential in potentialsDict[pot_container]: top.convert_potential_styles({pot_container: potential}) diff --git a/gmso/lib/potential_templates.py b/gmso/lib/potential_templates.py index 82585e220..cce5a999b 100644 --- a/gmso/lib/potential_templates.py +++ b/gmso/lib/potential_templates.py @@ -16,7 +16,6 @@ UnknownParameterError, ) from gmso.utils.expression import PotentialExpression -from gmso.utils.singleton import Singleton POTENTIAL_JSONS = list(Path(__file__).parent.glob("jsons/*.json")) JSON_DIR = Path.joinpath(Path(__file__).parent, "jsons") @@ -160,7 +159,7 @@ def assert_can_parameterize_with( ) -class PotentialTemplateLibrary(Singleton): +class PotentialTemplateLibrary: """A singleton collection of all the potential templates.""" def __init__(self): diff --git a/gmso/tests/test_potential_templates.py b/gmso/tests/test_potential_templates.py index 1c93747dc..2f990e327 100644 --- a/gmso/tests/test_potential_templates.py +++ b/gmso/tests/test_potential_templates.py @@ -14,10 +14,6 @@ class TestPotentialTemplates(BaseTest): def templates(self): return PotentialTemplateLibrary() - def test_singleton_behavior(self, templates): - assert id(templates) == id(PotentialTemplateLibrary()) - assert id(PotentialTemplateLibrary()) == id(PotentialTemplateLibrary()) - def test_lennard_jones_potential(self, templates): lennard_jones_potential = templates["LennardJonesPotential"] assert lennard_jones_potential.name == "LennardJonesPotential" diff --git a/gmso/tests/test_specific_ff_to_residue.py b/gmso/tests/test_specific_ff_to_residue.py index 3f0361a1c..61637787d 100644 --- a/gmso/tests/test_specific_ff_to_residue.py +++ b/gmso/tests/test_specific_ff_to_residue.py @@ -525,6 +525,7 @@ def test_specific_ff_params_benzene_aa_grouped(self): gmso_match_ff_by="group", boxes_for_simulation=1, ) + assert test_topology.n_sites == 4 assert test_topology.n_bonds == 0 assert test_topology.n_angles == 0 diff --git a/gmso/utils/conversions.py b/gmso/utils/conversions.py index f2a1a4c01..df4ee7bfc 100644 --- a/gmso/utils/conversions.py +++ b/gmso/utils/conversions.py @@ -92,7 +92,8 @@ def _conversion_from_template_name( current_expression, new_potential ) if modified_connection_parametersDict: # try sympy conversions - current_expression.name = new_potential.name + if not conn_typeStr == "atom_types": + current_expression.name = new_potential.name current_expression.expression = new_potential.expression current_expression.parameters.update(modified_connection_parametersDict) @@ -112,7 +113,9 @@ def _conversion_from_template_obj( current_expression, potential_template ) if modified_connection_parametersDict: # try sympy conversions - current_expression.name = potential_template.name + if not conn_typeStr == "atom_types": + # only change name for connections, leave atom_type name + current_expression.name = potential_template.name current_expression.expression = potential_template.expression current_expression.parameters.update(modified_connection_parametersDict) diff --git a/gmso/utils/specific_ff_to_residue.py b/gmso/utils/specific_ff_to_residue.py index 6639b778f..a8ab0d338 100644 --- a/gmso/utils/specific_ff_to_residue.py +++ b/gmso/utils/specific_ff_to_residue.py @@ -238,7 +238,7 @@ def specific_ff_to_residue( {'expression': confirmed singular improper types expression or equation, 'improper_types': gmso Topology.improper_types}. combining_rule: str - The possible mixing/combining rules are 'geometric' or 'lorentz', + The possible mixing/combining rules are 'geometric' or 'lorentz', which provide the geometric and arithmetic mixing rule, respectively. NOTE: Arithmetic means the 'lorentz' combining or mixing rule. NOTE: GMSO default to the 'lorentz' mixing rule if none is provided,