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 1 commit
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
Next Next commit
initially create charmm parameter writer files and tests
  • Loading branch information
CalCraven committed Sep 16, 2024
commit 3327de217947493af7033168a8f6dda3330c985d
181 changes: 181 additions & 0 deletions gmso/formats/prm_writer.py
Original file line number Diff line number Diff line change
@@ -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")
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": "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"
}
}
49 changes: 49 additions & 0 deletions gmso/tests/test_par.py
Original file line number Diff line number Diff line change
@@ -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