Page MenuHomec4science
Files F62173412

calc_table.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 11, 09:52

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

c4science · Help