From ee67c5d92fd367e2f06e91dc244f15652a1ff7fd Mon Sep 17 00:00:00 2001 From: jonathan Date: Wed, 19 Dec 2018 19:24:46 +0000 Subject: [PATCH 01/13] first commit. The linking of residues is not fully implemented. --- test/data/.sugar.xtc_offsets.npz | Bin 0 -> 8948 bytes test/data/.water.xtc_offsets.npz | Bin 0 -> 1028 bytes test/data/global.bnd | 8 ++++++++ test/data/global.gro | 9 +++++++++ test/test_bondset.py | 24 ++++++++++++++++++++++++ 5 files changed, 41 insertions(+) create mode 100644 test/data/.sugar.xtc_offsets.npz create mode 100644 test/data/.water.xtc_offsets.npz create mode 100644 test/data/global.bnd create mode 100644 test/data/global.gro diff --git a/test/data/.sugar.xtc_offsets.npz b/test/data/.sugar.xtc_offsets.npz new file mode 100644 index 0000000000000000000000000000000000000000..7e45d60db3eddd787e9be79c268d0ac3c09fb7b9 GIT binary patch literal 8948 zcmb{20gz*Z9mnz9U9bp(AXwxG7X&$iAXqHIt&?iBOVi3#?&R9I+uo%FtAeGWX>BVd zG-Z)$=(f_@wrQKT+cq=^LPOIwZPRvV=Xt=|w{y|De1Rnd*B=eKW+tQuK8-;U)j*K+6n zPTS?S|8?0)Ga54b;9A6oiVQ8L=j8Yj+F`_LIlLJ?E0dekwk@zI{e5&>Cbp!lnxr6o z3+HKmY08+1Hf@7B>D`vD%gA=L!HgN{+@3DW;0|8cFP-#QD+oE*2kkp}xxEXly$bXz9&p{;f- zNZ-D6Q^xkAP0KJRy$-r2Bm2{a127|92htT8a?*2h{2YCyXJzsX z+O`sl(tjr1mI*&?eGe9-?<~3{qh--$xtYk6Gy&rYkaZ4m~Hw&!rt7 zz-bu>(ruaeAZ`5+7UbZ0bVH7W=)P4rA>ALQt1^5(-IM7s?O2V|a(E5hk%H9R@lJU>brW-IPz2kIUMn6j% zKZjZAzLBoVP@0~Tq^mM~8{L!X z+iB+nmSyk`x+_z6()PQsBm+6RBgei(Tkpn#^nICb$ylB?-Ge#ly_c@Z$X95?Bxa=R zt8_&Mzedl>EVdY1;WCEX&}J>8?!ugtnEjDE$x6ZJBtGw*C|g z(pRBdGX4;4ei-v|;1RkZN2+w+qc|ZwkI^+5d7L)XFe6=0&=nbalAe>}KcgLWoR-5+ z(X%r7G;RAi7Nx&Iw`JlNwDp%*kiK8hEg5gp=3ir8dVfRLW%ReSv4vUb`W;=7q2JS8 znfe25o53kL^bFmS@n>oCbC{FfHeHv|KhlOjVMaRtOqXS_Lw9BBdD{L0mSo^Xx+4>_ zwB;qt%Ym2ahK#;K8@rg5?pNum48KPAWcn|(a}LWg_*c3sQ-7oFf5(yx^yrQp`v=|s zPn?oNuhTPf)UfD{4VFbXDF+wRO&ME4n~a!~-p%N`jBZXFx4^7)_t8}u-jeRg^j5TU zYb?v4iSEkeHneS9EK2`&bX$&^Y0LJQmjgS{4LPzSZM0xUx^|)~GPpC{m8o55yA?|^ zuq)k>iQQ=H?pTn%e!3-Nd(fsmF(A zy>7ZDBZt$5BQPUfN75A;^3Zc~{3zOSG)~LmW9V6#^wPFtu_*n=(QTPHp0*BPLHbUh zTQYtkZ9WO}a$t~d$dQw2<0+Vx?o;Wi4EyMwOrJ(OPsg$huAsXzHALIbz>*BCq&qTk zCT;QKq#S$?-ITGjXw!Q!Cp`hWD#K^fJ(+$V?R-C$;6rp>M$e;-A-IDROwD}{L zlb(;#H5rN0Jvn|MZC{5)>0eLJ$oL3t{ut(^cLQCQk&Sduren18<5-r#i|DRQT}<2J zSd{)t=(bE;N?R_&yc`&%8!~!1ZTJLcr1J{8EQ1MpR*qdsTd%@`9K4!t%Gelfx(0L7 zb1hwykx$ZwBxa=RI=Ui5*VA2@`V?(XVMzu)O?PDCGqm*vEJ)ut-IDRo(&o=$UJl$y zH{?i~?)yAWNY71lO-62}4H?Wx*DZ8e25+UiGW7-8p2d<3e398I1JtxP%LOUjLS`L4eo|Vb3 z(YE`rDE$SxEfZg-t>3_c9K4@y%GeZb`X)|D&$sBBjC`9m6frHG-=WJg_+7dyQ{SWQ zB`nFn_vwxt`vKkmL!6RB)AWoS{Sj^XF;2?CpU_PiE7PV2a6)<>q-!$rQ`%6$jC4Ij zS7hj6dQOf%LOZHBEr%baXJzs++V(gWTwzExILRzoSjR$DH*3fv(Hw z3~hJ@Gt%`eU6G;Z=&nq*Y1<#MDE)t;+cNQI+S{sk!!OZ2 zIsP*3cm+$+-=$|{{8ifY8s?<;FLYf-=V;?!F)Q7FqpLFXcY02a_h|b+uqcQANw;MD zb=quL{Ko%bc{#9%Zpe|vv~dY$rOQZHWN}G_wRyJ(r2YxGQKNq z-VO88yE|Q%(SF*v2WF*vPr4#Qd(mB)veCAsSd{*~>9$PlLtE`wkiLECmW=I3o0eft zdL49KM)#+U2VhpZ52ULy?4)}#eGu(D7|Swv2;G$_7j0jTB^fxB?#Qvj=zcd&$)UsP z85uuFRg=wEtZ_=(jB_Tlvpq>UU0iuKTIJH>cfle%kY6YcKo{`%Wr?$-0;t~+m4Z2g7r{%wl?!@%3mkF&7&S1ej= zSz;j!|Na7qb{GmZvRxap1b)`2hW_10~*=~0_r(6^4d-S)Q2yXj$@p4hhz cTe2TA^-pY literal 0 HcmV?d00001 diff --git a/test/data/.water.xtc_offsets.npz b/test/data/.water.xtc_offsets.npz new file mode 100644 index 0000000000000000000000000000000000000000..514b7280ac8de1ad4be73da58d8fdc5e749d9964 GIT binary patch literal 1028 zcmWIWW@Zs#fB;2?tWUF*ZvZ(U%nii(X=%l&CB=Gq1(lQiLVW`w85znLs?}3ci<67g ztrXO4GA-0~6x7r5i%Nzy)z4 z&}keH)&fBY9UujjQHIc0v>|kh350gAh0s4-A@m$S2<;LMp#x9?;6Jm`mJXl<2!kvt z&a6tM*!L8A9T;d_PvWt;Hp=8G2h6pcKn2MqnYnnJn?_t@QNg+Xn7#V3~J=Te@!sNmq~OtA1}WD;S~IZRi1u&=v>u7g8_>c(byBBw2v)K9J^T2Jrye>cGkX literal 0 HcmV?d00001 diff --git a/test/data/global.bnd b/test/data/global.bnd new file mode 100644 index 0000000..1ccadef --- /dev/null +++ b/test/data/global.bnd @@ -0,0 +1,8 @@ +; comments begin with a semicolon +[ global ] +O1_1_GLX O1_2_GLY +O1_1_GLX O1_2_GLY C1_2_GLY + + + + diff --git a/test/data/global.gro b/test/data/global.gro new file mode 100644 index 0000000..24324b1 --- /dev/null +++ b/test/data/global.gro @@ -0,0 +1,9 @@ +Alpha-allose in water + 6 + 1GLX C1 1 2.712 2.740 1.280 + 1GLX O1 2 2.653 2.609 1.281 + 1GLX HO1 3 2.563 2.621 1.324 + 2GLY C1 4 2.744 2.804 1.414 + 2GLY O1 5 2.843 2.731 1.490 + 2GLY HO1 6 2.839 2.770 1.582 + 6.00000 6.00000 6.00000 diff --git a/test/test_bondset.py b/test/test_bondset.py index c2e53cd..cdab447 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -276,3 +276,27 @@ def test_get_lines_for_bond_dump_sample(self): self.assertIn(line, expected) self.assertNotIn(line, seen) seen.add(line) + + def test_global_bond_create(self): + DummyOptions.generate_angles=False + measure = BondSet("test/data/global.bnd", DummyOptions) + self.assertEqual(len(measure.global_connections), 2) + names = ["O1", "O1"] + resids = [1, 2] + resnames = ["GLX", "GLY"] + self.assertListEqual(resnames, measure.global_connections[0].resnames) + self.assertListEqual(resids, measure.global_connections[0].resids) + self.assertListEqual(names, measure.global_connections[0].atoms) + + def test_global_bond_get_atoms(self): + DummyOptions.generate_angles = False + measure = BondSet("test/data/global.bnd", DummyOptions) + frame = Frame("test/data/global.gro") + target_bond = [frame.residues[0]["O1"], frame.residues[1]["O1"]] + target_angle = [frame.residues[0]["O1"], frame.residues[1]["O1"], frame.residues[1]["C1"]] + atoms = measure.global_connections[0].get_atoms(frame) + self.assertListEqual(target_bond, atoms) + atoms = measure.global_connections[1].get_atoms(frame) + self.assertListEqual(target_angle, atoms) + + From 655c6d72ce40e56dd6848125d66444c00ea3cf52 Mon Sep 17 00:00:00 2001 From: jonathan Date: Wed, 19 Dec 2018 19:39:44 +0000 Subject: [PATCH 02/13] first commit. The linking of residues is not fully implemented. --- pycgtool/bondset.py | 122 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 101 insertions(+), 21 deletions(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index f78a39d..fcf6701 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -100,6 +100,47 @@ def __repr__(self): except (AttributeError, TypeError): return "".format(", ".join(self.atoms)) +class GlobalBond(Bond): + __slots__ = ["atoms", "atom_numbers", "values", "eqm", "fconst", "gromacs_type_id", "_func_form", + "resids", "resnames", "labels"] + def __init__(self, atoms, atom_numbers=None, func_form=None): + """ + Class for bonds between not necessarily adjacent residues + :param List[str] atoms: List of atom names defining the bond + :param List[int] atom_numbers: List of atom numbers defining the bond + :param func_form: Functional form to use for Boltzmann Inversion + """ + self.labels = atoms + atom_names = [] + resids = [] + resnames = [] + for atom in atoms: + try: + name, resid, resname = atom.split("_") + resid = int(resid) + atom_names.append(name) + resids.append(resid) + resnames.append(resname) + + except ValueError: + print("incorrect syntax for entry '{}' in [global] section".format(atom)) + raise SyntaxError + Bond.__init__(self, atom_names, atom_numbers=atom_numbers, func_form=func_form) + self.resids = resids + self.resnames = resnames + + def get_atoms(self, residues): + atoms = [] + for resid, resname, name in zip(self.resids, self.resnames, self.atoms): + for residue in residues: + if residue.num == resid: + if residue.name == resname: + try: + atom = residue[name] + atoms.append(atom) + except KeyError: + pass + return atoms class BondSet: """ @@ -115,6 +156,7 @@ def __init__(self, filename, options): :return: Instance of BondSet """ self._molecules = {} + self.global_connections = [] self._fconst_constr_threshold = options.constr_threshold @@ -154,8 +196,12 @@ def __init__(self, filename, options): with CFG(filename) as cfg: for mol_name, mol_section in cfg.items(): - self._molecules[mol_name] = [] - mol_bonds = self._molecules[mol_name] + is_global = False + if mol_name == "global": + is_global = True + else: + self._molecules[mol_name] = [] + mol_bonds = self._molecules[mol_name] angles_defined = False for atomlist in mol_section: @@ -180,31 +226,50 @@ def __init__(self, filename, options): if {x for x in atomlist if atomlist.count(x) > 1}: raise ValueError("Defined bond '{0}' contains duplicate atoms".format(atomlist)) - - mol_bonds.append(Bond(atoms=atomlist, func_form=func_form)) + if is_global: + self.global_connections.append(GlobalBond(atoms=atomlist, func_form=func_form)) + else: + mol_bonds.append(Bond(atoms=atomlist, func_form=func_form)) if len(atomlist) > 2: angles_defined = True if not angles_defined: - angles, dihedrals = self._create_angles(mol_bonds) - - if options.generate_angles: - for atomlist in angles: - mol_bonds.append( - Bond( - atoms=atomlist, - func_form=self._functional_forms[3](circular_mean, circular_variance) + if is_global: + if options.generate_angles or options.generate_dihedrals: + logger.warning("Automated generation of angles or dihedrals between " + "residues not implemented! Please specifiy angles and dihedrals in [global] " + "section of *.bnd file") + else: + + angles, dihedrals = self._create_angles(mol_bonds) + + if options.generate_angles: + for atomlist in angles: + mol_bonds.append( + Bond( + atoms=atomlist, + func_form=self._functional_forms[3](circular_mean, circular_variance) + ) ) - ) - - if options.generate_dihedrals: - for atomlist in dihedrals: - mol_bonds.append( - Bond( - atoms=atomlist, - func_form=self._functional_forms[4](circular_mean, circular_variance) + + if options.generate_dihedrals: + for atomlist in dihedrals: + mol_bonds.append( + Bond( + atoms=atomlist, + func_form=self._functional_forms[4](circular_mean, circular_variance) + ) ) - ) + if is_global: + self._molecules["global"] = self.global_connections + + def join_residues(self, residues): + """ + join together residues connected by global bonds into new molecules + :param residues: Iterator of Residue objects + :return: + """ + pass @staticmethod def _create_angles(mol_bonds): @@ -455,6 +520,16 @@ def calc_dihedral(atoms): except ZeroDivisionError as e: e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.atoms)),) raise e + #calculate values for global bonds# + if len(self.global_connections) > 0: + for bond in self.global_connections: + try: + atoms = bond.get_atoms(frame) + val = calc[len(atoms)](atoms) + bond.values.append(val) + except ZeroDivisionError as e: + e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.labels)),) + raise e def boltzmann_invert(self, progress=False): """ @@ -463,6 +538,8 @@ def boltzmann_invert(self, progress=False): :param progress: Display a progress bar using tqdm if available """ bond_iter = itertools.chain(*self._molecules.values()) + if len(self.global_connections) > 0: + bond_iter = itertools.chain(bond_iter, self.global_connections) if progress: total = sum(map(len, self._molecules.values())) bond_iter = tqdm(bond_iter, total=total, ncols=80) @@ -514,6 +591,9 @@ def dump_values(self, target_number=10000): if bonds: lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) file_write_lines("{0}_dihedral.dat".format(mol), lines) + if len(self.global_connections) > 0: + global_bonds = self.global_connections + lines = BondSet._get_lines_for_bond_dump(glo) def __len__(self): return len(self._molecules) From 26ce80d879277a1577811929539c47c446522254 Mon Sep 17 00:00:00 2001 From: Jon Shearer Date: Sat, 22 Dec 2018 18:58:02 +0000 Subject: [PATCH 03/13] Now can connect residues using the [ global ] section. Will only work for one molecule in water at the moment --- README.md | 1 + pycgtool/bondset.py | 267 ++++++++++++++++++++++++------- pycgtool/frame.py | 16 ++ pycgtool/pycgtool.py | 10 +- pycgtool/util.py | 37 +++++ test/data/.sugar.xtc_offsets.npz | Bin 8948 -> 0 bytes test/data/.water.xtc_offsets.npz | Bin 1028 -> 0 bytes test/data/global-cg.gro | 7 + test/data/global.bnd | 4 +- test/data/global.map | 8 + test/test_bondset.py | 37 ++++- test/test_util.py | 7 +- 12 files changed, 322 insertions(+), 72 deletions(-) delete mode 100644 test/data/.sugar.xtc_offsets.npz delete mode 100644 test/data/.water.xtc_offsets.npz create mode 100644 test/data/global-cg.gro create mode 100644 test/data/global.map diff --git a/README.md b/README.md index ff69d2f..f800150 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,7 @@ # PyCGTOOL +WARNING THIS IS A VERY EXPERIMENTAL BRANCH - EXERCISE CAUTION! Please see http://pycgtool.readthedocs.io/en/master/ for full documentation. diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index fcf6701..01758d7 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -8,6 +8,7 @@ import math import logging import numpy as np +from copy import deepcopy try: @@ -17,6 +18,7 @@ from .mapping import VirtualMap from .functionalforms import FunctionalForms +from .frame import Molecule from .parsers.cfg import CFG from .util import ( circular_mean, @@ -29,7 +31,8 @@ vector_angle, vector_angle_signed, vector_cross, - vector_len + vector_len, + merge_list_of_lists ) logger = logging.getLogger(__name__) @@ -142,6 +145,31 @@ def get_atoms(self, residues): pass return atoms + def get_residue_ids(self, frame): + residue_ids = [] + for resid, resname, name in zip(self.resids, self.resnames, self.atoms): + for ind, residue in enumerate(frame): + if residue.num == resid: + if residue.name == resname: + residue_ids.append(ind) + return residue_ids + + def populate_ids(self, mol_beads): + ids = [] + for resid, resname, name in zip(self.resids, self.resnames, self.atoms): + for ind, beads in enumerate(mol_beads, start=1): + index = [bead.name for bead in beads] + if ind == resid: + try: + ids.extend([ beads[index.index(atom)].num for atom in index if atom == name]) + except ValueError as e: + missing = [atom for atom in self.atoms if atom.lstrip("+-") not in index] + e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, resname),) + raise + except KeyError: + pass + self.atom_numbers = ids + class BondSet: """ Class used to perform bond measurements in a Frame. @@ -260,16 +288,6 @@ def __init__(self, filename, options): func_form=self._functional_forms[4](circular_mean, circular_variance) ) ) - if is_global: - self._molecules["global"] = self.global_connections - - def join_residues(self, residues): - """ - join together residues connected by global bonds into new molecules - :param residues: Iterator of Residue objects - :return: - """ - pass @staticmethod def _create_angles(mol_bonds): @@ -293,9 +311,14 @@ def get_bonds(self, mol, natoms, select=lambda x: True): :param function select: Optional lambda, return only bonds for which this is True :return List[Bond]: List of bonds """ + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + else: + bonds = self._molecules[mol] + if natoms == -1: - return [bond for bond in self._molecules[mol] if select(bond)] - return [bond for bond in self._molecules[mol] if len(bond.atoms) == natoms and select(bond)] + return [bond for bond in bonds if select(bond)] + return [bond for bond in bonds if len(bond.atoms) == natoms and select(bond)] def get_bond_lengths(self, mol, with_constr=False): """ @@ -305,10 +328,15 @@ def get_bond_lengths(self, mol, with_constr=False): :param with_constr: Include constraints? :return: List of bonds """ + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + else: + bonds = self._molecules[mol] + if with_constr: - return [bond for bond in self._molecules[mol] if len(bond.atoms) == 2] + return [bond for bond in bonds if len(bond.atoms) == 2] else: - return [bond for bond in self._molecules[mol] if len(bond.atoms) == 2 and bond.fconst < self._fconst_constr_threshold] + return [bond for bond in bonds if len(bond.atoms) == 2 and bond.fconst < self._fconst_constr_threshold] def get_bond_length_constraints(self, mol): """ @@ -317,7 +345,11 @@ def get_bond_length_constraints(self, mol): :param mol: Molecule name to return bonds for :return: List of bonds """ - return [bond for bond in self._molecules[mol] if len(bond.atoms) == 2 and bond.fconst >= self._fconst_constr_threshold] + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + else: + bonds = self._molecules[mol] + return [bond for bond in bonds if len(bond.atoms) == 2 and bond.fconst >= self._fconst_constr_threshold] def get_bond_angles(self, mol, exclude_triangle=True): """ @@ -327,7 +359,11 @@ def get_bond_angles(self, mol, exclude_triangle=True): :param exclude_triangle: Exclude angles that are part of a triangle? :return: List of bonds """ - angles = [bond for bond in self._molecules[mol] if len(bond.atoms) == 3] + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + else: + bonds = self._molecules[mol] + angles = [bond for bond in bonds if len(bond.atoms) == 3] if exclude_triangle: edges = [tuple(bond.atoms) for bond in self.get_bond_lengths(mol, with_constr=True)] @@ -350,7 +386,11 @@ def get_bond_dihedrals(self, mol): :param mol: Molecule name to return bonds for :return: List of bonds """ - return [bond for bond in self._molecules[mol] if len(bond.atoms) == 4] + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + else: + bonds = self._molecules[mol] + return [bond for bond in bonds if len(bond.atoms) == 4] def get_virtual_beads(self, mapping): """ @@ -385,6 +425,61 @@ def _populate_atom_numbers(self, mapping): e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, mol),) raise + def connect_residues(self, frame, mapping): + """ + connects residues together to form new molecules + :param frame: Frame from which to determine inter residue connections + """ + if len(self.global_connections) > 0: + fragments = [] + for bond in self.global_connections: + fragments.append(bond.get_residue_ids(frame)) + + self._populate_atom_numbers(mapping) + + molecule_internal_ids = merge_list_of_lists(fragments) + if len(molecule_internal_ids) > 1: + print("More than one multi residue molecule detected - this is not yet supported by pycgtool") + raise NotImplementedError + molecule_internal_ids = molecule_internal_ids[0] + resnames = [frame[mol_id].name for mol_id in molecule_internal_ids] + all_bonds = [] + all_beads = [] + start = 0 + for resid, resname in enumerate(resnames): + if resname in mapping: + beads = list(map(deepcopy, mapping[resname])) + if len(self._molecules) > 0: + bonds = list(map(deepcopy, self._molecules[resname])) + all_bonds.append(bonds) + for i, bead in enumerate(beads): + bead.num = start + i + + if len(self._molecules) > 0: + index = [bead.name for bead in beads] + for bond in bonds: + try: + bond.atom_numbers = [index.index(atom.lstrip("+-")) for atom in bond.atoms] + except ValueError as e: + missing = [atom for atom in bond.atoms if atom.lstrip("+-") not in index] + e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, resname),) + raise + bond.atom_numbers = [start + ind for ind in bond.atom_numbers] + + all_beads.append(beads) + start = beads[-1].num + 1 + + + for bond in self.global_connections: + bond.populate_ids(all_beads) + + all_bonds.extend(self.global_connections) + molecule = Molecule(resnames, all_bonds, all_beads) + molecule_name = "_".join(resnames) + self._molecules[molecule_name] = molecule + else: + return + def write_itp(self, filename, mapping): """ Output a GROMACS .itp file containing atoms/beads and bonded terms. @@ -427,47 +522,96 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= # Print molecule not_calc = " Parameters have not been calculated." for mol in self._molecules: - if mol not in mapping: - logger.warning("Molecule '{0}' present in bonding file, but not in mapping.".format(mol) + not_calc) - continue - if not all((bond.fconst is not None for bond in self._molecules[mol])): - logger.warning("Molecule '{0}' has no measured bond values.".format(mol) + not_calc) - continue - - ret_lines.append("\n[ moleculetype ]") - ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) - - ret_lines.append("\n[ atoms ]") - - for i, bead in enumerate(mapping[mol], start=1): - # atnum type resnum resname atname c-group charge (mass) - if isinstance(bead, VirtualMap): - ret_lines.append(atom_template["mass"].format( - i, bead.type, 1, mol, bead.name, i, bead.charge, bead.mass - )) - else: - ret_lines.append(atom_template["nomass"].format( - i, bead.type, 1, mol, bead.name, i, bead.charge - )) - - virtual_beads = self.get_virtual_beads(mapping[mol]) - if len(virtual_beads) != 0: - ret_lines.append("\n[ virtual_sitesn ]") - excl_lines = ["\n[ exclusions ]"] #exlusions section for virtual sites - for vbead in virtual_beads: - CGids = [bead.num + 1 for bead in mapping[mol] if bead.name in vbead.atoms] - CGids.sort() - CGids_string = " ".join(map(str, CGids)) - ret_lines.append("{0:^6d} {1:^6d} {2}".format(vbead.num+1, vbead.gromacs_type_id, CGids_string)) - vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string - excl_lines.append(vsite_exclusions) - ret_lines.extend(excl_lines) - - ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) - ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) - + if not isinstance(self._molecules[mol], Molecule): + if mol not in mapping: + logger.warning("Molecule '{0}' present in bonding file, but not in mapping.".format(mol) + not_calc) + continue + if not all((bond.fconst is not None for bond in self._molecules[mol])): + logger.warning("Molecule '{0}' has no measured bond values.".format(mol) + not_calc) + continue + + ret_lines.append("\n[ moleculetype ]") + ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) + + ret_lines.append("\n[ atoms ]") + for i, bead in enumerate(mapping[mol], start=1): + # atnum type resnum resname atname c-group charge (mass) + if isinstance(bead, VirtualMap): + ret_lines.append(atom_template["mass"].format( + i, bead.type, 1, mol, bead.name, i, bead.charge, bead.mass + )) + else: + ret_lines.append(atom_template["nomass"].format( + i, bead.type, 1, mol, bead.name, i, bead.charge + )) + + virtual_beads = self.get_virtual_beads(mapping[mol]) + if len(virtual_beads) != 0: + ret_lines.append("\n[ virtual_sitesn ]") + excl_lines = ["\n[ exclusions ]"] #exlusions section for virtual sites + for vbead in virtual_beads: + CGids = [bead.num + 1 for bead in mapping[mol] if bead.name in vbead.atoms] + CGids.sort() + CGids_string = " ".join(map(str, CGids)) + ret_lines.append("{0:^6d} {1:^6d} {2}".format(vbead.num+1, vbead.gromacs_type_id, CGids_string)) + vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string + excl_lines.append(vsite_exclusions) + ret_lines.extend(excl_lines) + + ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) + ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) + + else: + mol_count = 0 + for res in self._molecules[mol].resnames: + if res not in mapping: + logger.warning( + "Residue '{0}' present in bonding file, but not in mapping.".format(res) + not_calc) + continue + if res not in self._molecules: + logger.warning(("Residue '{0}' has no internal bonds").format(res)) + + ret_lines.append("\n[ moleculetype ]") + mol_name = "mol_" + str(mol_count) + ret_lines.append("{0:4s} {1:4d}".format(mol_name, 1)) + + ret_lines.append("\n[ atoms ]") + + for resid, beads in self._molecules[mol].resid_to_beads.items(): + for bead in beads: + # atnum type resnum resname atname c-group charge (mass) + if isinstance(bead, VirtualMap): + ret_lines.append(atom_template["mass"].format( + bead.num + 1, bead.type, resid, mol, bead.name, bead.num + 1, bead.charge, bead.mass + )) + else: + ret_lines.append(atom_template["nomass"].format( + bead.num + 1, bead.type, resid, mol, bead.name, bead.num + 1, bead.charge + )) + + virtual_beads = self.get_virtual_beads(self._molecules[mol].beads) + if len(virtual_beads) != 0: + ret_lines.append("\n[ virtual_sitesn ]") + excl_lines = ["\n[ exclusions ]"] # exlusions section for virtual sites + for vbead in virtual_beads: + CGids = [bead.num + 1 for bead in self._molecules[mol].beads if bead.name in vbead.atoms] + CGids.sort() + CGids_string = " ".join(map(str, CGids)) + ret_lines.append( + "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) + vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string + excl_lines.append(vsite_exclusions) + ret_lines.extend(excl_lines) + + ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) + ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) + ret_lines.extend( + write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) + ret_lines.extend( + write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) + mol_count += 1 return ret_lines def apply(self, frame): @@ -577,6 +721,8 @@ def dump_values(self, target_number=10000): for mol in self._molecules: if mol == "SOL": continue + if isinstance(self._molecules[mol], Molecule): + continue bonds = self.get_bond_lengths(mol, with_constr=True) if bonds: lines = BondSet._get_lines_for_bond_dump(bonds, target_number) @@ -593,7 +739,8 @@ def dump_values(self, target_number=10000): file_write_lines("{0}_dihedral.dat".format(mol), lines) if len(self.global_connections) > 0: global_bonds = self.global_connections - lines = BondSet._get_lines_for_bond_dump(glo) + lines = BondSet._get_lines_for_bond_dump(global_bonds, target_number) + file_write_lines("global_length.dat".format(mol), lines) def __len__(self): return len(self._molecules) diff --git a/pycgtool/frame.py b/pycgtool/frame.py index bcf2dda..ea9e2d3 100644 --- a/pycgtool/frame.py +++ b/pycgtool/frame.py @@ -66,6 +66,19 @@ def add_missing_data(self, other): setattr(self, attr, getattr(other, attr)) +class Molecule: + """ + Holds data for a molecule comprised of multiple residues + """ + __slots__ = ["resnames", "bonds", "beads", "resid_to_beads"] + + def __init__(self, resnames, bonds, beads): + self.resnames = resnames + self.resid_to_beads = dict(zip(range(len(resnames)), beads)) + self.bonds = bonds + self.beads = np.array(beads).flatten().tolist() + + class Residue: """ Hold data for a residue - list of atoms @@ -96,6 +109,9 @@ def __getitem__(self, item): def __len__(self): return len(self.atoms) + def __contains__(self, item): + return item in self.name_to_num + def add_atom(self, atom): """ Add an Atom to this Residue and store location in index diff --git a/pycgtool/pycgtool.py b/pycgtool/pycgtool.py index dc086aa..dd40703 100755 --- a/pycgtool/pycgtool.py +++ b/pycgtool/pycgtool.py @@ -63,10 +63,14 @@ def main_loop(): if args.map: logger.info("Beginning Boltzmann inversion") bonds.boltzmann_invert(progress=(not args.quiet)) + bonds.connect_residues(cgframe, mapping) if config.output_forcefield: - logger.info("Creating GROMACS forcefield directory") - ForceField(config.output_name).write(config.output_name, mapping, bonds) - logger.info("GROMACS forcefield directory created") + if len(bonds.global_connections) > 0: + logger.warning("Cannot create GROMACS forcefield when using [ global ] connections!") + else: + logger.info("Creating GROMACS forcefield directory") + ForceField(config.output_name).write(config.output_name, mapping, bonds) + logger.info("GROMACS forcefield directory created") else: bonds.write_itp(config.output_name + ".itp", mapping=mapping) diff --git a/pycgtool/util.py b/pycgtool/util.py index c4d4492..6230043 100644 --- a/pycgtool/util.py +++ b/pycgtool/util.py @@ -67,6 +67,21 @@ def vector_cross(u, v): res[2] = u[0] * v[1] - u[1] * v[0] return res +def shift_func_value(func_to_decorate, shift, periodic=None): + """ + shift value that a function returns + :param func_to_decorate: + :param shift: value to shift return value by + :return: new_func: new function + """ + + def new_func(*original_args, **original_kwargs): + value = func_to_decorate(*original_args, **original_kwargs) - shift + if periodic is not None: + value -= np.rint(value / two_pi) * two_pi + return value + return new_func + def circular_mean(values): """ return average of values on a cirle @@ -293,6 +308,28 @@ def backup_file(name): logger.warning("Existing file {0} backed up as {1}".format(name, new_name)) return new_name +def merge_list_of_lists(l): + out = [] + while len(l) > 0: + first, *rest = l + first = set(first) + + lf = -1 + while len(first) > lf: + lf = len(first) + + rest2 = [] + for r in rest: + if len(first.intersection(set(r))) > 0: + first |= set(r) + else: + rest2.append(r) + rest = rest2 + + out.append(first) + l = rest + return [list(sorted(x)) for x in out] + def sliding(vals): """ diff --git a/test/data/.sugar.xtc_offsets.npz b/test/data/.sugar.xtc_offsets.npz deleted file mode 100644 index 7e45d60db3eddd787e9be79c268d0ac3c09fb7b9..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8948 zcmb{20gz*Z9mnz9U9bp(AXwxG7X&$iAXqHIt&?iBOVi3#?&R9I+uo%FtAeGWX>BVd zG-Z)$=(f_@wrQKT+cq=^LPOIwZPRvV=Xt=|w{y|De1Rnd*B=eKW+tQuK8-;U)j*K+6n zPTS?S|8?0)Ga54b;9A6oiVQ8L=j8Yj+F`_LIlLJ?E0dekwk@zI{e5&>Cbp!lnxr6o z3+HKmY08+1Hf@7B>D`vD%gA=L!HgN{+@3DW;0|8cFP-#QD+oE*2kkp}xxEXly$bXz9&p{;f- zNZ-D6Q^xkAP0KJRy$-r2Bm2{a127|92htT8a?*2h{2YCyXJzsX z+O`sl(tjr1mI*&?eGe9-?<~3{qh--$xtYk6Gy&rYkaZ4m~Hw&!rt7 zz-bu>(ruaeAZ`5+7UbZ0bVH7W=)P4rA>ALQt1^5(-IM7s?O2V|a(E5hk%H9R@lJU>brW-IPz2kIUMn6j% zKZjZAzLBoVP@0~Tq^mM~8{L!X z+iB+nmSyk`x+_z6()PQsBm+6RBgei(Tkpn#^nICb$ylB?-Ge#ly_c@Z$X95?Bxa=R zt8_&Mzedl>EVdY1;WCEX&}J>8?!ugtnEjDE$x6ZJBtGw*C|g z(pRBdGX4;4ei-v|;1RkZN2+w+qc|ZwkI^+5d7L)XFe6=0&=nbalAe>}KcgLWoR-5+ z(X%r7G;RAi7Nx&Iw`JlNwDp%*kiK8hEg5gp=3ir8dVfRLW%ReSv4vUb`W;=7q2JS8 znfe25o53kL^bFmS@n>oCbC{FfHeHv|KhlOjVMaRtOqXS_Lw9BBdD{L0mSo^Xx+4>_ zwB;qt%Ym2ahK#;K8@rg5?pNum48KPAWcn|(a}LWg_*c3sQ-7oFf5(yx^yrQp`v=|s zPn?oNuhTPf)UfD{4VFbXDF+wRO&ME4n~a!~-p%N`jBZXFx4^7)_t8}u-jeRg^j5TU zYb?v4iSEkeHneS9EK2`&bX$&^Y0LJQmjgS{4LPzSZM0xUx^|)~GPpC{m8o55yA?|^ zuq)k>iQQ=H?pTn%e!3-Nd(fsmF(A zy>7ZDBZt$5BQPUfN75A;^3Zc~{3zOSG)~LmW9V6#^wPFtu_*n=(QTPHp0*BPLHbUh zTQYtkZ9WO}a$t~d$dQw2<0+Vx?o;Wi4EyMwOrJ(OPsg$huAsXzHALIbz>*BCq&qTk zCT;QKq#S$?-ITGjXw!Q!Cp`hWD#K^fJ(+$V?R-C$;6rp>M$e;-A-IDROwD}{L zlb(;#H5rN0Jvn|MZC{5)>0eLJ$oL3t{ut(^cLQCQk&Sduren18<5-r#i|DRQT}<2J zSd{)t=(bE;N?R_&yc`&%8!~!1ZTJLcr1J{8EQ1MpR*qdsTd%@`9K4!t%Gelfx(0L7 zb1hwykx$ZwBxa=RI=Ui5*VA2@`V?(XVMzu)O?PDCGqm*vEJ)ut-IDRo(&o=$UJl$y zH{?i~?)yAWNY71lO-62}4H?Wx*DZ8e25+UiGW7-8p2d<3e398I1JtxP%LOUjLS`L4eo|Vb3 z(YE`rDE$SxEfZg-t>3_c9K4@y%GeZb`X)|D&$sBBjC`9m6frHG-=WJg_+7dyQ{SWQ zB`nFn_vwxt`vKkmL!6RB)AWoS{Sj^XF;2?CpU_PiE7PV2a6)<>q-!$rQ`%6$jC4Ij zS7hj6dQOf%LOZHBEr%baXJzs++V(gWTwzExILRzoSjR$DH*3fv(Hw z3~hJ@Gt%`eU6G;Z=&nq*Y1<#MDE)t;+cNQI+S{sk!!OZ2 zIsP*3cm+$+-=$|{{8ifY8s?<;FLYf-=V;?!F)Q7FqpLFXcY02a_h|b+uqcQANw;MD zb=quL{Ko%bc{#9%Zpe|vv~dY$rOQZHWN}G_wRyJ(r2YxGQKNq z-VO88yE|Q%(SF*v2WF*vPr4#Qd(mB)veCAsSd{*~>9$PlLtE`wkiLECmW=I3o0eft zdL49KM)#+U2VhpZ52ULy?4)}#eGu(D7|Swv2;G$_7j0jTB^fxB?#Qvj=zcd&$)UsP z85uuFRg=wEtZ_=(jB_Tlvpq>UU0iuKTIJH>cfle%kY6YcKo{`%Wr?$-0;t~+m4Z2g7r{%wl?!@%3mkF&7&S1ej= zSz;j!|Na7qb{GmZvRxap1b)`2hW_10~*=~0_r(6^4d-S)Q2yXj$@p4hhz cTe2TA^-pY diff --git a/test/data/.water.xtc_offsets.npz b/test/data/.water.xtc_offsets.npz deleted file mode 100644 index 514b7280ac8de1ad4be73da58d8fdc5e749d9964..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1028 zcmWIWW@Zs#fB;2?tWUF*ZvZ(U%nii(X=%l&CB=Gq1(lQiLVW`w85znLs?}3ci<67g ztrXO4GA-0~6x7r5i%Nzy)z4 z&}keH)&fBY9UujjQHIc0v>|kh350gAh0s4-A@m$S2<;LMp#x9?;6Jm`mJXl<2!kvt z&a6tM*!L8A9T;d_PvWt;Hp=8G2h6pcKn2MqnYnnJn?_t@QNg+Xn7#V3~J=Te@!sNmq~OtA1}WD;S~IZRi1u&=v>u7g8_>c(byBBw2v)K9J^T2Jrye>cGkX diff --git a/test/data/global-cg.gro b/test/data/global-cg.gro new file mode 100644 index 0000000..ba30404 --- /dev/null +++ b/test/data/global-cg.gro @@ -0,0 +1,7 @@ +Alpha-allose in water + 4 + 1GLX S1 1 2.712 2.740 1.280 + 1GLX S2 2 2.653 2.609 1.281 + 2GLY S1 3 2.744 2.804 1.414 + 2GLY S2 4 2.843 2.731 1.490 + 6.00000 6.00000 6.00000 diff --git a/test/data/global.bnd b/test/data/global.bnd index 1ccadef..cab2579 100644 --- a/test/data/global.bnd +++ b/test/data/global.bnd @@ -1,7 +1,7 @@ ; comments begin with a semicolon [ global ] -O1_1_GLX O1_2_GLY -O1_1_GLX O1_2_GLY C1_2_GLY +S2_1_GLX S1_2_GLY +S2_1_GLX S1_2_GLY S2_2_GLY diff --git a/test/data/global.map b/test/data/global.map new file mode 100644 index 0000000..dd22119 --- /dev/null +++ b/test/data/global.map @@ -0,0 +1,8 @@ +; these are fake bead types +[ GLY ] +S1 TC3 C1 +S2 TC3 O1 HO1 + +[ GLX ] +S1 TC3 C1 +S2 TC3 O1 HO1 diff --git a/test/test_bondset.py b/test/test_bondset.py index cdab447..56877aa 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -281,7 +281,7 @@ def test_global_bond_create(self): DummyOptions.generate_angles=False measure = BondSet("test/data/global.bnd", DummyOptions) self.assertEqual(len(measure.global_connections), 2) - names = ["O1", "O1"] + names = ["S2", "S1"] resids = [1, 2] resnames = ["GLX", "GLY"] self.assertListEqual(resnames, measure.global_connections[0].resnames) @@ -291,12 +291,37 @@ def test_global_bond_create(self): def test_global_bond_get_atoms(self): DummyOptions.generate_angles = False measure = BondSet("test/data/global.bnd", DummyOptions) - frame = Frame("test/data/global.gro") - target_bond = [frame.residues[0]["O1"], frame.residues[1]["O1"]] - target_angle = [frame.residues[0]["O1"], frame.residues[1]["O1"], frame.residues[1]["C1"]] - atoms = measure.global_connections[0].get_atoms(frame) + cgframe = Frame("test/data/global-cg.gro") + target_bond = [cgframe.residues[0]["S2"], cgframe.residues[1]["S1"]] + target_angle = [cgframe.residues[0]["S2"], cgframe.residues[1]["S1"], cgframe.residues[1]["S2"]] + atoms = measure.global_connections[0].get_atoms(cgframe) self.assertListEqual(target_bond, atoms) - atoms = measure.global_connections[1].get_atoms(frame) + atoms = measure.global_connections[1].get_atoms(cgframe) self.assertListEqual(target_angle, atoms) + def test_connect_residues(self): + DummyOptions.generate_angles = False + mapping = Mapping("test/data/global.map", DummyOptions) + measure = BondSet("test/data/global.bnd", DummyOptions) + frame = Frame("test/data/global-cg.gro") + measure.global_connections[0].eqm = 2. + measure.global_connections[0].fconst = 1000. + measure.global_connections[1].eqm = 90. + measure.global_connections[1].fconst = 100. + measure.connect_residues(frame, mapping) + self.assertEqual(measure["GLX_GLY"].beads[-1].num, 3) + self.assertListEqual(measure["GLX_GLY"].bonds[0].atom_numbers, [1, 2]) + + def test_global_itp(self): + DummyOptions.generate_angles = False + mapping = Mapping("test/data/global.map", DummyOptions) + measure = BondSet("test/data/global.bnd", DummyOptions) + frame = Frame("test/data/global-cg.gro") + measure.global_connections[0].eqm = 2. + measure.global_connections[0].fconst = 1000. + measure.global_connections[1].eqm = 3. + measure.global_connections[1].fconst = 100. + measure.connect_residues(frame, mapping) + measure.write_itp("global.itp", mapping) + diff --git a/test/test_util.py b/test/test_util.py index 84983a6..54a4ef3 100644 --- a/test/test_util.py +++ b/test/test_util.py @@ -11,7 +11,7 @@ from pycgtool.util import dir_up, backup_file, sliding, r_squared, dist_with_pbc from pycgtool.util import SimpleEnum, FixedFormatUnpacker from pycgtool.util import file_write_lines, cmp_whitespace_float -from pycgtool.util import circular_mean, circular_variance +from pycgtool.util import circular_mean, circular_variance, merge_list_of_lists class UtilTest(unittest.TestCase): @@ -272,6 +272,11 @@ def test_file_write_lines_append(self): with open(filename) as f: self.assertListEqual(lines, f.read().splitlines()) + def test_merge_list_of_lists(self): + l = [['a', 'b', 'c'], ['b', 'd', 'e'], ['k'], ['o', 'p'], ['e', 'f'], ['p', 'a'], ['d', 'g']] + merged_list = [['a','b','c','d','e','f','g','o','p'],['k']] + self.assertListEqual(merged_list,merge_list_of_lists(l)) + if __name__ == '__main__': unittest.main() From a9d67ad1f3d03a5f8255339add074a1c2c18be39 Mon Sep 17 00:00:00 2001 From: Jon Shearer Date: Sat, 22 Dec 2018 19:29:34 +0000 Subject: [PATCH 04/13] Updated documentation --- doc/source/index.rst | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/doc/source/index.rst b/doc/source/index.rst index d6b31bd..e100b5e 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -191,6 +191,16 @@ generate_angles Generate angles from bonds **True**, Fals generate_dihedrals Generate dihedrals from bonds **False**, True ================== ========================================== ======================= +Connect Residues +~~~~~~~~~~~~~~~~ +In the .bnd file the [ global ] section can be used to form bonds, angles or dihedrals between any residue in a system. +The syntax for a bond in the global section is:: + [atom_name1]_[resid1]_[resname1] [atom_name2]_[resid2]_[resname2] + or angle:: + [atom_name1]_[resid1]_[resname1] [atom_name2]_[resid2]_[resname2] [atom_name3]_[resid3]_[resname3] +Currently this feature only works for a single molecule and the name of the molecule outputted will always be 'mol_0'. +Also itp reading is not currently supported for multi-residue itps + Indexes ======= From a67d29e4e9ef3449741f2055e1c01f09d58909b9 Mon Sep 17 00:00:00 2001 From: jonathan Date: Fri, 18 Jan 2019 11:47:48 +0000 Subject: [PATCH 05/13] Fixed multiple bugs in the output of multi-residue information. At the moment only one multi-residue molecule is supported. --- pycgtool/bondset.py | 88 +++++++++++++++++++++++++++------------------ pycgtool/frame.py | 5 +-- 2 files changed, 57 insertions(+), 36 deletions(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 01758d7..27be9ec 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -451,7 +451,7 @@ def connect_residues(self, frame, mapping): beads = list(map(deepcopy, mapping[resname])) if len(self._molecules) > 0: bonds = list(map(deepcopy, self._molecules[resname])) - all_bonds.append(bonds) + all_bonds.extend(bonds) for i, bead in enumerate(beads): bead.num = start + i @@ -475,7 +475,7 @@ def connect_residues(self, frame, mapping): all_bonds.extend(self.global_connections) molecule = Molecule(resnames, all_bonds, all_beads) - molecule_name = "_".join(resnames) + molecule_name = "mol_01" self._molecules[molecule_name] = molecule else: return @@ -540,6 +540,7 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= ret_lines.append(atom_template["mass"].format( i, bead.type, 1, mol, bead.name, i, bead.charge, bead.mass )) + else: ret_lines.append(atom_template["nomass"].format( i, bead.type, 1, mol, bead.name, i, bead.charge @@ -574,35 +575,39 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= logger.warning(("Residue '{0}' has no internal bonds").format(res)) ret_lines.append("\n[ moleculetype ]") - mol_name = "mol_" + str(mol_count) - ret_lines.append("{0:4s} {1:4d}".format(mol_name, 1)) - ret_lines.append("\n[ atoms ]") + ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) - for resid, beads in self._molecules[mol].resid_to_beads.items(): + ret_lines.append("\n[ atoms ]") + for resid in range(1, len(self._molecules[mol].resnames) + 1): + beads = self._molecules[mol].resid_to_beads[resid] + resname = self._molecules[mol].resid_to_resname[resid] for bead in beads: # atnum type resnum resname atname c-group charge (mass) if isinstance(bead, VirtualMap): ret_lines.append(atom_template["mass"].format( - bead.num + 1, bead.type, resid, mol, bead.name, bead.num + 1, bead.charge, bead.mass + bead.num + 1, bead.type, resid, resname, bead.name, bead.num + 1, bead.charge, bead.mass )) else: ret_lines.append(atom_template["nomass"].format( - bead.num + 1, bead.type, resid, mol, bead.name, bead.num + 1, bead.charge + bead.num + 1, bead.type, resid, resname, bead.name, bead.num + 1, bead.charge )) virtual_beads = self.get_virtual_beads(self._molecules[mol].beads) if len(virtual_beads) != 0: ret_lines.append("\n[ virtual_sitesn ]") excl_lines = ["\n[ exclusions ]"] # exlusions section for virtual sites - for vbead in virtual_beads: - CGids = [bead.num + 1 for bead in self._molecules[mol].beads if bead.name in vbead.atoms] - CGids.sort() - CGids_string = " ".join(map(str, CGids)) - ret_lines.append( - "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) - vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string - excl_lines.append(vsite_exclusions) + for resid in range(1, len(self._molecules[mol].resnames) + 1): + beads = self._molecules[mol].resid_to_beads[resid] + virtual_beads = self.get_virtual_beads(beads) + for vbead in virtual_beads: + CGids = [bead.num + 1 for bead in beads if bead.name in vbead.atoms] + CGids.sort() + CGids_string = " ".join(map(str, CGids)) + ret_lines.append( + "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) + vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string + excl_lines.append(vsite_exclusions) ret_lines.extend(excl_lines) ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) @@ -721,26 +726,41 @@ def dump_values(self, target_number=10000): for mol in self._molecules: if mol == "SOL": continue - if isinstance(self._molecules[mol], Molecule): - continue - bonds = self.get_bond_lengths(mol, with_constr=True) - if bonds: - lines = BondSet._get_lines_for_bond_dump(bonds, target_number) - file_write_lines("{0}_length.dat".format(mol), lines) + if not isinstance(self._molecules[mol], Molecule): + bonds = self.get_bond_lengths(mol, with_constr=True) + if bonds: + lines = BondSet._get_lines_for_bond_dump(bonds, target_number) + file_write_lines("{0}_length.dat".format(mol), lines) + + bonds = self.get_bond_angles(mol) + if bonds: + lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) + file_write_lines("{0}_angle.dat".format(mol), lines) + + bonds = self.get_bond_dihedrals(mol) + if bonds: + lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) + file_write_lines("{0}_dihedral.dat".format(mol), lines) + else: + bonds = self.get_bond_lengths(mol, with_constr=True) + global_bonds = [bond for bond in bonds if isinstance(bond, GlobalBond)] + if global_bonds: + lines = BondSet._get_lines_for_bond_dump(global_bonds, target_number) + file_write_lines("{0}_length.dat".format(mol), lines) + + bonds = self.get_bond_angles(mol) + global_bonds = [bond for bond in bonds if isinstance(bond, GlobalBond)] + if global_bonds: + lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) + file_write_lines("{0}_angle.dat".format(mol), lines) + + bonds = self.get_bond_dihedrals(mol) + global_bonds = [bond for bond in bonds if isinstance(bond, GlobalBond)] + if global_bonds: + lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) + file_write_lines("{0}_dihedral.dat".format(mol), lines) - bonds = self.get_bond_angles(mol) - if bonds: - lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) - file_write_lines("{0}_angle.dat".format(mol), lines) - bonds = self.get_bond_dihedrals(mol) - if bonds: - lines = BondSet._get_lines_for_bond_dump(bonds, target_number, rad2deg=True) - file_write_lines("{0}_dihedral.dat".format(mol), lines) - if len(self.global_connections) > 0: - global_bonds = self.global_connections - lines = BondSet._get_lines_for_bond_dump(global_bonds, target_number) - file_write_lines("global_length.dat".format(mol), lines) def __len__(self): return len(self._molecules) diff --git a/pycgtool/frame.py b/pycgtool/frame.py index ea9e2d3..6f8130a 100644 --- a/pycgtool/frame.py +++ b/pycgtool/frame.py @@ -70,11 +70,12 @@ class Molecule: """ Holds data for a molecule comprised of multiple residues """ - __slots__ = ["resnames", "bonds", "beads", "resid_to_beads"] + __slots__ = ["resnames", "bonds", "beads", "resid_to_beads", "resid_to_resname"] def __init__(self, resnames, bonds, beads): self.resnames = resnames - self.resid_to_beads = dict(zip(range(len(resnames)), beads)) + self.resid_to_resname = dict(zip(range(1, len(resnames)+1), resnames)) + self.resid_to_beads = dict(zip(range(1, len(resnames)+1), beads)) self.bonds = bonds self.beads = np.array(beads).flatten().tolist() From 2dfa5cde6302abed9b430a3b7500fab8f531fbb0 Mon Sep 17 00:00:00 2001 From: jonathan Date: Fri, 18 Jan 2019 11:54:30 +0000 Subject: [PATCH 06/13] updated tests --- test/test_bondset.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_bondset.py b/test/test_bondset.py index 56877aa..443d446 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -309,8 +309,8 @@ def test_connect_residues(self): measure.global_connections[1].eqm = 90. measure.global_connections[1].fconst = 100. measure.connect_residues(frame, mapping) - self.assertEqual(measure["GLX_GLY"].beads[-1].num, 3) - self.assertListEqual(measure["GLX_GLY"].bonds[0].atom_numbers, [1, 2]) + self.assertEqual(measure["mol_01"].beads[-1].num, 3) + self.assertListEqual(measure["mol_01"].bonds[0].atom_numbers, [1, 2]) def test_global_itp(self): DummyOptions.generate_angles = False From e51667c43fc5bf628764ee7cfbce075f372738ad Mon Sep 17 00:00:00 2001 From: jonathan Date: Mon, 21 Jan 2019 13:55:24 +0000 Subject: [PATCH 07/13] Global section removed. Now can enter inter residue connections under a molecule section [< mol >] and the name used for mol will apper in the itp file. --- pycgtool/bondset.py | 155 ++++++++++++++++++++++++------------------- pycgtool/frame.py | 16 +++-- test/test_bondset.py | 37 ++++++----- 3 files changed, 121 insertions(+), 87 deletions(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 27be9ec..7b3019f 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -184,7 +184,6 @@ def __init__(self, filename, options): :return: Instance of BondSet """ self._molecules = {} - self.global_connections = [] self._fconst_constr_threshold = options.constr_threshold @@ -222,11 +221,19 @@ def __init__(self, filename, options): except AttributeError: pass + molecule_delimiters = ('<', '>') + is_global = False with CFG(filename) as cfg: for mol_name, mol_section in cfg.items(): - is_global = False - if mol_name == "global": + if "<" in mol_name: is_global = True + for mol_name, mol_section in cfg.items(): + now_global = False + if "<" in mol_name: + now_global = True + mol_name = mol_name.strip(" ".join(molecule_delimiters)) + self._molecules[mol_name] = [] + mol_bonds = self._molecules[mol_name] else: self._molecules[mol_name] = [] mol_bonds = self._molecules[mol_name] @@ -254,8 +261,8 @@ def __init__(self, filename, options): if {x for x in atomlist if atomlist.count(x) > 1}: raise ValueError("Defined bond '{0}' contains duplicate atoms".format(atomlist)) - if is_global: - self.global_connections.append(GlobalBond(atoms=atomlist, func_form=func_form)) + if now_global: + mol_bonds.append(GlobalBond(atoms=atomlist, func_form=func_form)) else: mol_bonds.append(Bond(atoms=atomlist, func_form=func_form)) if len(atomlist) > 2: @@ -265,8 +272,8 @@ def __init__(self, filename, options): if is_global: if options.generate_angles or options.generate_dihedrals: logger.warning("Automated generation of angles or dihedrals between " - "residues not implemented! Please specifiy angles and dihedrals in [global] " - "section of *.bnd file") + "residues not implemented! Please specifiy angles and dihedrals in [< {0} >] " + "section of *.bnd file".format(mol_name)) else: angles, dihedrals = self._create_angles(mol_bonds) @@ -288,6 +295,8 @@ def __init__(self, filename, options): func_form=self._functional_forms[4](circular_mean, circular_variance) ) ) + if now_global: + self._molecules[mol_name] = Molecule(bonds=mol_bonds) @staticmethod def _create_angles(mol_bonds): @@ -430,55 +439,61 @@ def connect_residues(self, frame, mapping): connects residues together to form new molecules :param frame: Frame from which to determine inter residue connections """ - if len(self.global_connections) > 0: - fragments = [] - for bond in self.global_connections: - fragments.append(bond.get_residue_ids(frame)) - - self._populate_atom_numbers(mapping) - - molecule_internal_ids = merge_list_of_lists(fragments) - if len(molecule_internal_ids) > 1: - print("More than one multi residue molecule detected - this is not yet supported by pycgtool") - raise NotImplementedError - molecule_internal_ids = molecule_internal_ids[0] - resnames = [frame[mol_id].name for mol_id in molecule_internal_ids] - all_bonds = [] - all_beads = [] - start = 0 - for resid, resname in enumerate(resnames): - if resname in mapping: - beads = list(map(deepcopy, mapping[resname])) - if len(self._molecules) > 0: - bonds = list(map(deepcopy, self._molecules[resname])) - all_bonds.extend(bonds) - for i, bead in enumerate(beads): - bead.num = start + i - - if len(self._molecules) > 0: - index = [bead.name for bead in beads] - for bond in bonds: - try: - bond.atom_numbers = [index.index(atom.lstrip("+-")) for atom in bond.atoms] - except ValueError as e: - missing = [atom for atom in bond.atoms if atom.lstrip("+-") not in index] - e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, resname),) - raise - bond.atom_numbers = [start + ind for ind in bond.atom_numbers] - - all_beads.append(beads) - start = beads[-1].num + 1 - - - for bond in self.global_connections: - bond.populate_ids(all_beads) - - all_bonds.extend(self.global_connections) - molecule = Molecule(resnames, all_bonds, all_beads) - molecule_name = "mol_01" - self._molecules[molecule_name] = molecule - else: - return + count = 0 + not_multi = [mol for mol in self._molecules if not isinstance(self._molecules[mol], Molecule)] + for mol in self._molecules: + if isinstance(self._molecules[mol], Molecule): + if count > 1: + print("More than one multi residue molecule detected - this is not yet supported by pycgtool") + raise NotImplementedError + mol_bonds = self._molecules[mol].bonds + fragments = [] + for bond in mol_bonds: + fragments.append(bond.get_residue_ids(frame)) + + self._populate_atom_numbers(mapping) + + molecule_internal_ids = merge_list_of_lists(fragments) + if len(molecule_internal_ids) > 1: + print("All Fragments of Molecule '{}' are not connected - add missing bonds to .bnd".format(mol)) + raise SyntaxError + molecule_internal_ids = molecule_internal_ids[0] + resnames = [frame[mol_id].name for mol_id in molecule_internal_ids] + all_bonds = [] + all_beads = [] + start = 0 + for resid, resname in enumerate(resnames): + if resname in mapping: + beads = list(map(deepcopy, mapping[resname])) + if len(not_multi) > 0: + bonds = list(map(deepcopy, self._molecules[resname])) + all_bonds.extend(bonds) + for i, bead in enumerate(beads): + bead.num = start + i + + if len(not_multi) > 0: + index = [bead.name for bead in beads] + for bond in bonds: + try: + bond.atom_numbers = [index.index(atom.lstrip("+-")) for atom in bond.atoms] + except ValueError as e: + missing = [atom for atom in bond.atoms if atom.lstrip("+-") not in index] + e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, resname),) + raise + bond.atom_numbers = [start + ind for ind in bond.atom_numbers] + + all_beads.append(beads) + start = beads[-1].num + 1 + + + for bond in mol_bonds: + bond.populate_ids(all_beads) + + all_bonds.extend(mol_bonds) + molecule = Molecule(resnames, all_bonds, all_beads) + self._molecules[mol] = molecule + count += 1 + return def write_itp(self, filename, mapping): """ @@ -668,17 +683,21 @@ def calc_dihedral(atoms): pass except ZeroDivisionError as e: e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.atoms)),) + #print(type(bond), mol_meas) raise e + #calculate values for global bonds# - if len(self.global_connections) > 0: - for bond in self.global_connections: - try: - atoms = bond.get_atoms(frame) - val = calc[len(atoms)](atoms) - bond.values.append(val) - except ZeroDivisionError as e: - e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.labels)),) - raise e + for mol in self._molecules: + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + for bond in bonds: + try: + atoms = bond.get_atoms(frame) + val = calc[len(atoms)](atoms) + bond.values.append(val) + except ZeroDivisionError as e: + e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.labels)),) + raise e def boltzmann_invert(self, progress=False): """ @@ -687,8 +706,10 @@ def boltzmann_invert(self, progress=False): :param progress: Display a progress bar using tqdm if available """ bond_iter = itertools.chain(*self._molecules.values()) - if len(self.global_connections) > 0: - bond_iter = itertools.chain(bond_iter, self.global_connections) + for mol in self._molecules: + if isinstance(self._molecules[mol], Molecule): + bonds = self._molecules[mol].bonds + bond_iter = itertools.chain(bond_iter, bonds) if progress: total = sum(map(len, self._molecules.values())) bond_iter = tqdm(bond_iter, total=total, ncols=80) diff --git a/pycgtool/frame.py b/pycgtool/frame.py index 6f8130a..7e786b8 100644 --- a/pycgtool/frame.py +++ b/pycgtool/frame.py @@ -72,12 +72,20 @@ class Molecule: """ __slots__ = ["resnames", "bonds", "beads", "resid_to_beads", "resid_to_resname"] - def __init__(self, resnames, bonds, beads): + def __init__(self, resnames=None, bonds=None, beads=None): self.resnames = resnames - self.resid_to_resname = dict(zip(range(1, len(resnames)+1), resnames)) - self.resid_to_beads = dict(zip(range(1, len(resnames)+1), beads)) + if beads and resnames is not None: + self.resid_to_resname = dict(zip(range(1, len(resnames)+1), resnames)) + self.resid_to_beads = dict(zip(range(1, len(resnames)+1), beads)) + self.beads = np.array(beads).flatten().tolist() self.bonds = bonds - self.beads = np.array(beads).flatten().tolist() + + def __len__(self): + return len(self.bonds) + + def __iter__(self): + return iter(self.bonds) + class Residue: diff --git a/test/test_bondset.py b/test/test_bondset.py index 443d446..67f82fb 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -278,49 +278,54 @@ def test_get_lines_for_bond_dump_sample(self): seen.add(line) def test_global_bond_create(self): + mol = "mol_01" DummyOptions.generate_angles=False measure = BondSet("test/data/global.bnd", DummyOptions) - self.assertEqual(len(measure.global_connections), 2) + self.assertEqual(len(measure[mol].bonds), 2) names = ["S2", "S1"] resids = [1, 2] resnames = ["GLX", "GLY"] - self.assertListEqual(resnames, measure.global_connections[0].resnames) - self.assertListEqual(resids, measure.global_connections[0].resids) - self.assertListEqual(names, measure.global_connections[0].atoms) + self.assertListEqual(resnames, measure[mol].bonds[0].resnames) + self.assertListEqual(resids, measure[mol].bonds[0].resids) + self.assertListEqual(names, measure[mol].bonds[0].atoms) def test_global_bond_get_atoms(self): + mol = "mol_01" DummyOptions.generate_angles = False measure = BondSet("test/data/global.bnd", DummyOptions) cgframe = Frame("test/data/global-cg.gro") target_bond = [cgframe.residues[0]["S2"], cgframe.residues[1]["S1"]] target_angle = [cgframe.residues[0]["S2"], cgframe.residues[1]["S1"], cgframe.residues[1]["S2"]] - atoms = measure.global_connections[0].get_atoms(cgframe) + atoms = measure[mol].bonds[0].get_atoms(cgframe) self.assertListEqual(target_bond, atoms) - atoms = measure.global_connections[1].get_atoms(cgframe) + atoms = measure[mol].bonds[1].get_atoms(cgframe) self.assertListEqual(target_angle, atoms) def test_connect_residues(self): + mol = "mol_01" DummyOptions.generate_angles = False mapping = Mapping("test/data/global.map", DummyOptions) measure = BondSet("test/data/global.bnd", DummyOptions) frame = Frame("test/data/global-cg.gro") - measure.global_connections[0].eqm = 2. - measure.global_connections[0].fconst = 1000. - measure.global_connections[1].eqm = 90. - measure.global_connections[1].fconst = 100. + measure[mol].bonds[0].eqm = 2. + measure[mol].bonds[0].fconst = 1000. + measure[mol].bonds[1].eqm = 90. + measure[mol].bonds[1].fconst = 100. measure.connect_residues(frame, mapping) - self.assertEqual(measure["mol_01"].beads[-1].num, 3) - self.assertListEqual(measure["mol_01"].bonds[0].atom_numbers, [1, 2]) + self.assertEqual(measure[mol].beads[-1].num, 3) + self.assertListEqual(measure[mol].bonds[0].atom_numbers, [1, 2]) def test_global_itp(self): + mol = "mol_01" DummyOptions.generate_angles = False mapping = Mapping("test/data/global.map", DummyOptions) measure = BondSet("test/data/global.bnd", DummyOptions) frame = Frame("test/data/global-cg.gro") - measure.global_connections[0].eqm = 2. - measure.global_connections[0].fconst = 1000. - measure.global_connections[1].eqm = 3. - measure.global_connections[1].fconst = 100. + print(measure) + measure[mol].bonds[0].eqm = 2. + measure[mol].bonds[0].fconst = 1000. + measure[mol].bonds[1].eqm = 3. + measure[mol].bonds[1].fconst = 100. measure.connect_residues(frame, mapping) measure.write_itp("global.itp", mapping) From 4ad9e903707b7a80ca4fd20f49b4bf339995a864 Mon Sep 17 00:00:00 2001 From: jonathan Date: Mon, 21 Jan 2019 13:57:12 +0000 Subject: [PATCH 08/13] updated inter residue test files --- test/data/global.bnd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/data/global.bnd b/test/data/global.bnd index cab2579..7625925 100644 --- a/test/data/global.bnd +++ b/test/data/global.bnd @@ -1,5 +1,5 @@ ; comments begin with a semicolon -[ global ] +[< mol_01 >] S2_1_GLX S1_2_GLY S2_1_GLX S1_2_GLY S2_2_GLY From 87820c354fba887c883fc0e628e6755f87999a23 Mon Sep 17 00:00:00 2001 From: jonathan Date: Fri, 22 Feb 2019 11:01:48 +0000 Subject: [PATCH 09/13] updated all bonset workflow to go via Molecule objects. Added test for multi residue molecule itp output --- pycgtool/bondset.py | 141 +++++++++++++++++++++++++++++-------------- pycgtool/frame.py | 18 +----- test/data/global.itp | 21 +++++++ test/test_bondset.py | 9 ++- 4 files changed, 124 insertions(+), 65 deletions(-) create mode 100644 test/data/global.itp diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 7b3019f..3ac31c2 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -18,7 +18,6 @@ from .mapping import VirtualMap from .functionalforms import FunctionalForms -from .frame import Molecule from .parsers.cfg import CFG from .util import ( circular_mean, @@ -38,6 +37,50 @@ logger = logging.getLogger(__name__) +class Molecule: + """ + Holds data for a molecule comprised of multiple residues + """ + __slots__ = ["resnames", "bonds", "beads", "resid_to_beads", "resid_to_resname"] + + def __init__(self, resnames=None, bonds=None, beads=None): + self.resnames = resnames + if beads and resnames is not None: + self.resid_to_resname = dict(zip(range(1, len(resnames) + 1), resnames)) + self.resid_to_beads = dict(zip(range(1, len(resnames) + 1), beads)) + self.beads = np.array(beads).flatten().tolist() + + self.bonds = bonds + if bonds is None: + self.bonds = [] + + + def add_bond(self, bond): + self.bonds.append(bond) + + @property + def inter(self): + '''collects inter residue bonds''' + inter = [] + for bond in self.bonds: + if isinstance(bond, GlobalBond): + inter.append(bond) + return inter + + @property + def is_multiresidue(self): + '''checks if Molecule has multiple residues''' + return True if len(self.inter) > 0 else False + + def __len__(self): + return len(self.bonds) + + def __iter__(self): + return iter(self.bonds) + + def __getitem__(self, item): + return self.bonds[item] + class Bond: """ Class holding the properties of a single bonded term. @@ -223,6 +266,7 @@ def __init__(self, filename, options): molecule_delimiters = ('<', '>') is_global = False + multi_count = 0 with CFG(filename) as cfg: for mol_name, mol_section in cfg.items(): if "<" in mol_name: @@ -230,12 +274,18 @@ def __init__(self, filename, options): for mol_name, mol_section in cfg.items(): now_global = False if "<" in mol_name: + multi_count += 1 + if multi_count > 1: + print("More than one multi residue molecule detected - this is not yet supported by pycgtool") + raise NotImplementedError now_global = True mol_name = mol_name.strip(" ".join(molecule_delimiters)) - self._molecules[mol_name] = [] + self._molecules[mol_name] = Molecule() + mol_bonds = self._molecules[mol_name] + else: - self._molecules[mol_name] = [] + self._molecules[mol_name] = Molecule() mol_bonds = self._molecules[mol_name] angles_defined = False @@ -262,14 +312,14 @@ def __init__(self, filename, options): if {x for x in atomlist if atomlist.count(x) > 1}: raise ValueError("Defined bond '{0}' contains duplicate atoms".format(atomlist)) if now_global: - mol_bonds.append(GlobalBond(atoms=atomlist, func_form=func_form)) + mol_bonds.add_bond(GlobalBond(atoms=atomlist, func_form=func_form)) else: - mol_bonds.append(Bond(atoms=atomlist, func_form=func_form)) + mol_bonds.add_bond(Bond(atoms=atomlist, func_form=func_form)) if len(atomlist) > 2: angles_defined = True if not angles_defined: - if is_global: + if now_global: if options.generate_angles or options.generate_dihedrals: logger.warning("Automated generation of angles or dihedrals between " "residues not implemented! Please specifiy angles and dihedrals in [< {0} >] " @@ -280,7 +330,7 @@ def __init__(self, filename, options): if options.generate_angles: for atomlist in angles: - mol_bonds.append( + mol_bonds.add_bond( Bond( atoms=atomlist, func_form=self._functional_forms[3](circular_mean, circular_variance) @@ -289,14 +339,13 @@ def __init__(self, filename, options): if options.generate_dihedrals: for atomlist in dihedrals: - mol_bonds.append( + mol_bonds.add_bond( Bond( atoms=atomlist, func_form=self._functional_forms[4](circular_mean, circular_variance) ) ) - if now_global: - self._molecules[mol_name] = Molecule(bonds=mol_bonds) + @staticmethod def _create_angles(mol_bonds): @@ -320,10 +369,7 @@ def get_bonds(self, mol, natoms, select=lambda x: True): :param function select: Optional lambda, return only bonds for which this is True :return List[Bond]: List of bonds """ - if isinstance(self._molecules[mol], Molecule): - bonds = self._molecules[mol].bonds - else: - bonds = self._molecules[mol] + bonds = self._molecules[mol] if natoms == -1: return [bond for bond in bonds if select(bond)] @@ -337,10 +383,7 @@ def get_bond_lengths(self, mol, with_constr=False): :param with_constr: Include constraints? :return: List of bonds """ - if isinstance(self._molecules[mol], Molecule): - bonds = self._molecules[mol].bonds - else: - bonds = self._molecules[mol] + bonds = self._molecules[mol] if with_constr: return [bond for bond in bonds if len(bond.atoms) == 2] @@ -354,10 +397,7 @@ def get_bond_length_constraints(self, mol): :param mol: Molecule name to return bonds for :return: List of bonds """ - if isinstance(self._molecules[mol], Molecule): - bonds = self._molecules[mol].bonds - else: - bonds = self._molecules[mol] + bonds = self._molecules[mol] return [bond for bond in bonds if len(bond.atoms) == 2 and bond.fconst >= self._fconst_constr_threshold] def get_bond_angles(self, mol, exclude_triangle=True): @@ -368,10 +408,7 @@ def get_bond_angles(self, mol, exclude_triangle=True): :param exclude_triangle: Exclude angles that are part of a triangle? :return: List of bonds """ - if isinstance(self._molecules[mol], Molecule): - bonds = self._molecules[mol].bonds - else: - bonds = self._molecules[mol] + bonds = self._molecules[mol] angles = [bond for bond in bonds if len(bond.atoms) == 3] if exclude_triangle: @@ -395,10 +432,7 @@ def get_bond_dihedrals(self, mol): :param mol: Molecule name to return bonds for :return: List of bonds """ - if isinstance(self._molecules[mol], Molecule): - bonds = self._molecules[mol].bonds - else: - bonds = self._molecules[mol] + bonds = self._molecules[mol] return [bond for bond in bonds if len(bond.atoms) == 4] def get_virtual_beads(self, mapping): @@ -439,14 +473,11 @@ def connect_residues(self, frame, mapping): connects residues together to form new molecules :param frame: Frame from which to determine inter residue connections """ - count = 0 - not_multi = [mol for mol in self._molecules if not isinstance(self._molecules[mol], Molecule)] + #TODO this section is a bit verbose - simplify + not_multi = [mol for mol in self._molecules if not self._molecules[mol].is_multiresidue] for mol in self._molecules: - if isinstance(self._molecules[mol], Molecule): - if count > 1: - print("More than one multi residue molecule detected - this is not yet supported by pycgtool") - raise NotImplementedError - mol_bonds = self._molecules[mol].bonds + if self._molecules[mol].is_multiresidue: + mol_bonds = self._molecules[mol].inter fragments = [] for bond in mol_bonds: fragments.append(bond.get_residue_ids(frame)) @@ -489,10 +520,19 @@ def connect_residues(self, frame, mapping): for bond in mol_bonds: bond.populate_ids(all_beads) + for mol in self._molecules: + num = 0 + for res in frame: + if res.name == mol: + num += 1 + if num > 1: + logger.warning("More than one {0} residue found in gro, note that pycgtool " + "does not currently support the fitting of bonding parameters of multi-residue " + "molecules + single residue molecules simultaneously.") + all_bonds.extend(mol_bonds) molecule = Molecule(resnames, all_bonds, all_beads) self._molecules[mol] = molecule - count += 1 return def write_itp(self, filename, mapping): @@ -537,7 +577,7 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= # Print molecule not_calc = " Parameters have not been calculated." for mol in self._molecules: - if not isinstance(self._molecules[mol], Molecule): + if not self._molecules[mol].is_multiresidue: if mol not in mapping: logger.warning("Molecule '{0}' present in bonding file, but not in mapping.".format(mol) + not_calc) continue @@ -574,6 +614,7 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= excl_lines.append(vsite_exclusions) ret_lines.extend(excl_lines) + ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) @@ -662,6 +703,7 @@ def calc_dihedral(atoms): 3: calc_angle, 4: calc_dihedral} + # TODO tidy this section for prev_res, res, next_res in sliding(frame): try: mol_meas = self._molecules[res.name] @@ -674,26 +716,33 @@ def calc_dihedral(atoms): for bond in mol_meas: try: - # TODO tidy this atoms = [adj_res.get(name[0], res)[name.lstrip("-+")] for name in bond.atoms] - val = calc[len(atoms)](atoms) + num_atoms = len(atoms) + if num_atoms == 0: + raise Exception("Error: {0} residue in *.map and *.gro file, but atoms {1} not in gro".format( + res.name, " ".join(bond.atoms))) + val = calc[num_atoms](atoms) bond.values.append(val) except (NotImplementedError, TypeError): # TypeError is raised when residues on end of chain calc bond to next pass except ZeroDivisionError as e: e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.atoms)),) - #print(type(bond), mol_meas) raise e - #calculate values for global bonds# + #calculate values for global bonds for mol in self._molecules: - if isinstance(self._molecules[mol], Molecule): + if self._molecules[mol].is_multiresidue: bonds = self._molecules[mol].bonds for bond in bonds: try: atoms = bond.get_atoms(frame) - val = calc[len(atoms)](atoms) + num_atoms = len(atoms) + if num_atoms == 0: + raise Exception( + "Error: {0} Molecule in *.map, but atoms {1} not in gro".format( + mol, " ".join(bond.atoms))) + val = calc[num_atoms](atoms) bond.values.append(val) except ZeroDivisionError as e: e.args = ("Zero division in calculation of <{0}>".format(" ".join(bond.labels)),) @@ -707,7 +756,7 @@ def boltzmann_invert(self, progress=False): """ bond_iter = itertools.chain(*self._molecules.values()) for mol in self._molecules: - if isinstance(self._molecules[mol], Molecule): + if self._molecules[mol].is_multiresidue: bonds = self._molecules[mol].bonds bond_iter = itertools.chain(bond_iter, bonds) if progress: diff --git a/pycgtool/frame.py b/pycgtool/frame.py index 7e786b8..30dfee3 100644 --- a/pycgtool/frame.py +++ b/pycgtool/frame.py @@ -8,6 +8,7 @@ import logging import numpy as np +import collections from .util import backup_file, file_write_lines from .parsers.cfg import CFG @@ -66,25 +67,8 @@ def add_missing_data(self, other): setattr(self, attr, getattr(other, attr)) -class Molecule: - """ - Holds data for a molecule comprised of multiple residues - """ - __slots__ = ["resnames", "bonds", "beads", "resid_to_beads", "resid_to_resname"] - def __init__(self, resnames=None, bonds=None, beads=None): - self.resnames = resnames - if beads and resnames is not None: - self.resid_to_resname = dict(zip(range(1, len(resnames)+1), resnames)) - self.resid_to_beads = dict(zip(range(1, len(resnames)+1), beads)) - self.beads = np.array(beads).flatten().tolist() - self.bonds = bonds - def __len__(self): - return len(self.bonds) - - def __iter__(self): - return iter(self.bonds) diff --git a/test/data/global.itp b/test/data/global.itp new file mode 100644 index 0000000..e13291f --- /dev/null +++ b/test/data/global.itp @@ -0,0 +1,21 @@ +; +; Topology prepared automatically using PyCGTOOL +; James Graham 2016 +; University of Southampton +; https://github.com/jag1g13/pycgtool +; + +[ moleculetype ] +mol_01 1 + +[ atoms ] + 1 TC3 1 GLX S1 1 0.000 + 2 TC3 1 GLX S2 2 0.000 + 3 TC3 2 GLY S1 3 0.000 + 4 TC3 2 GLY S2 4 0.000 + +[ bonds ] + 2 3 1 2.00000 1000.00000 + +[ angles ] + 2 3 4 2 30.00000 100.00000 diff --git a/test/test_bondset.py b/test/test_bondset.py index 67f82fb..ea10f4f 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -3,6 +3,7 @@ import logging import math import os +import numpy as np from pycgtool.bondset import BondSet from pycgtool.frame import Frame @@ -321,12 +322,16 @@ def test_global_itp(self): mapping = Mapping("test/data/global.map", DummyOptions) measure = BondSet("test/data/global.bnd", DummyOptions) frame = Frame("test/data/global-cg.gro") - print(measure) measure[mol].bonds[0].eqm = 2. measure[mol].bonds[0].fconst = 1000. - measure[mol].bonds[1].eqm = 3. + measure[mol].bonds[1].eqm = np.deg2rad(30.) measure[mol].bonds[1].fconst = 100. measure.connect_residues(frame, mapping) measure.write_itp("global.itp", mapping) + logging.disable(logging.NOTSET) + + self.assertTrue(cmp_file_whitespace_float("global.itp", "test/data/global.itp", + rtol=0.005, verbose=True)) + From 0ff2ca4f50897c0b326441e5c7b6bab253371867 Mon Sep 17 00:00:00 2001 From: jonathan Date: Fri, 22 Feb 2019 11:20:12 +0000 Subject: [PATCH 10/13] fixed bug in previous commit --- pycgtool/bondset.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 3ac31c2..0d7392c 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -796,7 +796,8 @@ def dump_values(self, target_number=10000): for mol in self._molecules: if mol == "SOL": continue - if not isinstance(self._molecules[mol], Molecule): + if not self._molecules[mol].is_multiresidue: + bonds = self.get_bond_lengths(mol, with_constr=True) if bonds: lines = BondSet._get_lines_for_bond_dump(bonds, target_number) From bec6af4815abcf239cd20756396e2441d6cf43d8 Mon Sep 17 00:00:00 2001 From: Jon Shearer Date: Fri, 22 Feb 2019 18:35:32 +0000 Subject: [PATCH 11/13] refactored output to make more streamlined for molecules with any number of residues --- pycgtool/bondset.py | 136 +++++++++++++++++--------------------------- 1 file changed, 53 insertions(+), 83 deletions(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 0d7392c..0479856 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -41,14 +41,19 @@ class Molecule: """ Holds data for a molecule comprised of multiple residues """ - __slots__ = ["resnames", "bonds", "beads", "resid_to_beads", "resid_to_resname"] + __slots__ = ["resnames", "bonds", "beads", "resname_to_beads"] def __init__(self, resnames=None, bonds=None, beads=None): + ''' + + :param resnames: list of resname names + :param bonds: list of lists of BeadMap objects in the same order as 'resnames' + :param beads: list of bonds + ''' self.resnames = resnames if beads and resnames is not None: - self.resid_to_resname = dict(zip(range(1, len(resnames) + 1), resnames)) - self.resid_to_beads = dict(zip(range(1, len(resnames) + 1), beads)) - self.beads = np.array(beads).flatten().tolist() + self.resname_to_beads = dict(zip(resnames, beads)) + self.beads = [bead for res_beads in beads for bead in res_beads] self.bonds = bonds if bonds is None: @@ -56,6 +61,11 @@ def __init__(self, resnames=None, bonds=None, beads=None): def add_bond(self, bond): + ''' + Add bond to object + + :param bond: instance of class with 'Bond' as their base class + ''' self.bonds.append(bond) @property @@ -285,7 +295,7 @@ def __init__(self, filename, options): mol_bonds = self._molecules[mol_name] else: - self._molecules[mol_name] = Molecule() + self._molecules[mol_name] = Molecule(resnames=[mol_name]) mol_bonds = self._molecules[mol_name] angles_defined = False @@ -319,7 +329,7 @@ def __init__(self, filename, options): angles_defined = True if not angles_defined: - if now_global: + if is_global: if options.generate_angles or options.generate_dihedrals: logger.warning("Automated generation of angles or dihedrals between " "residues not implemented! Please specifiy angles and dihedrals in [< {0} >] " @@ -577,102 +587,62 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= # Print molecule not_calc = " Parameters have not been calculated." for mol in self._molecules: - if not self._molecules[mol].is_multiresidue: - if mol not in mapping: - logger.warning("Molecule '{0}' present in bonding file, but not in mapping.".format(mol) + not_calc) - continue - if not all((bond.fconst is not None for bond in self._molecules[mol])): - logger.warning("Molecule '{0}' has no measured bond values.".format(mol) + not_calc) - continue + molecule = self._molecules[mol] + + if not all((bond.fconst is not None for bond in self._molecules[mol])): + logger.warning("Molecule '{0}' has no measured bond values.".format(mol) + not_calc) + continue + + ret_lines.append("\n[ moleculetype ]") + ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) + + ret_lines.append("\n[ atoms ]") + + for resid, res in enumerate(molecule.resnames, start=1): + beads = molecule.resname_to_beads[res] if molecule.is_multiresidue else mapping[mol] - ret_lines.append("\n[ moleculetype ]") - ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) - ret_lines.append("\n[ atoms ]") - for i, bead in enumerate(mapping[mol], start=1): + if res not in mapping: + logger.warning("Residue '{0}' present in bonding file, but not in mapping.".format(mol) + not_calc) + continue + + for i, bead in enumerate(beads, start=1): # atnum type resnum resname atname c-group charge (mass) if isinstance(bead, VirtualMap): ret_lines.append(atom_template["mass"].format( - i, bead.type, 1, mol, bead.name, i, bead.charge, bead.mass + bead.num + 1, bead.type, resid, res, bead.name, bead.num + 1, bead.charge, bead.mass )) else: ret_lines.append(atom_template["nomass"].format( - i, bead.type, 1, mol, bead.name, i, bead.charge + bead.num + 1, bead.type, resid, res, bead.name, bead.num + 1, bead.charge )) - virtual_beads = self.get_virtual_beads(mapping[mol]) - if len(virtual_beads) != 0: - ret_lines.append("\n[ virtual_sitesn ]") - excl_lines = ["\n[ exclusions ]"] #exlusions section for virtual sites + all_beads = molecule.beads if molecule.is_multiresidue else mapping[mol] + virtual_beads = self.get_virtual_beads(all_beads) + if len(virtual_beads) != 0: + ret_lines.append("\n[ virtual_sitesn ]") + excl_lines = ["\n[ exclusions ]"] #exlusions section for virtual sites + for res in molecule.resnames: + beads = molecule.resname_to_beads[res] if molecule.is_multiresidue else mapping[mol] + virtual_beads = self.get_virtual_beads(beads) for vbead in virtual_beads: - CGids = [bead.num + 1 for bead in mapping[mol] if bead.name in vbead.atoms] + CGids = [bead.num + 1 for bead in beads if bead.name in vbead.atoms] CGids.sort() CGids_string = " ".join(map(str, CGids)) - ret_lines.append("{0:^6d} {1:^6d} {2}".format(vbead.num+1, vbead.gromacs_type_id, CGids_string)) + ret_lines.append( + "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string excl_lines.append(vsite_exclusions) - ret_lines.extend(excl_lines) + ret_lines.extend(excl_lines) - ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) - ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) + ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) + - else: - mol_count = 0 - for res in self._molecules[mol].resnames: - if res not in mapping: - logger.warning( - "Residue '{0}' present in bonding file, but not in mapping.".format(res) + not_calc) - continue - if res not in self._molecules: - logger.warning(("Residue '{0}' has no internal bonds").format(res)) - - ret_lines.append("\n[ moleculetype ]") - - ret_lines.append("{0:4s} {1:4d}".format(mol, 1)) - - ret_lines.append("\n[ atoms ]") - for resid in range(1, len(self._molecules[mol].resnames) + 1): - beads = self._molecules[mol].resid_to_beads[resid] - resname = self._molecules[mol].resid_to_resname[resid] - for bead in beads: - # atnum type resnum resname atname c-group charge (mass) - if isinstance(bead, VirtualMap): - ret_lines.append(atom_template["mass"].format( - bead.num + 1, bead.type, resid, resname, bead.name, bead.num + 1, bead.charge, bead.mass - )) - else: - ret_lines.append(atom_template["nomass"].format( - bead.num + 1, bead.type, resid, resname, bead.name, bead.num + 1, bead.charge - )) - - virtual_beads = self.get_virtual_beads(self._molecules[mol].beads) - if len(virtual_beads) != 0: - ret_lines.append("\n[ virtual_sitesn ]") - excl_lines = ["\n[ exclusions ]"] # exlusions section for virtual sites - for resid in range(1, len(self._molecules[mol].resnames) + 1): - beads = self._molecules[mol].resid_to_beads[resid] - virtual_beads = self.get_virtual_beads(beads) - for vbead in virtual_beads: - CGids = [bead.num + 1 for bead in beads if bead.name in vbead.atoms] - CGids.sort() - CGids_string = " ".join(map(str, CGids)) - ret_lines.append( - "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) - vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string - excl_lines.append(vsite_exclusions) - ret_lines.extend(excl_lines) - - ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) - ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) - ret_lines.extend( - write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) - ret_lines.extend( - write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) - mol_count += 1 return ret_lines def apply(self, frame): From 46ab34ef65a6d7b6653fff10045b6266d62c0ce2 Mon Sep 17 00:00:00 2001 From: Jon Shearer Date: Sat, 23 Feb 2019 15:04:38 +0000 Subject: [PATCH 12/13] removed any overlap between Molecule instances and Bead instances --- pycgtool/bondset.py | 120 ++++++++++++++++++++++--------------------- pycgtool/util.py | 15 ------ test/test_bondset.py | 4 +- 3 files changed, 65 insertions(+), 74 deletions(-) diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 0479856..d789037 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -41,36 +41,30 @@ class Molecule: """ Holds data for a molecule comprised of multiple residues """ - __slots__ = ["resnames", "bonds", "beads", "resname_to_beads"] + __slots__ = ["resnames", "bonds"] - def __init__(self, resnames=None, bonds=None, beads=None): - ''' - - :param resnames: list of resname names - :param bonds: list of lists of BeadMap objects in the same order as 'resnames' - :param beads: list of bonds - ''' + def __init__(self, resnames=None, bonds=None): + """ + :param resnames: list of residue names + :param bonds: list of lists of Bond objects in the same order as 'resnames' + """ self.resnames = resnames - if beads and resnames is not None: - self.resname_to_beads = dict(zip(resnames, beads)) - self.beads = [bead for res_beads in beads for bead in res_beads] - self.bonds = bonds if bonds is None: self.bonds = [] def add_bond(self, bond): - ''' + """ Add bond to object :param bond: instance of class with 'Bond' as their base class - ''' + """ self.bonds.append(bond) @property def inter(self): - '''collects inter residue bonds''' + """collects inter residue bonds""" inter = [] for bond in self.bonds: if isinstance(bond, GlobalBond): @@ -79,7 +73,7 @@ def inter(self): @property def is_multiresidue(self): - '''checks if Molecule has multiple residues''' + """checks if Molecule has multiple residues""" return True if len(self.inter) > 0 else False def __len__(self): @@ -185,10 +179,15 @@ def __init__(self, atoms, atom_numbers=None, func_form=None): self.resids = resids self.resnames = resnames - def get_atoms(self, residues): + def get_atoms(self, frame): + """ + get atoms involved in bond + :param frame: Frame object + :return: list of Atom objects + """ atoms = [] for resid, resname, name in zip(self.resids, self.resnames, self.atoms): - for residue in residues: + for residue in frame: if residue.num == resid: if residue.name == resname: try: @@ -199,6 +198,11 @@ def get_atoms(self, residues): return atoms def get_residue_ids(self, frame): + """ + Get internal resids of residues involved with bond + :param frame: + :return: list of resids + """ residue_ids = [] for resid, resname, name in zip(self.resids, self.resnames, self.atoms): for ind, residue in enumerate(frame): @@ -208,6 +212,10 @@ def get_residue_ids(self, frame): return residue_ids def populate_ids(self, mol_beads): + """ + populate internal indices of bond + :param mol_beads: beads in the molecule + """ ids = [] for resid, resname, name in zip(self.resids, self.resnames, self.atoms): for ind, beads in enumerate(mol_beads, start=1): @@ -332,8 +340,8 @@ def __init__(self, filename, options): if is_global: if options.generate_angles or options.generate_dihedrals: logger.warning("Automated generation of angles or dihedrals between " - "residues not implemented! Please specifiy angles and dihedrals in [< {0} >] " - "section of *.bnd file".format(mol_name)) + "residues not implemented! Please specifiy angles and dihedrals in [< {0} >]" + "section of *.bnd file".format(mol_name)) else: angles, dihedrals = self._create_angles(mol_bonds) @@ -356,7 +364,6 @@ def __init__(self, filename, options): ) ) - @staticmethod def _create_angles(mol_bonds): """ @@ -483,23 +490,25 @@ def connect_residues(self, frame, mapping): connects residues together to form new molecules :param frame: Frame from which to determine inter residue connections """ - #TODO this section is a bit verbose - simplify + # TODO this section is a bit verbose - simplify not_multi = [mol for mol in self._molecules if not self._molecules[mol].is_multiresidue] for mol in self._molecules: if self._molecules[mol].is_multiresidue: mol_bonds = self._molecules[mol].inter - fragments = [] + fragments_resid = [] for bond in mol_bonds: - fragments.append(bond.get_residue_ids(frame)) + fragments_resid.append(bond.get_residue_ids(frame)) self._populate_atom_numbers(mapping) - molecule_internal_ids = merge_list_of_lists(fragments) - if len(molecule_internal_ids) > 1: - print("All Fragments of Molecule '{}' are not connected - add missing bonds to .bnd".format(mol)) + molecule_internal_resids = merge_list_of_lists(fragments_resid) + if len(molecule_internal_resids) > 1: + print("All fragments of Molecule '{}' are not connected - add missing bonds to .bnd".format(mol)) raise SyntaxError - molecule_internal_ids = molecule_internal_ids[0] - resnames = [frame[mol_id].name for mol_id in molecule_internal_ids] + + # populate residue bond ids + molecule_internal_resids = molecule_internal_resids[0] + resnames = [frame[resid].name for resid in molecule_internal_resids] all_bonds = [] all_beads = [] start = 0 @@ -509,6 +518,7 @@ def connect_residues(self, frame, mapping): if len(not_multi) > 0: bonds = list(map(deepcopy, self._molecules[resname])) all_bonds.extend(bonds) + for i, bead in enumerate(beads): bead.num = start + i @@ -526,22 +536,12 @@ def connect_residues(self, frame, mapping): all_beads.append(beads) start = beads[-1].num + 1 - + # populate inter-residue bonds ids for bond in mol_bonds: bond.populate_ids(all_beads) - for mol in self._molecules: - num = 0 - for res in frame: - if res.name == mol: - num += 1 - if num > 1: - logger.warning("More than one {0} residue found in gro, note that pycgtool " - "does not currently support the fitting of bonding parameters of multi-residue " - "molecules + single residue molecules simultaneously.") - all_bonds.extend(mol_bonds) - molecule = Molecule(resnames, all_bonds, all_beads) + molecule = Molecule(resnames, all_bonds) self._molecules[mol] = molecule return @@ -557,7 +557,7 @@ def write_itp(self, filename, mapping): def itp_text(self, mapping): atom_template = {"nomass": "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f}", - "mass": "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f} {7:8.3f}"} + "mass": "{0:4d} {1:4s} {2:4d} {3:4s} {4:4s} {5:4d} {6:8.3f} {7:8.3f}"} def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity=None, rad2deg=False): ret_lines = [] @@ -598,8 +598,10 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= ret_lines.append("\n[ atoms ]") + # print residues + start = 1 for resid, res in enumerate(molecule.resnames, start=1): - beads = molecule.resname_to_beads[res] if molecule.is_multiresidue else mapping[mol] + beads = mapping[res] if res not in mapping: @@ -610,38 +612,40 @@ def write_bond_angle_dih(bonds, section_header, print_fconst=True, multiplicity= # atnum type resnum resname atname c-group charge (mass) if isinstance(bead, VirtualMap): ret_lines.append(atom_template["mass"].format( - bead.num + 1, bead.type, resid, res, bead.name, bead.num + 1, bead.charge, bead.mass + start + bead.num, bead.type, resid, res, bead.name, start + bead.num, bead.charge, bead.mass )) else: ret_lines.append(atom_template["nomass"].format( - bead.num + 1, bead.type, resid, res, bead.name, bead.num + 1, bead.charge + start + bead.num, bead.type, resid, res, bead.name, start + bead.num, bead.charge )) + start = beads[-1].num + 2 - all_beads = molecule.beads if molecule.is_multiresidue else mapping[mol] - virtual_beads = self.get_virtual_beads(all_beads) + # print virtual sites + virtual_beads = [self.get_virtual_beads(mapping[res]) for res in self._molecules[mol].resnames] + virtual_beads = [bead for res_beads in virtual_beads for bead in res_beads] if len(virtual_beads) != 0: ret_lines.append("\n[ virtual_sitesn ]") - excl_lines = ["\n[ exclusions ]"] #exlusions section for virtual sites + excl_lines = ["\n[ exclusions ]"] for res in molecule.resnames: - beads = molecule.resname_to_beads[res] if molecule.is_multiresidue else mapping[mol] + beads = mapping[res] virtual_beads = self.get_virtual_beads(beads) for vbead in virtual_beads: - CGids = [bead.num + 1 for bead in beads if bead.name in vbead.atoms] - CGids.sort() - CGids_string = " ".join(map(str, CGids)) + cg_ids = [bead.num + 1 for bead in beads if bead.name in vbead.atoms] + cg_ids.sort() + cg_ids_string = " ".join(map(str, cg_ids)) ret_lines.append( - "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, CGids_string)) - vsite_exclusions = "{} ".format(vbead.num + 1) + CGids_string + "{0:^6d} {1:^6d} {2}".format(vbead.num + 1, vbead.gromacs_type_id, cg_ids_string)) + vsite_exclusions = "{} ".format(vbead.num + 1) + cg_ids_string excl_lines.append(vsite_exclusions) ret_lines.extend(excl_lines) - ret_lines.extend(write_bond_angle_dih(self.get_bond_lengths(mol), "bonds")) ret_lines.extend(write_bond_angle_dih(self.get_bond_angles(mol), "angles", rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, rad2deg=True)) - ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", print_fconst=False)) - + ret_lines.extend(write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", multiplicity=1, + rad2deg=True)) + ret_lines.extend(write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", + print_fconst=False)) return ret_lines diff --git a/pycgtool/util.py b/pycgtool/util.py index 6230043..19f5c01 100644 --- a/pycgtool/util.py +++ b/pycgtool/util.py @@ -67,21 +67,6 @@ def vector_cross(u, v): res[2] = u[0] * v[1] - u[1] * v[0] return res -def shift_func_value(func_to_decorate, shift, periodic=None): - """ - shift value that a function returns - :param func_to_decorate: - :param shift: value to shift return value by - :return: new_func: new function - """ - - def new_func(*original_args, **original_kwargs): - value = func_to_decorate(*original_args, **original_kwargs) - shift - if periodic is not None: - value -= np.rint(value / two_pi) * two_pi - return value - return new_func - def circular_mean(values): """ return average of values on a cirle diff --git a/test/test_bondset.py b/test/test_bondset.py index ea10f4f..73253c4 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -313,8 +313,8 @@ def test_connect_residues(self): measure[mol].bonds[1].eqm = 90. measure[mol].bonds[1].fconst = 100. measure.connect_residues(frame, mapping) - self.assertEqual(measure[mol].beads[-1].num, 3) self.assertListEqual(measure[mol].bonds[0].atom_numbers, [1, 2]) + self.assertListEqual(measure[mol].bonds[1].atom_numbers, [1, 2, 3]) def test_global_itp(self): mol = "mol_01" @@ -327,6 +327,8 @@ def test_global_itp(self): measure[mol].bonds[1].eqm = np.deg2rad(30.) measure[mol].bonds[1].fconst = 100. measure.connect_residues(frame, mapping) + + logging.disable(logging.WARNING) measure.write_itp("global.itp", mapping) logging.disable(logging.NOTSET) From dfadc0ab535f266aae3e86c3196534cb1ac62c02 Mon Sep 17 00:00:00 2001 From: jonathan Date: Mon, 25 Feb 2019 12:26:12 +0000 Subject: [PATCH 13/13] updated documentation --- doc/source/index.rst | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/doc/source/index.rst b/doc/source/index.rst index e100b5e..492ed14 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -193,13 +193,12 @@ generate_dihedrals Generate dihedrals from bonds **False**, Tru Connect Residues ~~~~~~~~~~~~~~~~ -In the .bnd file the [ global ] section can be used to form bonds, angles or dihedrals between any residue in a system. -The syntax for a bond in the global section is:: +In the .bnd file the [< mol >] section can be used to form bonds, angles or dihedrals between any residue in a system. +The syntax for an inter residue bond is:: [atom_name1]_[resid1]_[resname1] [atom_name2]_[resid2]_[resname2] or angle:: [atom_name1]_[resid1]_[resname1] [atom_name2]_[resid2]_[resname2] [atom_name3]_[resid3]_[resname3] -Currently this feature only works for a single molecule and the name of the molecule outputted will always be 'mol_0'. -Also itp reading is not currently supported for multi-residue itps +Currently this feature only works for a single molecule. Note that itp reading is not currently supported for multi-residue itps Indexes =======