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 new file mode 100644 index 000000000..bbb47c999 --- /dev/null +++ b/gmso/formats/prm_writer.py @@ -0,0 +1,250 @@ +"""CHARMM Par format.""" + +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): + """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"] + lennard_jones_potential.expression /= 4 # no 4*epsilon term + 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") +def write_prm(topology, filename, strict_potentials=False): + """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" + """ + # Validation + # TODO: Use strict_x, (e.g. x=bonds) to validate what topology attrs to convert + if not strict_potentials: + 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": [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, + unit_system.convert_parameter( + site.atom_type.mass, n_decimals=3, name="mass" + ), + ) + ) + for atom in unique_atoms: + f.write("MASS -1 {:8s} {:8s}\n".format(atom[0], atom[1])) + + # TODO: Works for harmonic bonds only + f.write("\nBONDS\n") + for bond in topology.bond_types: + 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), + unit_system.convert_parameter( + bond.parameters["r_eq"], n_decimals=3 + ), + ) + ) + + f.write("\nANGLES\n") + for angle in topology.angle_types: + atom1, atom2, atom3 = angle.member_types + if angle.name == "UreyBradleyAnglePotential": + f.write( + "{:8s} {:8s} {:8s} {:7s} {:7s} {:7s} {:7s}\n".format( + atom1, + atom2, + atom3, + 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} {:7s} {:7s}\n".format( + atom1, + atom2, + atom3, + unit_system.convert_parameter( + angle.parameters["k"], n_decimals=3 + ), + unit_system.convert_parameter( + angle.parameters["theta_eq"], n_decimals=3, name="theta_eq" + ), + ) + ) + + # These dihedrals need to be PeriodicTorsion Style (Charmm style) + f.write("\nDIHEDRALS\n") + for dihedral in topology.dihedral_types: + # works for PeriodicTorsion + atom1, atom2, atom3, atom4 = dihedral.member_types + 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" + ), + ) + ) + + f.write("\nIMPROPER\n") + for improper in topology.improper_types: + atom1, atom2, atom3, atom4 = improper.member_types + f.write( + "{:8s} {:8s} {:8s} {:8s} {:.5s} {:5.3f} {:.5s}\n".format( + atom1, + atom2, + atom3, + atom4, + unit_system.convert_parameter( + improper.parameters["k"], n_decimals=3 + ), + 0.0, + unit_system.convert_parameter( + improper.parameters["phi_eq"], n_decimals=3 + ), + ) + ) + + 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: + f.write( + "{:8s} {:8.3f} {:8s} {:8s} {:8.3f} {:8s} {:8s}\n".format( + atype.name, + 0.0, # ignored, + unit_system.convert_parameter( + atype.parameters["epsilon"], n_decimals=3 + ), + unit_system.convert_parameter( + atype.parameters["sigma"], n_decimals=3 + ), + 0.0, + unit_system.convert_parameter( + atype.parameters["epsilon"], n_decimals=3 + ), + unit_system.convert_parameter( + atype.parameters["sigma"], n_decimals=3 + ), + ) + ) + # TODO: Add NONBONDED section: NONBONDED nbxmod 5 atom cdiel shift vatom vdistance vswitch - + # !adm jr., 5/08/91, suggested cutoff scheme + 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: + 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}) + 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 new file mode 100644 index 000000000..ab0246b7a --- /dev/null +++ b/gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json @@ -0,0 +1,9 @@ +{ + "name": "LAMMPSHarmonicImproperPotential", + "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 new file mode 100644 index 000000000..55c896913 --- /dev/null +++ b/gmso/lib/jsons/UreyBradleyAnglePotential.json @@ -0,0 +1,11 @@ +{ + "name": "UreyBradleyAnglePotential", + "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/lib/potential_templates.py b/gmso/lib/potential_templates.py index 9a029b59a..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") @@ -70,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: @@ -156,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_prm.py b/gmso/tests/test_prm.py new file mode 100644 index 000000000..a45aaca72 --- /dev/null +++ b/gmso/tests/test_prm.py @@ -0,0 +1,62 @@ +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 + +# 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 + + +class TestPar(BaseTest): + def test_save_charmm(self): + top = Topology.load(get_fn("charmm_dihedral.mol2")) + for site in top.sites: + site.name = f"_{site.name}" + + ff = ForceField(get_fn("charmm_truncated.xml")) + 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.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.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 + # 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 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/files/charmm_amoeba.xml b/gmso/utils/files/charmm_amoeba.xml new file mode 100644 index 000000000..c3d464e67 --- /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_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_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 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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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,