diff --git a/launch.json b/launch.json index c12767ba..2f71cc15 100644 --- a/launch.json +++ b/launch.json @@ -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/" } ] } diff --git a/src/kimmdy/constants.py b/src/kimmdy/constants.py index 29af47ba..8a3f4787 100644 --- a/src/kimmdy/constants.py +++ b/src/kimmdy/constants.py @@ -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", diff --git a/src/kimmdy/coordinates.py b/src/kimmdy/coordinates.py index aa78f35a..6b514d27 100644 --- a/src/kimmdy/coordinates.py +++ b/src/kimmdy/coordinates.py @@ -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 @@ -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 @@ -220,7 +218,6 @@ def _get_explicit_MultipleDihedrals( else: dihedrals_in = self.mol_a.proper_dihedrals.get(key) - if not dihedrals_in: return None @@ -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) @@ -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 @@ -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 @@ -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: @@ -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() ) @@ -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"], ) ) @@ -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"], ) )