← 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/protein_folding_examples/1bead+chaperone/frustrated/moltemplate_files/generate_tables/calc_dihedral_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/protein_folding_examples/1bead+chaperone/frustrated/moltemplate_files/generate_tables/calc_dihedral_table.py
1
 
#!/usr/bin/env python
2
 
 
3
 
# Calculate a table of dihedral angle interactions used in the alpha-helix
4
 
# and beta-sheet regions of the frustrated protein model described in
5
 
#  provided in figure 8 of the supplemental materials section of:
6
 
# AI Jewett, A Baumketner and J-E Shea, PNAS, 101 (36), 13192-13197, (2004)
7
 
# Note that the "A" and "B" parameters were incorrectly reported to be 
8
 
# 5.4*epsilon and 6.0*epsilon.  The values used were 5.6 and 6.0 epsilon.
9
 
# The phiA and phiB values were 57.29577951308232 degrees (1 rad) 
10
 
# and 180 degrees, respectively.  Both expA and expB were 6.0.
11
 
#
12
 
# To generate the table used for the alpha-helix (1 degree resolution) use this:
13
 
#  ./calc_dihedral_table.py 6.0 57.29577951308232 6  5.6 180 6  0.0 359  360
14
 
# To generate the table used for the beta-sheets (1 degree resolution) use this:
15
 
#  ./calc_dihedral_table.py 5.6 57.29577951308232 6  6.0 180 6  0.0 359  360
16
 
#
17
 
# (If you're curious as to why I set the location of the minima at phi_alpha
18
 
#  to 1.0 radians (57.2957795 degrees), there was no particularly good reason.
19
 
#  I think the correct value turns out to be something closer to 50 degrees.)
20
 
 
21
 
 
22
 
from math import *
23
 
import sys
24
 
 
25
 
 
26
 
# The previous version included the repulsive core term
27
 
def U(phi, A, phiA, expA, B, phiB, expB, use_radians=False):
28
 
    conv_units = pi/180.0
29
 
    if use_radians:
30
 
        conv_units = 1.0
31
 
    termA = pow(cos(0.5*(phi-phiA)*conv_units), expA)
32
 
    termB = pow(cos(0.5*(phi-phiB)*conv_units), expB)
33
 
    return -A*termA - B*termB
34
 
 
35
 
# The previous version included the repulsive core term
36
 
def F(phi, A, phiA, expA, B, phiB, expB, use_radians=False):
37
 
    conv_units = pi/180.0
38
 
    if use_radians:
39
 
        conv_units = 1.0
40
 
    termA = (0.5*sin(0.5*(phi-phiA)*conv_units) * 
41
 
             expA * pow(cos(0.5*(phi-phiA)*conv_units), expA-1.0))
42
 
    termB = (0.5*sin(0.5*(phi-phiB)*conv_units) * 
43
 
             expB * pow(cos(0.5*(phi-phiB)*conv_units), expB-1.0))
44
 
    return -conv_units*(A*termA + B*termB)
45
 
 
46
 
if len(sys.argv) != 10:
47
 
    sys.stderr.write("Error: expected 9 arguments:\n"
48
 
                     "\n"
49
 
                     "Usage: "+sys.argv[0]+" A phiA expA B phiB expB phiMin phiMax N\n\n")
50
 
    sys.exit(-1)
51
 
 
52
 
A       = float(sys.argv[1])
53
 
phiA    = float(sys.argv[2])
54
 
expA    = float(sys.argv[3])
55
 
B       = float(sys.argv[4])
56
 
phiB    = float(sys.argv[5])
57
 
expB    = float(sys.argv[6])
58
 
phi_min = float(sys.argv[7])
59
 
phi_max = float(sys.argv[8])
60
 
N       =   int(sys.argv[9])
61
 
 
62
 
for i in range(0,N):
63
 
    phi = phi_min + i*(phi_max - phi_min)/(N-1)
64
 
    U_phi = U(phi, A, phiA, expA, B, phiB, expB, use_radians=False)
65
 
    F_phi = F(phi, A, phiA, expA, B, phiB, expB, use_radians=False)
66
 
    print(str(i+1)+' '+str(phi)+' '+str(U_phi)+' '+str(F_phi))
67