-
Notifications
You must be signed in to change notification settings - Fork 0
/
frcmod2omm.py
158 lines (125 loc) · 4.26 KB
/
frcmod2omm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
'''
Script contains conversion factors and functions for FRCMOD to OMM
conversion. Following parameters can be converted with this script.
1. HarmonicBondForce
2. HarmonicAngleForce
3. PeriodicTorsionForce: Proper
4. PeriodicTorsionForce: Improper
5. NonbondedForce
@author: Madhuranga Rathnayake
'''
import math
radian2degree = 57.2957795130 # 1 rad = 57.2957795130 degree
kcal2kj = 4.184 # 1 kcal = 4.184 kj
ang2nm = 0.1
def _bond(line):
"""(list:str) -> str
Parameters: list: processed line from .frcmod file
Return : Force Constant (K) and min.dist (r)
K: kcal/mol/(A**2) -> K/2: 2*kj/mol/nm**2
r: Ang -> nm
----
c3-h1 330.60 1.097
"""
atoms = line[:5].replace('-', ' ').split()
params = line[5:].split()
k = float(params[0])
r = float(params[1])
omm_t1 = atoms[0]
omm_t2 = atoms[1]
omm_k = 2*k*kcal2kj/(ang2nm*ang2nm)
omm_r = r*ang2nm
omm_out = '<Bond type1="{}" type2="{}" length="{}" k="{}"/>'.format(omm_t1, omm_t2, omm_r, omm_k)
print(omm_out)
return(omm_out)
def _angle(line):
"""(list:str) -> str
Parameters: list: processed line from .frcmod file
Return : Force Constant (K) and min.angle (a)
K: kcal/mol/(rad**2) -> K/2: 2*kj/mol/(rad**2)
a: degrees -> rad
----
c3-n -c3 63.030 115.640
"""
atoms = line[:8].replace('-', ' ').split()
params = line[8:].split()
k = float(params[0])
a = float(params[1])
omm_t1 = atoms[0]
omm_t2 = atoms[1]
omm_t3 = atoms[2]
omm_k = 2*k*kcal2kj
omm_a = math.radians(a)
omm_out = '<Angle type1="{}" type2="{}" type3="{}" angle="{}" k="{}"/>'.format(omm_t1, omm_t2, omm_t3, omm_a, omm_k)
print(omm_out)
return(omm_out)
def _dihedral(line):
"""(list:str) -> str
Parameters: list: processed line from .frcmod file
Return : Force Constant (K) and min.angle (a)
K(PK/IDIVF): kcal -> K: kj/mol
phase: degrees -> rad
----
h1-c3-n -c3 6 0.000 0.000 2.000
"""
atoms = line[:11].replace('-', ' ').split()
params = line[11:].split()
idivf = float(params[0])
pk = float(params[1])
phase = float(params[2])
pn = float(params[3])
omm_t1 = atoms[0]
omm_t2 = atoms[1]
omm_t3 = atoms[2]
omm_t4 = atoms[3]
omm_k = (pk/idivf) * kcal2kj
omm_phase = math.radians(phase)
omm_pn = int(pn)
omm_out = '<Proper type1="{}" type2="{}" type3="{}" type4="{}" periodicity1="{}" phase1="{}" k1="{}"/>'.format(omm_t1, omm_t2, omm_t3, omm_t4, omm_pn, omm_phase, omm_k)
print(omm_out)
return(omm_out)
def _improper(line):
"""(list:str) -> str
Parameters: list: processed line from .frcmod file
Return : Force Constant (K) and min.angle (a)
K: kcal -> K: kj/mol
phase: degrees -> rad
----
c -c3-n -c3 1.1 180.0 2.0
^
out of plane
"""
atoms = line[:11].replace('-', ' ').split()
params = line[11:].split()
k = float(params[0])
phase = float(params[1])
pn = float(params[2])
omm_t1 = atoms[2] # this is the one in out-of-plane. must go 1st in OMM
omm_t2 = atoms[0]
omm_t3 = atoms[1]
omm_t4 = atoms[3]
omm_k = k*kcal2kj
omm_p = math.radians(phase)
omm_pn = int(pn)
omm_out = '<Improper type1="{}" type2="{}" type3="{}" type4="{}" periodicity1="{}" phase1="{}" k1="{}"/>'.format(omm_t1, omm_t2, omm_t3, omm_t4, omm_pn, omm_p, omm_k)
print(omm_out)
return(omm_out)
def _nonbonding(line):
"""(list:str) -> str
Parameters: list: processed line from .frcmod file
Return : Force Constant (K) and min.angle (a)
rmin/2: Ang -> nm
epsilon: kcal/mol -> kj/mol
----
c3 1.9080 0.1094
"""
llist = line.split()
sqrt62 = math.pow(2, 1/6)
atom_t1 = llist[0]
rmin2 = float(llist[1])
epsilon = float(llist[2])
omm_sigma = ang2nm * 2 * rmin2 / sqrt62
omm_epsilon = kcal2kj * epsilon
omm_out = '<Atom type="{}" charge="XXXX" sigma="{}" epsilon="{}"/>'.format(atom_t1, omm_sigma, omm_epsilon)
print(omm_out)
return(omm_out)