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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
|
# Copyright 2007-2008 Nanorex, Inc. See LICENSE file for details.
"""
YukawaPotential.py
Create an .xvg file suitable for passing as the argument to -table for
mdrun.
@author: Eric M
@version: $Id$
@copyright: 2007-2008 Nanorex, Inc. See LICENSE file for details.
The table consists of 7 columns:
r f f'' g g'' h h''
In other words, three functions evaluated at evenly spaced r values,
and their second derivatives. For the standard Coulomb and
Lennard-Jones potentials, the functions are:
f(r) = 1/r
g(r) = -1/r^6
h(r) = 1/r^12
This table is read by the GROMACS mdrun process when user specified
potential functions are used. We use this for the non-bonded
interactions between PAM5 model DNA helices.
"""
# DELTA_R should be 0.002 for a gromacs compiled to use single
# precision floats, and 0.0005 for doubles.
DELTA_R = 0.0005
# values from Physics, Halliday and Resnick, Third Edition, 1978
ELECTRON_CHARGE = 1.6021892e-19 # Coulombs (A s)
VACUUM_PERMITTIVITY = 8.854187818e-12 # Farads/meter (A^2 s^4 kg^-1 m^-3)
BOLTZMANN = 1.380662e-23 # Joules/Kelvin (kg m^2 s^-2 K^-1)
AVOGADRO = 6.022045e23 # particles/mol
PHOSPHATE_DIAMETER = 0.4 # nm
import math
class YukawaPotential(object):
def __init__(self, simulatorParameters):
self.rCutoff = simulatorParameters.getYukawaRCutoff()
self.rSwitch = simulatorParameters.getYukawaRSwitch()
self.shift = simulatorParameters.getYukawaShift()
self.charge = simulatorParameters.getYukawaCounterionCharge()
self.molarity = simulatorParameters.getYukawaCounterionMolarity()
self.temperature = simulatorParameters.getYukawaTemperatureKelvin()
self.dielectric = simulatorParameters.getYukawaDielectric()
self.fudge = simulatorParameters.getYukawaConstantMultiple()
# Bjerrum length (nm)
L_B = 1e9 * ELECTRON_CHARGE * ELECTRON_CHARGE / (4.0 * math.pi *
VACUUM_PERMITTIVITY *
self.dielectric *
BOLTZMANN *
self.temperature)
# (A^2 s^2) (A^-2 s^-4 kg m^3) (kg^-1 m^-2 s^2 K) K^-1 = m
# about 0.7nm for default conditions
# Debye length (nm)
Lambda_D = 1.0 / math.sqrt(4.0 * math.pi * L_B * self.molarity *
self.charge * self.charge)
# 1/sqrt(nm / litre) --> nm
# about 1.2nm for default conditions
mol_K_per_kJ = 1000 / (AVOGADRO * BOLTZMANN)
t1 = 1.0 + PHOSPHATE_DIAMETER / (2.0 * Lambda_D)
E_p_p = self.temperature * L_B / \
(mol_K_per_kJ * PHOSPHATE_DIAMETER * t1 * t1)
self.YA = E_p_p * PHOSPHATE_DIAMETER * self.fudge
self.YB = PHOSPHATE_DIAMETER
self.YC = Lambda_D
#print "L_B = %e, Lambda_D = %e" % (L_B, Lambda_D)
#print "t1 = %e, E_p_p = %e" % (t1, E_p_p)
#print "mol_K_per_kJ = %e" % mol_K_per_kJ
#print "YA = %e, YB = %e. YC = %e" % (self.YA, self.YB, self.YC)
self.YShift = 0.0
if (self.shift):
self.YShift = self.yukawa(self.rCutoff)
# Switching function to smoothly reduce a function to zero.
# Transition begins when r == start, below which the value is 1. The
# value smoothly changes to 0 by the time r == end, after which it
# remains 0.
#
# Potential functions are multiplied by the switching function S:
#
# (r/start - 1)^4
# S = (1 - ----------------- ) ^4
# (end/start - 1)^4
#
# in the range start < r < end.
self.switch_len = self.rCutoff - self.rSwitch
self.switch_len_4 = self.switch_len * self.switch_len * self.switch_len * self.switch_len
self.switch_d_1 = ((self.rCutoff / self.rSwitch) - 1)
self.switch_d_2 = -16.0 / (self.switch_d_1 * self.switch_d_1 * self.switch_d_1 * self.switch_d_1 * self.rSwitch)
self.switch_d2_1 = (self.rCutoff / self.rSwitch) - 1
self.switch_d2_1_4 = self.switch_d2_1 * self.switch_d2_1 * self.switch_d2_1 * self.switch_d2_1
self.switch_d2_1_8 = self.switch_d2_1_4 * self.switch_d2_1_4
self.switch_d2_1_8_start_2 = self.switch_d2_1_8 * self.rSwitch * self.rSwitch
if (self.switch_len <= 0.0):
self.func = self.yukawa
self.d2_func = self.d2_yukawa
else:
self.func = self.switch_yukawa
self.d2_func = self.d2_switch_yukawa
def writeToFile(self, pathName):
tableFile = open(pathName, 'w')
r = 0.0
# We go to 2 * self.rCutoff because GROMACS reuses the mdp
# option table-extension (how far past self.rCutoff we need to
# extend the table) as the length of the 1-4 interaction
# table. Since we want 1-4 interactions to go to self.rCutoff
# as well, we need table-extension to be self.rCutoff, which
# results in the normal table being 2 * self.rCutoff. Silly,
# really.
while (r < 2 * self.rCutoff + (DELTA_R / 2.0)):
print >>tableFile, \
"%8.4f %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e" % (r,
self.r_1(r),
self.d2_r_1(r),
self.func(r),
self.d2_func(r),
0,
0)
r += DELTA_R
tableFile.close()
def r_1(self, r):
if (r < 0.04):
return 0.0
return 1.0 / r
def d2_r_1(self, r):
if (r < 0.04):
return 0.0
return 2.0 / (r * r * r)
def r_6(self, r):
if (r < 0.04):
return 0.0
return -1.0 / (r * r * r * r * r * r)
def d2_r_6(self, r):
if (r < 0.04):
return 0.0
return -42.0 / (r * r * r * r * r * r * r * r)
def r_12(self, r):
if (r < 0.04):
return 0.0
return 1.0 / (r * r * r * r * r * r * r * r * r * r * r * r)
def d2_r_12(self, r):
if (r < 0.04):
return 0.0
return 156.0 / (r * r * r * r * r * r * r * r * r * r * r * r * r * r)
def yukawa(self, r):
if (r < 0.04):
return self.yukawa(0.04) + 0.04 - r ;
return (self.YA / r) * math.exp(-(r - self.YB) / self.YC) - self.YShift
def d_yukawa(self, r):
if (r < 0.04):
return 0.0
return self.YA * math.exp(-(r - self.YB) / self.YC) * (-1.0/(r * r) - 1.0/(self.YC * r))
def d2_yukawa(self, r):
if (r < 0.04):
return 0.0
return self.YA * math.exp(-(r - self.YB) / self.YC) * (2.0/(self.YC * r * r) + 1.0/(self.YC * self.YC * r) + 2.0/(r * r * r))
def switch(self, r):
if (r <= self.rSwitch):
return 1.0
if (r >= self.rCutoff):
return 0.0
rDiff = r - self.rSwitch
S1 = ((rDiff * rDiff * rDiff * rDiff) / self.switch_len_4) - 1
return S1 * S1 * S1 * S1
def d_switch(self, r):
if (r <= self.rSwitch):
return 0.0
if (r >= self.rCutoff):
return 0.0
t1 = r - self.rSwitch
t1_4 = t1 * t1 * t1 * t1
t2 = 1 - t1_4 / self.switch_len_4
t3 = (r / self.rSwitch) - 1
t4 = t2 * t2 * t2 * t3 * t3 * t3
return self.switch_d_2 * t4
def d2_switch(self, r):
if (r <= self.rSwitch):
return 0.0
if (r >= self.rCutoff):
return 0.0
t1 = r - self.rSwitch
t1_4 = t1 * t1 * t1 * t1
t2 = t1_4 / self.switch_len_4
t3 = (r / self.rSwitch) - 1
t3_2 = t3 * t3
t3_4 = t3_2 * t3_2
return (48 * (-(1 - t2) * self.switch_d2_1_4 + 4 * t3_4) * (t2 - 1) * (t2 - 1) * t3_2) / self.switch_d2_1_8_start_2
def switch_yukawa(self, r):
return self.switch(r) * self.yukawa(r)
def d2_switch_yukawa(self, r):
return self.d2_switch(r) * self.yukawa(r) + 2.0 * self.d_switch(r) * self.d_yukawa(r) + self.switch(r) * self.d2_yukawa(r)
|
