← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to tools/moltemplate/examples/CG_biomolecules/vesicle_Brannigan2005+Bellesia2010/moltemplate_files/calc_table/version_charmm_cutoff/calc_table.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
tools/moltemplate/examples/CG_biomolecules/vesicle_Brannigan2005+Bellesia2010/moltemplate_files/calc_table/version_charmm_cutoff/calc_table.py
1
 
#!/usr/bin/env python
2
 
 
3
 
# Calculate a table of pairwise energies and forces between "INT" atoms
4
 
# in the lipid membrane model described in 
5
 
#   Brannigan et al, Phys Rev E, 72, 011915 (2005)
6
 
# The energy of this interaction U(r) = eps*(0.4*(sigma/r)^12 - 3.0*(sigma/r)^2)
7
 
# I realized later this is not what we want because although energy is conserved
8
 
# all enrgies are shifted with respect to energies used in the Brannigan paper
9
 
# (by 0.27 kCal/mole) and the later Watson JCP 2011 paper (by 0.224 kCal/mole).
10
 
# (So don't use this.)
11
 
 
12
 
# Calculate and print a
13
 
 
14
 
def S(r, rc1, rc2, derivative=False):
15
 
    """ 
16
 
    Calculate the switching function S(r) which decays continuously
17
 
    between 1 and 0 in the range from rc1 to rc2 (rc2>rc1):
18
 
       S(r) = (rc2^2 - r^2)^2 * (rc2^2 + 2*r^2 - 3*rc1^2) / (rc2^2-rc1^2)^3
19
 
    I'm using the same smoothing/switching cutoff function used by the CHARMM
20
 
    force-fields.  (I'm even using the same code to implement it, taken 
21
 
    from lammps charmm/coul/charmm pair style, rewritten in python.)
22
 
 
23
 
    """
24
 
    assert(rc2>rc1)
25
 
    rsq           = r*r
26
 
    rc1sq         = rc1*rc1
27
 
    rc2sq         = rc2*rc2
28
 
    denom_lj_inv  = (1.0 / ((rc2sq-rc1sq)*
29
 
                            (rc2sq-rc1sq)*
30
 
                            (rc2sq-rc1sq)))
31
 
    if rsq > rc2sq:
32
 
        return 0.0
33
 
    elif rsq < rc1sq:
34
 
        if derivative:
35
 
            return 0.0
36
 
        else:
37
 
            return 1.0
38
 
    else:
39
 
        rc2sq_minus_rsq    = (rc2sq - rsq)
40
 
        rc2sq_minus_rsq_sq = rc2sq_minus_rsq * rc2sq_minus_rsq
41
 
        if derivative:
42
 
            return (12.0 * rsq * rc2sq_minus_rsq * (rsq-rc1sq) * denom_lj_inv)
43
 
        else:
44
 
            return (rc2sq_minus_rsq_sq *
45
 
                    (rc2sq + 2.0*rsq - 3.0*rc1sq) * denom_lj_inv)
46
 
 
47
 
 
48
 
def U(r, eps, sigma):
49
 
    return eps*   (0.4*pow((sigma/r),12) -   3.0*sigma*sigma/(r*r))
50
 
 
51
 
def F(r, eps, sigma):
52
 
    return eps*(12*0.4*pow((sigma/r),13)/sigma - 2*3.0*sigma*sigma/(r*r*r))
53
 
 
54
 
epsilon = 2.75/4.184 # kCal/mole
55
 
sigma   = 7.5
56
 
Rmin    = 0.02
57
 
Rmax    = 22.6
58
 
Rc1     = 22.0
59
 
Rc2     = 22.5
60
 
N       = 1130
61
 
 
62
 
for i in range(0,N):
63
 
    r = Rmin + i*(Rmax-Rmin)/(N-1)
64
 
    U_r = U(r, epsilon, sigma)
65
 
    F_r = F(r, epsilon, sigma)
66
 
    # Multiply U(r) & F(r) by the smoothing/switch function
67
 
    U_r = U_r * S(r, Rc1, Rc2)
68
 
    F_r = U_r * S(r, Rc1, Rc2, True)  +  F_r * S(r, Rc1, Rc2, False)
69
 
    print(str(i+1)+' '+str(r)+' '+str(U_r)+' '+str(F_r))
70