Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Charmm parameters writer #848

Open
wants to merge 17 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
3327de2
initially create charmm parameter writer files and tests
CalCraven Sep 16, 2024
25f893f
function charmm file writers
CalCraven Sep 16, 2024
13e0d6e
Added support and testing files for urey bradley potentials
CalCraven Sep 27, 2024
2c7df9c
Add extra checking for reading independent variables from forcefield …
CalCraven Sep 27, 2024
7bfbb3f
Add compatibility checks for charmm parm files
CalCraven Sep 30, 2024
ea5ec28
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Sep 30, 2024
2f10015
Fix unused variables in charmm writer and update epsilon in lj equati…
CalCraven Sep 30, 2024
362f747
Merge branch 'main' of https://github.com/mosdef-hub/gmso into charmm…
CalCraven Sep 30, 2024
86b6d06
Merge branch 'charmm-writer' of https://github.com/CalCraven/gmso int…
CalCraven Sep 30, 2024
2f7657e
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Sep 30, 2024
9b4c3e7
Merge branch 'main' into charmm-writer
chrisjonesBSU Oct 15, 2024
d107c21
Merge branch 'main' of https://github.com/mosdef-hub/gmso into charmm…
CalCraven Oct 15, 2024
d150a31
Merge branch 'main' into charmm-writer
chrisjonesBSU Oct 18, 2024
2fb6ffe
Merge branch 'main' of https://github.com/mosdef-hub/gmso into charmm…
CalCraven Oct 30, 2024
2fab3c0
fix failing tests with singularity templates
CalCraven Nov 5, 2024
c325111
Merge branch 'charmm-writer' of https://github.com/CalCraven/gmso int…
CalCraven Nov 5, 2024
a3c0d42
Merge branch 'main' into charmm-writer
CalCraven Dec 18, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions gmso/formats/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
250 changes: 250 additions & 0 deletions gmso/formats/prm_writer.py
Original file line number Diff line number Diff line change
@@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to make sure we're using this everywhere we write out parameters? I see we use it for the atom mass, but not anywhere else.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep, for in order to handle any input unyts this will do the filtering. I was just too lazy to do it every time at first. Was wondering if you think we should make a system that isn't specific to LAMMPS for this one. It is the exact same as the real units, but maybe a touch confusing to copy them across to other engines. What do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if unit_system should be a parameter of write_prm that defaults to None. If nothing is passed, then all units are used as-is from the toplogy? If we set unit_system = topology.unit_system as default behavior, would that work as well?

def write_prm(topology, filename, unit_system=None):
    if not unit_system:
        unit_system = topology.unit_system

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The PRM file looks to me like it has default assumed units of kCal/mol, angstrom, amu, radians, etc for it's units. So it can only convert to one defined system.


# 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
Comment on lines +87 to +91

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.
# sigma 1-4
# epsilon 1-4
# else:
# vals = uniqueTypesDict[key]
Comment on lines +94 to +95

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

# 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)}"
)
9 changes: 9 additions & 0 deletions gmso/lib/jsons/LAMMPSHarmonicImproperPotential.json
Original file line number Diff line number Diff line change
@@ -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"
}
}
11 changes: 11 additions & 0 deletions gmso/lib/jsons/UreyBradleyAnglePotential.json
Original file line number Diff line number Diff line change
@@ -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"
}
}
9 changes: 6 additions & 3 deletions gmso/lib/potential_templates.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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):
Expand Down
4 changes: 0 additions & 4 deletions gmso/tests/test_potential_templates.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
62 changes: 62 additions & 0 deletions gmso/tests/test_prm.py
Original file line number Diff line number Diff line change
@@ -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
1 change: 1 addition & 0 deletions gmso/tests/test_specific_ff_to_residue.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading
Loading