Skip to content

Commit

Permalink
topology transistion should be accurate now
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmbuhr
jmbuhr committed Nov 15, 2024
1 parent 5d8d6c6 commit abe31f7
Show file tree
Hide file tree
Showing 3 changed files with 79 additions and 96 deletions.
7 changes: 7 additions & 0 deletions launch.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,13 @@
"request": "launch",
"program": "${workspaceFolder}/src/kimmdy/cmd.py",
"cwd": "${workspaceFolder}/example/alanine_hat_naive/"
},
{
"name": "KIMMDY custom run",
"type": "python",
"request": "launch",
"program": "${workspaceFolder}/src/kimmdy/cmd.py",
"cwd": "${workspaceFolder}/../../examples/gly-hydrolysis/"
}
]
}
14 changes: 14 additions & 0 deletions src/kimmdy/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,20 @@
"dihedral_restraints": [0, 1, 2, 3],
}

# from gmx data dir
DEFAULT_EDISSOC: dict[tuple[str, str], float] = { # type: ignore
tuple(sorted(list(k))): v for k, v in {
("C", "N"): 500.0,
("CA", "C"): 341.0,
("CA", "N"): 379.0,
("CA", "CB"): 400.0,
("CB", "CG"): 400.0,
("CG", "CD"): 400.0,
("CD", "CE"): 400.0,
("CE", "NZ"): 400.0
}.items()
}

# see https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html
FFFUNC = {
"mult_proper_dihedral": "9",
Expand Down
154 changes: 58 additions & 96 deletions src/kimmdy/coordinates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
import MDAnalysis as mda
import numpy as np

from kimmdy.constants import REACTIVE_MOLECULEYPE, FFFUNC
from kimmdy.constants import DEFAULT_EDISSOC, REACTIVE_MOLECULEYPE, FFFUNC
from kimmdy.parsing import read_plumed, write_plumed
from kimmdy.recipe import Place
from kimmdy.tasks import TaskFiles
Expand Down Expand Up @@ -200,8 +200,6 @@ def _get_explicit_MultipleDihedrals(
"""Takes a valid dihedral key and returns explicit
dihedral parameters for a given topology
"""


if use_improper:
funct = FFFUNC["mult_improper_dihedral"]
dihedraltypes =self.ff.improper_dihedraltypes
Expand All @@ -220,7 +218,6 @@ def _get_explicit_MultipleDihedrals(
else:
dihedrals_in = self.mol_a.proper_dihedrals.get(key)


if not dihedrals_in:
return None

Expand All @@ -239,14 +236,12 @@ def _get_explicit_MultipleDihedrals(

type_key.append(t)



multiple_dihedrals = MultipleDihedrals(
*key, funct=funct, dihedrals={}
*key, funct=funct, dihedrals={}
)
for periodicity in range(1, periodicity_max+1):
match_obj = match_atomic_item_to_atomic_type(
type_key, dihedraltypes, str(periodicity)
id=type_key, types=dihedraltypes, periodicity=str(periodicity)
)
if match_obj:
assert isinstance(match_obj, DihedralType)
Expand Down Expand Up @@ -338,103 +333,64 @@ def _make_pair(self, id1: str, id2: str, is_bind: bool) -> Pair:

def _interpolate_two_dihedrals(
self,
dihedral_key: tuple[str, str, str, str],
key: tuple[str, str, str, str],
dihedral_a: Optional[Dihedral],
dihedral_b: Optional[Dihedral],
funct: str,
periodicity: str,
) -> Dihedral:
"""Merge one to two Dihedrals or -Types into a Dihedral in free-energy syntax"""
"""Merge one to two Dihedrals into a Dihedral in free-energy syntax.
if funct == FFFUNC["mult_proper_dihedral"]:
dihedral_types = self.ff.proper_dihedraltypes
elif funct == FFFUNC["mult_improper_dihedral"]:
dihedral_types = self.ff.improper_dihedraltypes
else:
Only one of the dihedrals can be None at any given time.
"""
if funct not in [FFFUNC["mult_proper_dihedral"], FFFUNC["mult_improper_dihedral"]]:
m = f"Can only interpolate between proper (type {FFFUNC['mult_proper_dihedral']}) or improper (type {FFFUNC['mult_proper_dihedral']}) dihedrals"
logger.error(m)
raise ValueError(m)

# convert implicit standard ff parameters to explicit, if necessary

if dihedral_a:
parameterizedA = get_explicit_or_type(
dihedral_key,
dihedral_a,
dihedral_types,
self.mol_a,
periodicity,
)
else:
parameterizedA = None

if dihedral_b:
parameterizedB = get_explicit_or_type(
dihedral_key,
dihedral_b,
dihedral_types,
self.mol_b,
periodicity,
)
else:
parameterizedB = None

# construct parameterized Dihedral
if parameterizedA is not None and parameterizedB is not None:
# same

assert (
type(parameterizedA) == Dihedral or type(parameterizedA) == DihedralType
)
assert (
type(parameterizedB) == Dihedral or type(parameterizedB) == DihedralType
)

if dihedral_a is not None and dihedral_b is not None:
# both parameters are given
# transition between parameters
# this can be the case when both dihedrals exist
# but the atomtypes of participating atoms change from A to B state
dihedralmerge = Dihedral(
*dihedral_key,
*key,
funct=funct,
c0=parameterizedA.c0,
c1=parameterizedA.c1,
periodicity=parameterizedA.periodicity,
c3=parameterizedB.c0,
c4=parameterizedB.c1,
c5=parameterizedB.periodicity,
c0=dihedral_a.c0,
c1=dihedral_a.c1,
periodicity=dihedral_a.periodicity,
c3=dihedral_b.c0,
c4=dihedral_b.c1,
c5=dihedral_b.periodicity,
)
elif parameterizedA is not None:
# breaking
if not (
type(parameterizedA) == Dihedral or type(parameterizedA) == DihedralType
):
m = f"parameterizedA {parameterizedA} is not a Dihedral or DihedralType"
logger.warning(m)
raise ValueError(m)
elif dihedral_a is not None:
# dihedral only in A
# vanish it by transitiononing to 0
dihedralmerge = Dihedral(
*dihedral_key,
*key,
funct=funct,
c0=parameterizedA.c0,
c1=parameterizedA.c1,
periodicity=parameterizedA.periodicity,
c3=parameterizedA.c0,
c0=dihedral_a.c0,
c1=dihedral_a.c1,
periodicity=dihedral_a.periodicity,
c3=dihedral_a.c0,
c4="0.00",
c5=parameterizedA.periodicity,
)
elif parameterizedB is not None:
# binding
assert (
type(parameterizedB) == Dihedral or type(parameterizedB) == DihedralType
c5=dihedral_a.periodicity,
)
elif dihedral_b is not None:
# dihedral only in B
# create it by transitioning from 0
dihedralmerge = Dihedral(
*dihedral_key,
*key,
funct=funct,
c0=parameterizedB.c0,
c0=dihedral_b.c0,
c1="0.00",
periodicity=parameterizedB.periodicity,
c3=parameterizedB.c0,
c4=parameterizedB.c1,
c5=parameterizedB.periodicity,
periodicity=dihedral_b.periodicity,
c3=dihedral_b.c0,
c4=dihedral_b.c1,
c5=dihedral_b.periodicity,
)
else:
m = f"Tried to merge two dihedrals of {dihedral_key} but no parameterized dihedrals found!"
m = f"Tried to merge two dihedrals of {key} but no parameterized dihedrals found!"
logger.error(m)
raise ValueError(m)
return dihedralmerge
Expand Down Expand Up @@ -515,11 +471,16 @@ def merge_bonds(self):
self.break_bind_atoms[key[1]] = atomtypes[1]
sigmaij_b, epsilonij_b = self._get_LJ_parameters(*key, use_state_b=True)

# TODO: get accurate well depth from gmx edissoc file
# get accurate well depth from gmx edissoc file
well_depth = self.default_morse_well_depth
type_key: tuple[str, str] = tuple(sorted(atomtypes)) # type: ignore
d = DEFAULT_EDISSOC.get(type_key)
if d is not None:
well_depth = d

k = parameterizedA.c1
if k is not None:
# like in gmx /src/gromacs/gmxpreprocess/tomorse.cpp
beta = sqrt(float(k)/2*well_depth)
else:
beta = self.default_beta_for_lj
Expand Down Expand Up @@ -574,8 +535,12 @@ def merge_bonds(self):
self.break_bind_atoms[key[1]] = atomtypes[1]
sigmaij_a, epsilonij_a = self._get_LJ_parameters(*key, use_state_b=True)

# TODO: get accurate well depth from gmx edissoc file
# get accurate well depth from gmx edissoc file
well_depth = self.default_morse_well_depth
type_key: tuple[str, str] = tuple(sorted(atomtypes)) # type: ignore
d = DEFAULT_EDISSOC.get(type_key)
if d is not None:
well_depth = d

k = parameterizedB.c1
if k is not None:
Expand Down Expand Up @@ -701,7 +666,6 @@ def merge_dihedrals(self):

# proper dihedrals
# proper dihedrals have a nested structure and need a different treatment from bonds, angles and improper dihedrals
# if indices change atomtypes and parameters change because of that, it will ignore these parameter change
keys = set(self.mol_a.proper_dihedrals.keys()) | set(
self.mol_b.proper_dihedrals.keys()
)
Expand Down Expand Up @@ -746,11 +710,10 @@ def merge_dihedrals(self):

self.mol_b.proper_dihedrals[key].dihedrals[periodicity] = (
self._interpolate_two_dihedrals(
key,
interactionA,
interactionB,
FFFUNC["mult_proper_dihedral"],
periodicity,
key=key,
dihedral_a=interactionA,
dihedral_b=interactionB,
funct=FFFUNC["mult_proper_dihedral"],
)
)

Expand Down Expand Up @@ -795,11 +758,10 @@ def merge_dihedrals(self):

self.mol_b.improper_dihedrals[key].dihedrals[periodicity] = (
self._interpolate_two_dihedrals(
key,
interactionA,
interactionB,
FFFUNC["mult_improper_dihedral"],
periodicity,
key=key,
dihedral_a=interactionA,
dihedral_b=interactionB,
funct=FFFUNC["mult_improper_dihedral"],
)
)

Expand Down

0 comments on commit abe31f7

Please sign in to comment.