-
Notifications
You must be signed in to change notification settings - Fork 27
/
CMD.py
201 lines (173 loc) · 9.46 KB
/
CMD.py
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
##############################
## 2 # COMMAND LINE PARSING ## -> @CMD <-
##############################
import sys, logging, os, inspect
import DOC
def str2atom(a):
""" Helper function to parse atom strings given on the command line:
resid, resname/resid, chain/resname/resid, resname/resid/atom,
chain/resname/resid/atom, chain//resid, chain/resname/atom """
a = a.split("/")
if len(a) == 1: # Only a residue number:
return (None, None, int(a[0]), None)
if len(a) == 2: # Residue name and number (CYS/123):
return (None, a[0], int(a[1]), None)
if len(a) == 3:
if a[2].isdigit(): # Chain, residue name, residue number
return (None, a[1], int(a[2]), a[0])
else: # Residue name, residue number, atom name
return (a[2], a[0], int(a[1]), None)
return (a[3], a[1], int(a[2]), a[0])
def option_parser(args, options, lists, version=0):
# Check whether there is a request for help
if '-h' in args or '--help' in args:
DOC.help()
# Convert the option list to a dictionary, discarding all comments
options = dict([i for i in options if not type(i) == str])
# This information we would like to print to some files,
# so let's put it in our information class
options['Version'] = version
options['Arguments'] = args[:]
while args:
ar = args.pop(0)
options[ar].setvalue([args.pop(0) for i in range(options[ar].num)])
## LOGGING ##
# Set the log level and communicate which options are set and what is happening
# If 'Verbose' is set, change the logger level
logLevel = options["-v"] and logging.DEBUG or logging.INFO
logging.basicConfig(format='%(levelname)-7s %(message)s', level=logLevel)
logging.info('MARTINIZE, script version %s'%version)
logging.info('If you use this script please cite:')
logging.info('de Jong et al., J. Chem. Theory Comput., 2013, DOI:10.1021/ct300646g')
# To make the program flexible, the forcefield parameters are defined
# for multiple forcefield. We first check for a FF file in the current directory.
# Next we check for the FF in globals (for catenated scripts).
# Next we check in at the location of the script and the subdiretory FF.
try:
options['ForceField'] = globals()[options['-ff'].value.lower()]()
except KeyError:
try:
_tmp = __import__(options['-ff'].value.lower()+"_ff")
options['ForceField'] = getattr(_tmp, options['-ff'].value.lower())()
except ImportError:
try:
# We add the directory where the script resides and a possible "ForceFields" directory to the search path
# realpath() will make your script run, even if you symlink it :)
cmd_folder = os.path.realpath(os.path.dirname(inspect.getfile(inspect.currentframe())))
if cmd_folder not in sys.path:
sys.path.insert(0, cmd_folder)
# use this if you want to include modules from a subfolder
cmd_subfolder = os.path.realpath(os.path.dirname(inspect.getfile(inspect.currentframe()))) + "/ForceFields"
if cmd_subfolder not in sys.path:
sys.path.insert(0, cmd_subfolder)
_tmp = __import__(options['-ff'].value.lower()+"_ff")
options['ForceField'] = getattr(_tmp, options['-ff'].value.lower())()
except:
logging.error("Forcefield '%s' can not be loaded." % (options['-ff']))
sys.exit()
# if os.path.exists(options['-ff'].value.lower()+'_ff.py'):
# _tmp = __import__(options['-ff'].value.lower()+"_ff")
# options['ForceField'] = getattr(_tmp, options['-ff'].value.lower())()
# elif os.path.exists('ForceFields/'+options['-ff'].value.lower()+'_ff.py'):
# _tmp = __import__("ForceFields."+options['-ff'].value.lower()+'_ff',fromlist="ForceFields")
# options['ForceField'] = getattr(_tmp, options['-ff'].value.lower())()
# elif options['-ff'].value.lower() in globals():
# options['ForceField'] = globals()[options['-ff'].value.lower()]()
# else:
# logging.error("Forcefield '%s' can not be found."%(options['-ff']))
# sys.exit()
#except:
# logging.error("Forcefield '%s' can not be loaded."%(options['-ff']))
# sys.exit()
# Process the raw options from the command line
# Boolean options are set to more intuitive variables
options['Collagen'] = options['-collagen']
options['chHIS'] = options['-his']
options['ChargesAtBreaks'] = options['-cb']
options['NeutralTermini'] = options['-nt']
options['ExtendedDihedrals'] = options['-ed']
options['RetainHETATM'] = False # options['-hetatm']
options['SeparateTop'] = options['-sep']
options['MixedChains'] = False # options['-mixed']
options['ElasticNetwork'] = options['-elastic']
# Parsing of some other options into variables
options['ElasticMaximumForce'] = options['-ef'].value
options['ElasticMinimumForce'] = options['-em'].value
options['ElasticLowerBound'] = options['-el'].value
options['ElasticUpperBound'] = options['-eu'].value
options['ElasticDecayFactor'] = options['-ea'].value
options['ElasticDecayPower'] = options['-ep'].value
options['ElasticBeads'] = options['-eb'].value.split(',')
options['PosResForce'] = options['-pf'].value
options['PosRes'] = [i.lower() for i in options['-p'].value.split(",")]
if "none" in options['PosRes']: options['PosRes'] = []
if "backbone" in options['PosRes']: options['PosRes'].append("BB")
if options['ForceField'].ElasticNetwork:
# Some forcefields, like elnedyn, always use an elatic network.
# This is set in the forcefield file, with the parameter ElasticNetwork.
options['ElasticNetwork'] = True
# Merges, links and cystines
options['mergeList'] = "all" in lists['merges'] and ["all"] or [i.split(",") for i in lists['merges']]
# Process links
linkList = []
linkListCG = []
for i in lists['links']:
ln = i.split(",")
a, b = str2atom(ln[0]), str2atom(ln[1])
if len(ln) > 3: # Bond with given length and force constant
bl, fc = (ln[2] and float(ln[2]) or None, float(ln[3]))
elif len(a) == 3: # Constraint at given distance
bl, fc = float(a[2]), None
else: # Constraint at distance in structure
bl, fc = None, None
# Store the link, but do not list the atom name in the
# atomistic link list. Otherwise it will not get noticed
# as a valid link when checking for merging chains
linkList.append(((None, a[1], a[2], a[3]), (None, b[1], b[2], b[3])))
linkListCG.append((a, b, bl, fc))
# Cystines
# This should be done for all special bonds listed in the _special_ dictionary
CystineCheckBonds = False # By default, do not detect cystine bridges
CystineMaxDist2 = (10*0.22)**2 # Maximum distance (A) for detection of SS bonds
for i in lists['cystines']:
if i.lower() == "auto":
CystineCheckBonds = True
elif i.replace(".", "").isdigit():
CystineCheckBonds = True
CystineMaxDist2 = (10*float(i))**2
else:
# This item should be a pair of cysteines
cysA, cysB = [str2atom(j) for j in i.split(",")]
# Internally we handle the residue number shifted by ord(' ')<<20.
# We have to add this to the cys-residue numbers given here as well.
constant = 32 << 20
linkList.append((("SG", "CYS", cysA[2]+constant, cysA[3]),
("SG", "CYS", cysB[2]+constant, cysB[3])))
linkListCG.append((("SC1", "CYS", cysA[2]+constant, cysA[3]),
("SC1", "CYS", cysB[2]+constant, cysB[3]), -1, -1))
# Now we have done everything to it, we can add Link/cystine related stuff to options
# 'multi' is not stored anywhere else, so that we also add
options['linkList'] = linkList
options['linkListCG'] = linkListCG
options['CystineCheckBonds'] = CystineCheckBonds
options['CystineMaxDist2'] = CystineMaxDist2
options['multi'] = lists['multi']
logging.info("Chain termini will%s be charged"%(options['NeutralTermini'] and " not" or ""))
logging.info("Residues at chain brakes will%s be charged"%((not options['ChargesAtBreaks']) and " not" or ""))
if 'ForceField' in options:
logging.info("The %s forcefield will be used."%(options['ForceField'].name))
else:
logging.error("Forcefield '%s' has not been implemented."%(options['-ff']))
sys.exit()
if options['ExtendedDihedrals']:
logging.info('Dihedrals will be used for extended regions. (Elastic bonds may be more stable)')
else:
logging.info('Local elastic bonds will be used for extended regions.')
if options['PosRes']:
logging.info("Position restraints will be generated.")
logging.warning("Position restraints are only enabled if -DPOSRES is set in the MDP file")
if options['MixedChains']:
logging.warning("So far no parameters for mixed chains are available. This might crash the program!")
if options['RetainHETATM']:
logging.warning("I don't know how to handle HETATMs. This will probably crash the program.")
return options