Skip to content
Snippets Groups Projects
Commit b1facdc8 authored by James Graham's avatar James Graham
Browse files

Removed factor of 2 from Boltzmann Inversion 

Started adding r2b file to ffout
parent 8a8e2b34
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
clean.sh 0 → 100755
+ 2
− 0
View file @ b1facdc8
rm *.itp* *.gro* *.dat* *.json* *.xtc*
rm -r *.ff
pycgtool/bondset.py
+ 6
− 3
View file @ b1facdc8
......@@ -47,6 +47,9 @@ class Bond:
def __len__(self):
return len(self.atoms)
def __iter__(self):
return iter(self.atoms)
def boltzmann_invert(self, temp=310, default_fc=True):
"""
Perform Boltzmann Inversion using measured values of bond to calculate equilibrium value and force constant.
......@@ -58,9 +61,9 @@ class Bond:
rt = 8.314 * temp / 1000.
rad2 = np.pi * np.pi / (180. * 180.)
conv = {2: lambda: 1250. if default_fc else rt / (2 * var),
3: lambda: 25. if default_fc else rt / (2 * np.sin(np.radians(mean))**2 * var * rad2),
4: lambda: rt / (2 * var * rad2)}
conv = {2: lambda: 1250. if default_fc else rt / var,
3: lambda: 25. if default_fc else rt / (np.sin(np.radians(mean))**2 * var * rad2),
4: lambda: rt / (var * rad2)}
self.eqm = mean
try:
......
pycgtool/forcefield.py
+ 38
− 0
View file @ b1facdc8
......@@ -4,6 +4,7 @@ This module contains a single class ForceField used to output a GROMACS .ff forc
import os
import shutil
import functools
from .util import dir_up
from .parsers.cfg import CFG
......@@ -66,6 +67,27 @@ class ForceField:
line += " {0:4d}".format(multiplicity)
print(line, file=file)
def any_starts_with(iterable, char):
"""
Return True if any atoms of any bonds in molecule start with 'char'.
i.e. if char='-' or '+' is part of polymer.
:param iterable: Iterable of bond entries to check
:param char: Char to check each atom name for startswith, in '-+'
:return: True if any atom name in molecule bonds starts with char, else False
"""
recurse = functools.partial(any_starts_with, char=char)
if type(iterable) is str:
return iterable.startswith(char)
else:
return any(map(recurse, iterable))
def strip_polymer_bonds(bonds, char):
return [bond for bond in bonds if not any_starts_with(bond, char)]
# def write_residue(name):
with open(os.path.join(self.dirname, name), "w") as rtp:
print("[ bondedtypes ]", file=rtp)
print(("{:4d}" * 8).format(1, 1, 1, 1, 1, 1, 0, 0), file=rtp)
......@@ -84,9 +106,25 @@ class ForceField:
bead.name, bead.type, bead.charge, 0
), file=rtp)
needs_terminal_entry = [False, False]
bond_tmp = bonds.get_bond_lengths(mol, with_constr=True)
write_bond_angle_dih(bond_tmp, "bonds", rtp)
needs_terminal_entry[0] |= any_starts_with(bond_tmp, "-")
needs_terminal_entry[1] |= any_starts_with(bond_tmp, "+")
bond_tmp = bonds.get_bond_angles(mol)
write_bond_angle_dih(bond_tmp, "angles", rtp)
needs_terminal_entry[0] |= any_starts_with(bond_tmp, "-")
needs_terminal_entry[1] |= any_starts_with(bond_tmp, "+")
bond_tmp = bonds.get_bond_dihedrals(mol)
write_bond_angle_dih(bond_tmp, "dihedrals", rtp, multiplicity=1)
needs_terminal_entry[0] |= any_starts_with(bond_tmp, "-")
needs_terminal_entry[1] |= any_starts_with(bond_tmp, "+")
# print(mol, needs_terminal_entry)
# bond_tmp = bonds.get_bond_lengths(mol, with_constr=True)
# if needs_terminal_entry[0]:
# print(strip_polymer_bonds(bond_tmp, "-"))
test/data/sugar_out.itp
+ 17
− 17
View file @ b1facdc8
......@@ -17,27 +17,27 @@ ALLA 1
6 P4 1 ALLA O5 6 0.000
[ bonds ]
1 2 1 0.22209 88271.49992
2 3 1 0.21998 17581.05011
3 4 1 0.21648 32772.55156
4 5 1 0.26454 9511.08560
5 6 1 0.20452 26728.60329
2 3 1 0.21998 35162.10022
3 4 1 0.21648 65545.10313
4 5 1 0.26454 19022.17121
5 6 1 0.20452 53457.20659
[ angles ]
1 2 3 1 77.69557 210.24268
2 3 4 1 115.59599 198.52391
3 4 5 1 108.51038 431.97481
4 5 6 1 81.23827 575.61056
5 6 1 1 146.49898 259.05013
6 1 2 1 99.84045 833.86832
1 2 3 1 77.69557 420.48535
2 3 4 1 115.59599 397.04782
3 4 5 1 108.51038 863.94961
4 5 6 1 81.23827 1151.22112
5 6 1 1 146.49898 518.10027
6 1 2 1 99.84045 1667.73664
[ dihedrals ]
1 2 3 4 1 -85.45004 117.05203 1
2 3 4 5 1 65.21378 58.73066 1
3 4 5 6 1 -22.73796 47.10087 1
4 5 6 1 1 55.18087 18.96053 1
5 6 1 2 1 -96.52349 20.36670 1
6 1 2 3 1 73.82477 129.60063 1
1 2 3 4 1 -85.45004 234.10405 1
2 3 4 5 1 65.21378 117.46132 1
3 4 5 6 1 -22.73796 94.20175 1
4 5 6 1 1 55.18087 37.92107 1
5 6 1 2 1 -96.52349 40.73340 1
6 1 2 3 1 73.82477 259.20126 1
[ constraints ]
1 2 1 0.22209
6 1 1 0.17931
test/test_bondset.py
+ 18
− 18
View file @ b1facdc8
......@@ -48,24 +48,24 @@ class BondSetTest(unittest.TestCase):
self.assertEqual(0, len(angles))
def test_bondset_boltzmann_invert(self):
ref = [( 0.222817161647, 9560719),
( 0.216766804211, 27930081),
( 0.223815134299, 3099614),
( 0.239439104442, 322294637),
( 0.215497208685, 3683331),
( 0.179544192953, 10260492),
( 76.9617550012, 13048242),
(116.009786726, 927619),
(109.247878273, 195899),
( 84.4955293511, 871991),
(149.557001531, 471274),
( 99.0524708671, 100943),
(-89.5321545897, 91248),
( 70.040930903, 169002),
(-21.1194544981, 995309),
( 48.3727747848, 35712),
(-85.923133935, 62394),
( 70.3195444564, 778389)]
ref = [( 0.222817161647, 19121439),
( 0.216766804211, 55860162),
( 0.223815134299, 6199228),
( 0.239439104442, 644589275),
( 0.215497208685, 7366663),
( 0.179544192953, 20520985),
( 76.9617550012, 26096485),
(116.009786726, 1855239),
(109.247878273, 391799),
( 84.4955293511, 1743982),
(149.557001531, 942548),
( 99.0524708671, 201886),
(-89.5321545897, 182496),
( 70.040930903, 338005),
(-21.1194544981, 1990619),
( 48.3727747848, 71425),
(-85.923133935, 124789),
( 70.3195444564, 1556778)]
measure = BondSet("test/data/sugar.bnd", DummyOptions)
frame = Frame("test/data/sugar.gro", xtc="test/data/sugar.xtc")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment