Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76392865
calc_table.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Aug 7, 17:29
Size
2 KB
Mime Type
text/x-python
Expires
Fri, Aug 9, 17:29 (2 d)
Engine
blob
Format
Raw Data
Handle
19658117
Attached To
rLAMMPS lammps
calc_table.py
View Options
#!/usr/bin/env python
# Calculate a table of pairwise energies and forces between "INT" atoms
# in the lipid membrane model described in
# Brannigan et al, Phys Rev E, 72, 011915 (2005)
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
# I realized later this is not what we want because although energy is conserved
# all enrgies are shifted with respect to energies used in the Brannigan paper
# (by 0.27 kCal/mole) and the later Watson JCP 2011 paper (by 0.224 kCal/mole).
# (So don't use this.)
# Calculate and print a
def
S
(
r
,
rc1
,
rc2
,
derivative
=
False
):
"""
Calculate the switching function S(r) which decays continuously
between 1 and 0 in the range from rc1 to rc2 (rc2>rc1):
S(r) = (rc2^2 - r^2)^2 * (rc2^2 + 2*r^2 - 3*rc1^2) / (rc2^2-rc1^2)^3
I'm using the same smoothing/switching cutoff function used by the CHARMM
force-fields. (I'm even using the same code to implement it, taken
from lammps charmm/coul/charmm pair style, rewritten in python.)
"""
assert
(
rc2
>
rc1
)
rsq
=
r
*
r
rc1sq
=
rc1
*
rc1
rc2sq
=
rc2
*
rc2
denom_lj_inv
=
(
1.0
/
((
rc2sq
-
rc1sq
)
*
(
rc2sq
-
rc1sq
)
*
(
rc2sq
-
rc1sq
)))
if
rsq
>
rc2sq
:
return
0.0
elif
rsq
<
rc1sq
:
if
derivative
:
return
0.0
else
:
return
1.0
else
:
rc2sq_minus_rsq
=
(
rc2sq
-
rsq
)
rc2sq_minus_rsq_sq
=
rc2sq_minus_rsq
*
rc2sq_minus_rsq
if
derivative
:
return
(
12.0
*
rsq
*
rc2sq_minus_rsq
*
(
rsq
-
rc1sq
)
*
denom_lj_inv
)
else
:
return
(
rc2sq_minus_rsq_sq
*
(
rc2sq
+
2.0
*
rsq
-
3.0
*
rc1sq
)
*
denom_lj_inv
)
def
U
(
r
,
eps
,
sigma
):
return
eps
*
(
0.4
*
pow
((
sigma
/
r
),
12
)
-
3.0
*
sigma
*
sigma
/
(
r
*
r
))
def
F
(
r
,
eps
,
sigma
):
return
eps
*
(
12
*
0.4
*
pow
((
sigma
/
r
),
13
)
/
sigma
-
2
*
3.0
*
sigma
*
sigma
/
(
r
*
r
*
r
))
epsilon
=
2.75
/
4.184
# kCal/mole
sigma
=
7.5
Rmin
=
0.02
Rmax
=
22.6
Rc1
=
22.0
Rc2
=
22.5
N
=
1130
for
i
in
range
(
0
,
N
):
r
=
Rmin
+
i
*
(
Rmax
-
Rmin
)
/
(
N
-
1
)
U_r
=
U
(
r
,
epsilon
,
sigma
)
F_r
=
F
(
r
,
epsilon
,
sigma
)
# Multiply U(r) & F(r) by the smoothing/switch function
U_r
=
U_r
*
S
(
r
,
Rc1
,
Rc2
)
F_r
=
U_r
*
S
(
r
,
Rc1
,
Rc2
,
True
)
+
F_r
*
S
(
r
,
Rc1
,
Rc2
,
False
)
print
(
str
(
i
+
1
)
+
' '
+
str
(
r
)
+
' '
+
str
(
U_r
)
+
' '
+
str
(
F_r
))
Event Timeline
Log In to Comment