from hoobas import LinearChain
from hoobas.Units import SimulationUnits
import warnings
from copy import deepcopy
from copy import copy
from hoobas.CGBead import Bead as CGBead
import numpy as np
import math
from hoobas import Composite
from collections import OrderedDict
from hoobas import Util
import re
from typing import Union
try:
NETWORKX_FOUND=True
import networkx
except ImportError:
NETWORKX_FOUND=False
print('Networkx package not found, itp autobuilder will be disabled')
[docs]class MMaps(type):
"""
Class that allows for global remapping of martini types in case they are needed. By default, martini bead names
are charged, polar, nonpolar and apolar. By calling MartiniMapping.remap, they can be changed. For instance,
MartiniMappings.remap('charged_0', 'Q00') will make all charged beads with sub 0 to have the name 'Q00'. This is useful
for conflicting names like 'Na' which may conflict with sodium ions
"""
# default Martini mapping
__MARTINI_MAP = {
'charged_0': 'Q0',
'charged_a': 'Qa',
'charged_d': 'Qd',
'charged_da': 'Qda',
'polar_1': 'P1',
'polar_2': 'P2',
'polar_3': 'P3',
'polar_4': 'P4',
'polar_5': 'P5',
'nonpolar_0': 'N0',
'nonpolar_a': 'Na',
'nonpolar_d': 'Nd',
'nonpolar_da': 'Nda',
'apolar_1': 'C1',
'apolar_2': 'C2',
'apolar_3': 'C3',
'apolar_4': 'C4',
'apolar_5': 'C5',
'big_polar_4': 'BP4',
'ring_charged_0': 'SQ0',
'ring_charged_a': 'SQa',
'ring_charged_d': 'SQd',
'ring_charged_da': 'SQda',
'ring_polar_1': 'SP1',
'ring_polar_2': 'SP2',
'ring_polar_3': 'SP3',
'ring_polar_4': 'SP4',
'ring_polar_5': 'SP5',
'ring_nonpolar_0': 'SN0',
'ring_nonpolar_a': 'SNa',
'ring_nonpolar_d': 'SNd',
'ring_nonpolar_da': 'SNda',
'ring_apolar_1': 'SC1',
'ring_apolar_2': 'SC2',
'ring_apolar_3': 'SC3',
'ring_apolar_4': 'SC4',
'ring_apolar_5': 'SC5',
'ring_big_polar_4': 'SBP4',
'polarizable': 'POL',
'dummy': 'D'
}
# inverse lookup table
__INVERSE_MARTINI_MAP = {
'P2': 'polar_2',
'P3': 'polar_3',
'P1': 'polar_1',
'P4': 'polar_4',
'P5': 'polar_5',
'Nda': 'nonpolar_da',
'Na': 'nonpolar_a',
'Qda': 'charged_da',
'Nd': 'nonpolar_d',
'Q0': 'charged_0',
'Qa': 'charged_a',
'N0': 'nonpolar_0',
'Qd': 'charged_d',
'C3': 'apolar_3',
'C2': 'apolar_2',
'C1': 'apolar_1',
'C5': 'apolar_5',
'C4': 'apolar_4',
'BP4': 'big_polar_4',
'SP2': 'ring_polar_2',
'SP3': 'ring_polar_3',
'SP1': 'ring_polar_1',
'SP4': 'ring_polar_4',
'SP5': 'ring_polar_5',
'SNda': 'ring_nonpolar_da',
'SNa': 'ring_nonpolar_a',
'SQda': 'ring_charged_da',
'SNd': 'ring_nonpolar_d',
'SQ0': 'ring_charged_0',
'SQa': 'ring_charged_a',
'SN0': 'ring_nonpolar_0',
'SQd': 'ring_charged_d',
'SC3': 'ring_apolar_3',
'SC2': 'ring_apolar_2',
'SC1': 'ring_apolar_1',
'SC5': 'ring_apolar_5',
'SC4': 'ring_apolar_4',
'SBP4': 'ring_big_polar_4',
'POL': 'polarizable',
'D': 'dummy'
}
# create a default martini map to import external molecules so that we can remap to hoomd by calling
# map[default[item]]
__DEFAULT_MARTINI_MAP = {
'P2': 'polar_2',
'P3': 'polar_3',
'P1': 'polar_1',
'P4': 'polar_4',
'P5': 'polar_5',
'Nda': 'nonpolar_da',
'Na': 'nonpolar_a',
'Qda': 'charged_da',
'Nd': 'nonpolar_d',
'Q0': 'charged_0',
'Qa': 'charged_a',
'N0': 'nonpolar_0',
'Qd': 'charged_d',
'C3': 'apolar_3',
'C2': 'apolar_2',
'C1': 'apolar_1',
'C5': 'apolar_5',
'C4': 'apolar_4',
'BP4': 'big_polar_4',
'SP2': 'ring_polar_2',
'SP3': 'ring_polar_3',
'SP1': 'ring_polar_1',
'SP4': 'ring_polar_4',
'SP5': 'ring_polar_5',
'SNda': 'ring_nonpolar_da',
'SNa': 'ring_nonpolar_a',
'SQda': 'ring_charged_da',
'SNd': 'ring_nonpolar_d',
'SQ0': 'ring_charged_0',
'SQa': 'ring_charged_a',
'SN0': 'ring_nonpolar_0',
'SQd': 'ring_charged_d',
'SC3': 'ring_apolar_3',
'SC2': 'ring_apolar_2',
'SC1': 'ring_apolar_1',
'SC5': 'ring_apolar_5',
'SC4': 'ring_apolar_4',
'SBP4': 'ring_big_polar_4',
'POL': 'polarizable',
'D': 'dummy'
}
[docs] @staticmethod
def remap(key: str, beadtype: str) -> None:
"""
Remaps a martini type to a new beadtype
:param key: martini type to remap
:type key: str
:param beadtype: beadtype to map to
:type beadtype: str
:return: None
"""
if not isinstance(beadtype, str):
warnings.warn('Non-string beadtype inserted in Martini mappings, coercing to string; check results')
old_key = MMaps.__MARTINI_MAP[key]
MMaps.__MARTINI_MAP[key] = beadtype
MMaps.__INVERSE_MARTINI_MAP[beadtype] = MMaps.__INVERSE_MARTINI_MAP.pop(old_key)
[docs] @staticmethod
def map() -> dict:
"""
Returns the Martini mapping
:return: map
:rtype: dict
"""
return MMaps.__MARTINI_MAP
def __getitem__(self, beadtype: str) -> str:
return MMaps.__MARTINI_MAP[beadtype]
[docs] @staticmethod
def reverse_lookup(beadname: str) -> str:
"""
Search which martini type a given beadtype corresponds to
:param beadname: bead to search
:type beadname: str
:return: Martini type
:rtype: str
"""
return MMaps.__INVERSE_MARTINI_MAP[beadname]
[docs] @staticmethod
def map_from_default(default_name: str) -> str:
"""
Returns the current beadtype of a default martini beadtype. For instance if P2 was remapped to Po2, then
map_from_default(P2) will yield Po2
:param default_name: default Martini mapping
:type default_name: str
:return: mapped beadtype
:rtype: str
"""
return MMaps.__MARTINI_MAP[MMaps.__DEFAULT_MARTINI_MAP[default_name]]
[docs] @staticmethod
def type_from_default(default_name: str) -> str:
"""
Returns the martini type of a given default beadtype
:param default_name: beadtype
:type default_name: str
:return: martini type
:rtype: str
"""
return MMaps.__DEFAULT_MARTINI_MAP[default_name]
[docs]class MartiniMappings(object, metaclass=MMaps):
"""
Some syntaxic sugar for python to make it able to call MartiniMappings[key] directly
"""
def __init__(self):
warnings.warn('Building a Mappings class, this should not happen, check results')
[docs]class MartiniForceMappings(object):
"""
Martini force mapping; use the get_pair(beadA, beadB, units) to get the LJ potential.
"""
__FORCE_TYPES = {
'O': {'eps': 5.6, 'sigma': 4.7},
'I': {'eps': 5.0, 'sigma': 4.7},
'II': {'eps': 4.5, 'sigma': 4.7},
'III': {'eps': 4.0, 'sigma': 4.7},
'IV': {'eps': 3.5, 'sigma': 4.7},
'V': {'eps': 3.1, 'sigma': 4.7},
'VI': {'eps': 2.7, 'sigma': 4.7},
'VII': {'eps': 2.3, 'sigma': 4.7},
'VIII': {'eps': 2.0, 'sigma': 4.7},
'IX': {'eps': 2.0, 'sigma': 6.2},
'PO': {'eps': 5.6*0.95, 'sigma': 4.7},
'PI': {'eps': 5.0*0.95, 'sigma': 4.7},
'PII': {'eps': 4.5*0.95, 'sigma': 4.7},
'PIII': {'eps': 4.0*0.95, 'sigma': 4.7},
'PIV': {'eps': 3.5*0.95, 'sigma': 4.7},
'PV': {'eps': 3.1*0.95, 'sigma': 4.7},
'PVI': {'eps': 2.7*0.95, 'sigma': 4.7},
'PVII': {'eps': 2.3*0.95, 'sigma': 4.7},
'PVIII': {'eps': 2.0*0.95, 'sigma': 4.7},
'PIX': {'eps': 2.0*0.95, 'sigma': 6.2},
'BP4-P4': {'eps': 5.6, 'sigma': 5.7}
}
__PAIR_TABLES = [ # last column is big polar 4
['O', 'O', 'O', 'II', 'O', 'O', 'O', 'I', 'I', 'I', 'I', 'I', 'IV', 'V', 'VI', 'VII', 'IX', 'IX', 'O'],
['O', 'I', 'O', 'II', 'O', 'O', 'O', 'I', 'I', 'I', 'III', 'I', 'IV', 'V', 'VI', 'VII', 'IX', 'IX', 'O'],
['O', 'O', 'I', 'II', 'O', 'O', 'O', 'I', 'I', 'I', 'I', 'III', 'IV', 'V', 'VI', 'VII', 'IX', 'IX', 'O'],
['II', 'II', 'II', 'IV', 'I', 'O', 'I', 'II', 'III', 'III', 'III', 'III', 'IV', 'V', 'VI', 'VII', 'IX', 'IX', 'O'],
['O', 'O', 'O', 'I', 'O', 'O', 'O', 'O', 'O', 'I', 'I', 'I', 'IV', 'V', 'VI', 'VI', 'VII', 'VIII', 'O'],
['O', 'O', 'O', 'O', 'O', 'I', 'I', 'II', 'II', 'III', 'III', 'III', 'IV', 'V', 'VI', 'VI', 'VII', 'VIII', 'BP4-P4'],
['O', 'O', 'O', 'I', 'O', 'I', 'I', 'II', 'II', 'II', 'II', 'II', 'IV', 'IV', 'V', 'V', 'VI', 'VII', 'I'],
['I', 'I', 'I', 'II', 'O', 'II', 'II', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'V', 'VI', 'VII', 'II'],
['I', 'I', 'I', 'III', 'O', 'II', 'II', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'IV', 'V', 'VI', 'II'],
['I', 'I', 'I', 'III', 'I', 'III', 'II', 'II', 'II', 'II', 'II', 'II', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'III'],
['I', 'III', 'I', 'III', 'I', 'III', 'II', 'II', 'II', 'II', 'III', 'II', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'III'],
['I', 'I', 'III', 'III', 'I', 'III', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'III'],
['IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'III', 'III', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'V', 'VI', 'IV'],
['V', 'V', 'V', 'V', 'V', 'V', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'V', 'V', 'V'],
['VI', 'VI', 'VI', 'VI', 'VI', 'VI', 'V', 'IV', 'IV', 'V', 'V', 'V', 'IV', 'IV', 'IV', 'IV', 'V', 'V', 'VI'],
['VII', 'VII', 'VII', 'VII', 'VI', 'VI', 'V', 'V', 'IV', 'VI', 'VI', 'VI', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'VI'],
['IX', 'IX', 'IX', 'IX', 'VII', 'VII', 'VI', 'VI', 'V', 'VI', 'VI', 'VI', 'V', 'V', 'V', 'IV', 'IV', 'IV', 'VII'],
['IX', 'IX', 'IX', 'IX', 'VIII', 'VIII', 'VII', 'VII', 'VI', 'VI', 'VI', 'VI', 'VI', 'V', 'V', 'IV', 'IV', 'IV', 'VIII'],
['O', 'O', 'O', 'O', 'O', 'BP4-P4', 'I', 'II', 'II', 'III', 'III', 'III', 'IV', 'V', 'VI', 'VI', 'VII', 'VIII', 'I'] # big polar4 line
]
__PAIR_POL_TABLES = [
['O', 'I', 'I', 'IV', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'III', 'IV', 'V', 'VI', 'VII', 'VII', 'O'],
['I', 'IV', 'III', 'VII', 'O', 'O', 'O', 'O', 'O', 'O', 'II', 'O', 'III', 'IV', 'V', 'VI', 'VII', 'VII', 'I'],
['I', 'III', 'IV', 'VII', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VII', 'I'],
['IV', 'VII', 'VII', 'IV', 'O', 'O', 'O', 'I', 'II', 'II', 'II', 'II', 'III', 'IV', 'V', 'VI', 'VII', 'VII', 'II'],
['O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'I', 'I', 'I', 'IV', 'V', 'VI', 'VI', 'VII', 'VIII', 'PO'],
['O', 'O', 'O', 'O', 'O', 'I', 'I', 'II', 'II', 'III', 'III', 'III', 'IV', 'V', 'VI', 'VI', 'VII', 'VIII', 'PI'],
['O', 'O', 'O', 'O', 'O', 'I', 'I', 'II', 'II', 'II', 'II', 'II', 'IV', 'IV', 'V', 'V', 'VI', 'VII', 'PI'],
['O', 'O', 'O', 'I', 'O', 'II', 'II', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'V', 'VI', 'VII', 'PII'],
['O', 'O', 'O', 'II', 'O', 'II', 'II', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'IV', 'V', 'VI', 'PII'],
['O', 'O', 'O', 'II', 'I', 'III', 'II', 'II', 'II', 'II', 'II', 'II', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'PIII'],
['O', 'II', 'O', 'II', 'I', 'III', 'II', 'II', 'II', 'II', 'III', 'II', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'PIII'],
['O', 'O', 'II', 'II', 'I', 'III', 'II', 'II', 'II', 'II', 'II', 'III', 'IV', 'IV', 'V', 'VI', 'VI', 'VI', 'PIII'],
['III', 'III', 'III', 'III', 'IV', 'IV', 'IV', 'III', 'III', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'V', 'VI', 'PIV'],
['IV', 'IV', 'IV', 'IV', 'V', 'V', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'V', 'V', 'PV'],
['V', 'V', 'V', 'V', 'VI', 'VI', 'V', 'IV', 'IV', 'V', 'V', 'V', 'IV', 'IV', 'IV', 'IV', 'V', 'V', 'PVI'],
['VI', 'VI', 'VI', 'VI', 'VI', 'VI', 'V', 'V', 'IV', 'VI', 'VI', 'VI', 'IV', 'IV', 'IV', 'IV', 'IV', 'IV', 'PVI'],
['VII', 'VII', 'VII', 'VII', 'VII', 'VII', 'VI', 'VI', 'V', 'VI', 'VI', 'VI', 'V', 'V', 'V', 'IV', 'IV', 'IV', 'PVII'],
['VII', 'VII', 'VII', 'VII', 'VIII', 'VIII', 'VII', 'VII', 'VI', 'VI', 'VI', 'VI', 'VI', 'V', 'V', 'IV', 'IV', 'IV', 'PVIII'],
['O', 'I', 'I', 'II', 'PO', 'PI', 'PI', 'PII', 'PII', 'PIII', 'PIII', 'PIII', 'PIV', 'PV', 'PVI', 'PVI', 'PVII', 'PVIII', 'III']
]
__TYPE_TO_INDEX = {
'charged_da': 0,
'charged_d': 1,
'charged_a': 2,
'charged_0': 3,
'polar_5': 4,
'polar_4': 5,
'polar_3': 6,
'polar_2': 7,
'polar_1': 8,
'nonpolar_da': 9,
'nonpolar_d': 10,
'nonpolar_a': 11,
'nonpolar_0': 12,
'apolar_5': 13,
'apolar_4': 14,
'apolar_3': 15,
'apolar_2': 16,
'apolar_1': 17,
'big_polar_4': 18,
'polarizable': 18,
'ring_charged_da': 0,
'ring_charged_d': 1,
'ring_charged_a': 2,
'ring_charged_0': 3,
'ring_polar_5': 4,
'ring_polar_4': 5,
'ring_polar_3': 6,
'ring_polar_2': 7,
'ring_polar_1': 8,
'ring_nonpolar_da': 9,
'ring_nonpolar_d': 10,
'ring_nonpolar_a': 11,
'ring_nonpolar_0': 12,
'ring_apolar_5': 13,
'ring_apolar_4': 14,
'ring_apolar_3': 15,
'ring_apolar_2': 16,
'ring_apolar_1': 17,
'ring_big_polar_4': 18,
'ring_polarizable': 18
}
def __init__(self):
warnings.warn('This class should not be instanciated')
[docs] @staticmethod
def get_pair(beadtypeA: str, beadtypeB: str, units: Union[None, SimulationUnits] = None): # units defaults to Ang; kJ / mol and amu
"""
Returns the pair ineraction parameters between two beadtypes. Units can be supplied, otherwise it uses whichever
units were defined as default
:param beadtypeA: first beadtype
:type beadtypeA: str
:param beadtypeB: second beadtype
:type beadtypeB: str
:param units: units of the LJ parameters
:type units: SimulationUnits
:return: LJ parameters
:rtype: dict
"""
# get initial name for the bead
martini_beadnameA = MartiniMappings.reverse_lookup(beadtypeA)
martini_beadnameB = MartiniMappings.reverse_lookup(beadtypeB)
rA = 'ring' in martini_beadnameA
rB = 'ring' in martini_beadnameB
epsScale = 0.75 if rA and rB else 1.0
sigmaScale = 4.3/4.7 if rA and rB else 1.0
beadIndex1 = MartiniForceMappings.__TYPE_TO_INDEX[martini_beadnameA]
beadIndex2 = MartiniForceMappings.__TYPE_TO_INDEX[martini_beadnameB]
if beadIndex1 > beadIndex2:
beadIndex1, beadIndex2 = beadIndex2, beadIndex1
key_interaction = MartiniForceMappings.__PAIR_TABLES[beadIndex1][beadIndex2]
forces = MartiniForceMappings.__FORCE_TYPES[key_interaction]
if units is None:
units = SimulationUnits()
return {'eps': epsScale * units.get_energy(forces['eps'], 'kJ/mol'), 'sigma': sigmaScale* units.get_length(forces['sigma'], 'A')}
[docs] @staticmethod
def get_pair_polarizable(beadtypeA: str, beadtypeB: str, units: Union[None, SimulationUnits] = None):
"""
Returns the pair ineraction parameters between two beadtypes for the polarizable martini force-field. Units can
be supplied, otherwise it uses whichever units were defined as default
:param beadtypeA: first beadtype
:type beadtypeA: str
:param beadtypeB: second beadtype
:type beadtypeB: str
:param units: units of the LJ parameters
:type units: SimulationUnits
:return: LJ parameters
:rtype: dict
"""
# get initial name for the bead
martini_beadnameA = MartiniMappings.reverse_lookup(beadtypeA)
martini_beadnameB = MartiniMappings.reverse_lookup(beadtypeB)
if martini_beadnameA == 'dummy' or martini_beadnameB == 'dummy':
return {'eps': 0.0, 'sigma': 0.0}
rA = 'ring' in martini_beadnameA
rB = 'ring' in martini_beadnameB
epsScale = 0.75 if rA and rB else 1.0
sigmaScale = 4.3/4.7 if rA and rB else 1.0
beadIndex1 = MartiniForceMappings.__TYPE_TO_INDEX[martini_beadnameA]
beadIndex2 = MartiniForceMappings.__TYPE_TO_INDEX[martini_beadnameB]
if beadIndex1 > beadIndex2:
beadIndex1, beadIndex2 = beadIndex2, beadIndex1
key_interaction = MartiniForceMappings.__PAIR_POL_TABLES[beadIndex1][beadIndex2]
forces = MartiniForceMappings.__FORCE_TYPES[key_interaction]
if units is None:
units = SimulationUnits()
return {'eps': epsScale * units.get_energy(forces['eps'], 'kJ/mol'), 'sigma': sigmaScale* units.get_length(forces['sigma'], 'A')}
[docs]class MartiniChain(LinearChain.LinearChain):
"""
General implementation of most Martini chains (lipids or others). Sets the bond length, units and Martini mappings
"""
# indexes to generate different angles in hoomd
_MARTINI_BOND_INDEX = 0
_MARTINI_ANGLE_INDEX = 0
_MARTINI_DIHEDRAL_INDEX = 0
def __init__(self, **kwargs):
# default all martini to nm, kJ/mol and amu
_su = SimulationUnits()
_su.set_length('nm')
_su.set_energy('kJ/mol')
_su.set_mass('amu')
# kuhn_length set to 0.47 since this is the martini value
super(MartiniChain, self).__init__(kuhn_length=0.47, units=_su, **kwargs)
self.mapped = False
def __map(self) -> None:
if self.mapped:
warnings.warn('mapping method was called on object that was already mapped; check beadtypes; martini mapping'
'skipped for this object of class: '+str(self.__class__.__name__), UserWarning)
return
for bead in self.beads:
try:
bead.type = MartiniMappings[bead.beadtype]
except KeyError:
warnings.warn('Found with type that doesn''t match Martini types during remap, beadtype was : '
+ str(bead.beadtype) + '; check beadtypes', UserWarning)
self.mapped = True
def __deepcopy__(self, memodict): # override the deepcopy method so that types are mapped when added to the system
cls = self.__class__
cp = cls.__new__(cls)
memodict[id(self)] = cp
for k, v in self.__dict__.items():
setattr(cp, k, deepcopy(v, memodict))
if not cp.mapped:
cp.__map()
return cp
[docs]class MartiniPolarizableWaterMolecule(MartiniChain):
def __init__(self, **kwargs):
super(MartiniPolarizableWaterMolecule, self).__init__(**kwargs)
self.beads += CGBead(position=np.array([0., 0., -0.14]), beadtype='D', charge=-0.46, mass=72.0)
self.beads += CGBead(position=np.array([0., 0., 0.0]), beadtype='POL', charge=0.0, mass=72.0)
self.beads += CGBead(position=np.array([0., 0., 0.14]), beadtype='D', charge=0.46, mass=72.0)
self.p_num.append(0)
self.p_num.append(1)
self.p_num.append(2)
self.a_types += ['MartiniWater', 4.2, 0.0]
self.angles += ['MartiniWater', 0, 1, 2]
self.constraints += Composite.Constraint((0, 1), 0.14)
self.constraints += Composite.Constraint((2, 1), 0.14)
self.bonds += ['Excl', 0, 2]
[docs]class ParsingError(Exception):
pass
[docs]class ITPMartiniParser(MartiniChain): # TODO: add a parsing units options for gro & itp that defaults to nm
"""
Generate a topology from an ITP file
"""
def __init__(self, itp_path, gro_path=None, **kwargs):
super(ITPMartiniParser, self).__init__(**kwargs)
# parse the itp file
itp_file = open(itp_path, 'r')
self.INSANE_field = {}
self.gro_section = {}
self.internal_gro_names = []
current_section = ''
current_headers = []
headers = False
c_dict = OrderedDict()
try:
for line in itp_file.readlines():
if line == '\n' or line[0] == '#':
continue
line = line.strip('\n')
# top of file should be header, read until we hit either a [ character or we find @INSANE
if '@INSANE' in line and 'alhead' in line: # parse the INSANE building information and grab the headgroup
try:
stripped = re.split(',', line.strip(";@INSANE"))
for subitem in stripped:
subitem = subitem.split('=')
self.INSANE_field[subitem[0].strip(' ')] = subitem[1].lstrip().rstrip()
except:
warnings.warn('failed to parse INSANE topology hints')
continue
# remove left leading whitespace
line = line.lstrip(' ')
if line == '':
continue
# remove end of line comments
if line[0] != ';' and ';' in line:
line = line.split(';')
line = line[0]
line = line.strip(' ')
# found section head
if '[' in line and line[0] != ';':
if len(c_dict) > 0:
self.gro_section[current_section] = c_dict
current_section = line.strip('[').strip(']').strip(' ')
headers = True
if 'virtual_sites' in current_section:
current_section = 'virtual_sites'
continue
if current_section == 'exclusions':
line = re.split(' +', line.lstrip(' '))
for idx, j in enumerate(line):
if idx == 0:
continue
self.bonds += ['Excl', int(line[0]), int(j)]
std_head = {'atoms': 'i type resnr residue atom cgnr charge mass'.split(' '),
'bonds': 'i j funct length force'.split(' '),
'angles': 'i j k funct angle force'.split(' '),
'dihedrals': 'i j k l funct angle force'.split(' '),
'constraints': 'i j funct length'.split(' '),
'virtual_sites': 'i j k l func a b c'.split(' ')}
# the line following a new section is usually a header. This can be checked by with ;
if headers and (line[0] != ';' or current_section in std_head):
headers = False
# standard headers
if current_section in std_head:
current_headers = std_head[current_section]
c_dict = OrderedDict()
for head in current_headers:
c_dict[head] = []
else:
current_headers = None
# unknown section; skip it
if headers is False and current_headers is None and current_headers != 'exclusions':
warnings.warn('Found unknown section ' + current_section + ' in ITP file, skipping')
continue
if headers and current_section != ' ':
current_headers = re.split('[ \t]+', line.strip(';').lstrip().rstrip())
c_dict = OrderedDict()
for head in current_headers:
c_dict[head] = []
headers = False
continue
if ';' in line: # comment without section head
continue
listline = re.split('[ \t]+', line.lstrip().rstrip())
if len(listline) > len(current_headers):
warnings.warn('Found bad data/headers in section ' + str(current_section) + ' skipping')
headers = False
current_headers = None
continue
for head_idx in range(len(listline)):
c_dict[current_headers[head_idx]].append(listline[head_idx])
self.gro_section[current_section] = c_dict
except Exception as e:
itp_file.close()
raise e
itp_file.close()
# reparse some of the stuff that should have other types
for key in self.INSANE_field:
if key == 'charge':
self.INSANE_field[key] = float(self.INSANE_field[key])
continue
if 'name' not in key:
self.INSANE_field[key] = re.split('[ \t]+', self.INSANE_field[key])
for key in self.gro_section:
for subkey in self.gro_section[key]:
if subkey == 'i' or subkey == 'j' or subkey == 'k' or subkey == 'l' or subkey == 'id' or subkey == 'funct':
self.gro_section[key][subkey] = list(map(int, self.gro_section[key][subkey]))
if subkey == 'angle' or subkey == 'length' or 'force' in subkey or 'charge' in subkey:
self.gro_section[key][subkey] = list(map(float, self.gro_section[key][subkey]))
# append the beads and other topology information
idname = 'id' if 'id' in self.gro_section['atoms'] else list(self.gro_section['atoms'].keys())[0]
atomname = 'atom' if 'atom' in self.gro_section['atoms'] else list(self.gro_section['atoms'].keys())[4]
bforcename = 'force.c.' if 'force.c.' in self.gro_section['bonds'] else list(self.gro_section['bonds'].keys())[4]
aforcename = 'force.c.' if 'force.c.' in self.gro_section['angles'] else list(self.gro_section['angles'].keys())[5]
if 'dihedrals' in self.gro_section:
dforcename = 'force.c' if 'force.c' in self.gro_section['dihedrals'] else list(self.gro_section['dihedrals'].keys())[6]
for beadidx in range(len(self.gro_section['atoms'][atomname])):
self.p_num.append(self.gro_section['atoms'][idname][beadidx])
bt = MartiniMappings.type_from_default(self.gro_section['atoms']['type'][beadidx])
bch = self.gro_section['atoms']['charge'][beadidx]
bpos = np.array([0., 0., self.p_num[-1] * 4.7])
self.beads += CGBead(position=bpos, charge=bch, beadtype=bt, mass=72.0)
self.internal_gro_names.append(self.gro_section['atoms']['atom'][beadidx])
if 'bonds' in self.gro_section:
for bondidx in range(len(self.gro_section['bonds']['i'])):
self.b_types += ['MartiniBond-' + str(MartiniChain._MARTINI_BOND_INDEX),
self.gro_section['bonds'][bforcename][bondidx],
self.gro_section['bonds']['length'][bondidx]]
self.bonds += ['MartiniBond-' + str(MartiniChain._MARTINI_BOND_INDEX),
self.gro_section['bonds']['i'][bondidx],
self.gro_section['bonds']['j'][bondidx]]
MartiniChain._MARTINI_BOND_INDEX += 1
if 'angles' in self.gro_section:
for angleidx in range(len(self.gro_section['angles']['j'])):
self.a_types += ['MartiniAngle-' + str(MartiniChain._MARTINI_ANGLE_INDEX),
self.gro_section['angles'][aforcename][angleidx],
math.radians(self.gro_section['angles']['angle'][angleidx])]
self.angles += ['MartiniAngle-' + str(MartiniChain._MARTINI_ANGLE_INDEX),
self.gro_section['angles']['i'][angleidx],
self.gro_section['angles']['j'][angleidx],
self.gro_section['angles']['k'][angleidx]]
MartiniChain._MARTINI_ANGLE_INDEX += 1
if 'dihedrals' in self.gro_section:
for dieidx in range(len(self.gro_section['dihedrals']['k'])):
self.d_types += ['MartiniDihedral-' + str(MartiniChain._MARTINI_DIHEDRAL_INDEX),
self.gro_section['dihedrals'][dforcename][dieidx],
math.radians(self.gro_section['dihedrals']['angle'][dieidx])]
self.dihedrals += ['MartiniDihedral-' + str(MartiniChain._MARTINI_DIHEDRAL_INDEX),
self.gro_section['dihedrals']['i'][dieidx],
self.gro_section['dihedrals']['j'][dieidx],
self.gro_section['dihedrals']['k'][dieidx],
self.gro_section['dihedrals']['l'][dieidx]]
MartiniChain._MARTINI_DIHEDRAL_INDEX += 1
if 'constraints' in self.gro_section:
for constidx in range(len(self.gro_section['constraints']['i'])):
self.constraints += Composite.Constraint((self.gro_section['constraints']['i'][constidx],
self.gro_section['constraints']['j'][constidx]),
self.gro_section['constraints']['length'][constidx])
# parse coordinates if supplied
if gro_path is not None:
gro_file = None
try:
gro_file = open(gro_path)
all_gro_lines = gro_file.readlines()
numatoms = int(all_gro_lines[1])
if numatoms != len(self.beads):
raise ParsingError('Could not match .gro file to supplied .itp topology: mismatching atom numbers')
for line in all_gro_lines[2:2 + len(self.beads)]:
subline = re.split('[ \t]+', line.rstrip().lstrip())
try:
bidx = self.internal_gro_names.index(subline[1])
bidx2 = self.p_num.index(int(subline[2])) # both of these should match
if bidx != bidx2:
raise ParsingError('Linear index doesn''t match the atom name indexes')
pos = np.array([float(subline[3]), float(subline[4]), float(subline[5])])
self.beads[bidx].position = pos
except ValueError:
raise ParsingError('Could not match atom name in .gro file to names of .itp file')
except ParsingError:
warnings.warn('Could not parse full topology; .gro parser returned: ' + str(ParsingError))
finally:
if gro_file is not None:
gro_file.close()
# TODO: better parsing here, but all .gro files have put the head at the top while the layers head at bottom
self.rotate([0., 1., 0., 0.])
self.zero_pos = [0., 0., 0.]
elif NETWORKX_FOUND: # attempt a construction by placing things sequentialy in the bonds
self.pnum_offset = 0 # match indexes with nodes
built_particles = [0]
# goal is to iteratively add bonded particles, starting from index 0 since we don't know any better. The
# constraints take precedence over any bond and must always be obeyed
while len(built_particles) < len(self.p_num):
# check for constraints
constraint_to_build = None
for built in built_particles:
if built in self.constraints.graph and len(self.constraints.graph[built]) != 0:
constraint_to_build = None
for item in self.constraints.graph[built]:
if item not in built_particles:
constraint_to_build = item
break
if constraint_to_build is not None:
break
if constraint_to_build is not None:
# we have a constraint to build, find first 3 other constraints leading to this one
# (three constraints leads to a single point if they are not degenerate)
cons_list = dict(copy(self.constraints.graph[constraint_to_build]))
for key in cons_list.keys():
if key not in built_particles:
del cons_list[key]
numel_cons = len(cons_list)
if numel_cons == 1:
dist = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[0])).distance
rv = np.random.uniform(0., 1., 3)
rv *= dist / np.linalg.norm(rv)
self.positions[constraint_to_build, :] = rv + self.positions[list(cons_list.keys())[0], :]
if numel_cons == 2:
dist1 = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[0])).distance
dist2 = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[1])).distance
p1 = self.positions[list(cons_list.keys())[0], :]
p2 = self.positions[list(cons_list.keys())[1], :]
self.positions[constraint_to_build, :] = Util.bilaterate(p1, p2, dist1, dist2)
if numel_cons == 3:
dist1 = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[0])).distance
dist2 = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[1])).distance
dist3 = self.constraints.get_topology((constraint_to_build, list(cons_list.keys())[2])).distance
p1 = self.positions[list(cons_list.keys())[0], :]
p2 = self.positions[list(cons_list.keys())[1], :]
p3 = self.positions[list(cons_list.keys())[2], :]
self.positions[constraint_to_build, :] = Util.trilaterate(p1, p2, p3, dist1, dist2, dist3)
built_particles.append(constraint_to_build)
continue
# same procedure, but for bonds instead
bond_to_build = None
for built in built_particles:
if built in self.bonds.graph and len(self.bonds.graph[built]) != 0:
for item in self.bonds.graph[built]:
if item not in built_particles:
bond_to_build = item
break
if bond_to_build is not None:
break
if bond_to_build is None:
raise ParsingError('Unable to build the positions from topology!')
bond_list = dict(copy(self.bonds.graph[bond_to_build]))
for key in bond_list.keys():
if key not in built_particles:
del bond_list[key]
numel_bond = len(bond_list)
if numel_bond == 1:
topo_tuple = (bond_to_build, list(bond_list.keys())[0])
dist = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
rv = np.random.uniform(0., 1., 3)
rv *= dist / np.linalg.norm(rv)
self.positions[bond_to_build, :] = rv + self.positions[list(bond_list.keys())[0], :]
if numel_bond == 2:
topo_tuple = (bond_to_build, list(bond_list.keys())[0])
dist1 = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
topo_tuple = (bond_to_build, list(bond_list.keys())[1])
dist2 = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
p1 = self.positions[list(bond_list.keys())[0], :]
p2 = self.positions[list(bond_list.keys())[1], :]
try:
self.positions[bond_to_build, :] = Util.bilaterate(p1, p2, dist1, dist2)
except ValueError:
rv = np.random.uniform(0., 1., 3)
rv *= dist1 / np.linalg.norm(rv)
self.positions[bond_to_build, :] = rv + self.positions[list(bond_list.keys())[0], :]
if numel_bond >= 3:
topo_tuple = (bond_to_build, list(bond_list.keys())[0])
dist1 = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
topo_tuple = (bond_to_build, list(bond_list.keys())[1])
dist2 = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
topo_tuple = (bond_to_build, list(bond_list.keys())[2])
dist3 = self.bonds.get_topology(topo_tuple, self.bonds.graph.edges[topo_tuple[0], topo_tuple[1]])['distance_constant']
p1 = self.positions[list(bond_list.keys())[0], :]
p2 = self.positions[list(bond_list.keys())[1], :]
p3 = self.positions[list(bond_list.keys())[2], :]
try:
self.positions[bond_to_build, :] = Util.trilaterate(p1, p2, p3, dist1, dist2, dist3)
except ValueError:
try:
self.positions[bond_to_build, :] = Util.bilaterate(p1, p2, dist1, dist2)
except ValueError:
rv = np.random.uniform(0., 1., 3)
rv *= dist1 / np.linalg.norm(rv)
self.positions[bond_to_build, :] = rv + self.positions[list(bond_list.keys())[0], :]
built_particles.append(bond_to_build)