Thanks John & Charles,
Yes it is the problem of "print=verbose=loopinfo" inside the loop
main. I double checked and confirmed.
Anyway attached please see the script revised according to the
suggestion of Charles, i.e. moving the ivm setup block outside of the
main loop:
Please comment:
1 remarks Solvent refinement protocol from ARIA1.2 (Nilges and
Linge), modified for XPLOR-NIH
2 remarks RDC refinment protocol from XPLOR-NIH sa_tmv_bice_rgyr.inp
3 {* Water-refinement from Spronk's scripts *}
4 {* C.A.E.M. Spronk, S.B. Nabuurs, A.M.J.J. Bonvin, E. Krieger,
G.W. Vuister & G. Vriend. *}
5 {* The precision of NMR structure ensembles revisited.
*}
6 {* Journal of Biomolecular NMR, 25 (2003), 225-234.
*}
7 {* RDC from XPLOR-NIH/eginput/gb1_rdc/sa_tmv_bice_rgyr.inp
*}
8 {* W/ suggestions from CDS, 20060425 xplor-nih maillist
*}
9 set message on echo on end
10
{*==========================================================================*}
11 {*====================== SET FILENAMES AND VARIABLES
=======================*}
12
{*==========================================================================*}
13 {* type of non-bonded parameters *}
14 {* choice: "PROLSQ" "PARMALLH6" "PARALLHDG" "OPLSX" *}
15 {* The water refinement uses the OPLSX parameters *}
16 evaluate ( $par_nonbonded = "OPLSX" )
17 evaluate ( $maxResid = 122 )
18 {* parameter file(s) *}
19 evaluate ( $par.1 = "TOPPAR:parallhdg5.3.pro" )
20 evaluate ( $par.2 = "TOPPAR:parallhdg5.3.sol" )
21 evaluate ( $par.3 = "TOPPAR:axes.par" )
22 evaluate ( $par.4 = "" )
23 evaluate ( $par.5 = "" )
24 {* topology file *}
25 evaluate ( $protein_topology = "TOPPAR:topallhdg5.3.pro" )
26 evaluate ( $solvent_topology = "TOPPAR:topallhdg5.3.sol" )
27 evaluate ( $axis_topology = "TOPPAR:top_axis.pro" )
28 {* structure file(s) *}
29 evaluate ( $axis_struct = "TOPPAR:axis_500.psf" )
30 {* struct.1 is updated in nmv_xplor.py script *}
31 evaluate ( $struct.1 = "STRUCTURES:input.psf" )
32 evaluate ( $struct.2 = "" )
33 evaluate ( $struct.3 = "" )
34 evaluate ( $struct.4 = "" )
35 evaluate ( $struct.5 = "" )
36 {* input coordinate file(s) *}
37 {* pdb.in.file.1 is updated in nmv_xplor.py script *}
38 evaluate ( $pdb.in.file.1 = "INPUTCOORDINATES:analyzed_1.pdb" )
39 evaluate ( $pdb.in.file.2 = "axis_500.pdb" )
40 {* output coordinate file(s) *}
41 {* pdb.out.file.1 is updated in nmv_xplor.py script *}
42 evaluate ( $pdb.out.file.1 = "OUTPUTCOORDINATES:refined_1.pdb" )
43 {* input distance restraint file(s) *}
44 {* noe.file.1 is updated in nmv_xplor.py script *}
45 evaluate ( $noe.file.1 = "TABLES:mytest_noe.tbl" )
46 {* Averaging for NOE restraints *}
47 {* noe.ave is updated in nmv_xplor.py script *}
48 evaluate ( $noe.ave = sum )
49 {* input dihedral restraint file(s) *}
50 {* cdih.file.1 is updated in nmv_xplor.py script *}
51 evaluate ( $cdih.file.1 = "TABLES:mytest_dihedrals.tbl" )
52 {* input hbonding restraint file *}
53 {* hbda.file.1 is updated in nmv_xplor.py script *}
54 evaluate ( $hbda.file.1 = "TABLES:mytest_hdba.tbl" )
55 {* seed for velocity generation *}
56 evaluate ( $seed = 12396 )
57 {* mdsteps.heat|hot|cool are updated in nmv_xplor.py script *}
58 {* set number of md steps for the heating stage *}
59 evaluate ( $mdsteps.heat = 10 )
60 {* set number of md steps for the hot md stage *}
61 evaluate ( $mdsteps.hot = 40 )
62 {* set number of md steps for the cooling stage *}
63 evaluate ( $mdsteps.cool = 10 )
64 {* all acceptance criteria are updated in nmv_xplor.py script *}
65 {* set acceptance criteria *}
66 evaluate ( $accept.noe = 0.50 )
67 evaluate ( $accept.cdih = 5.00 )
68 evaluate ( $accept.coup = 1.00 )
69 evaluate ( $accept.sani = 0.00 )
70 evaluate ( $accept.vean = 5.00 )
71
{*==========================================================================*}
72 {*=========== READ THE STRUCTURE, TOPOLOGY AND PARAMETER FILES
=============*}
73
{*==========================================================================*}
74 !structure @$struct.1 @$axis_struct end
75 !topology @$solvent_topology end
76 topology @$solvent_topology @$protein_topology @$axis_topology end
77 parameter @$par.1 @$par.2 @$par.3 end
78 segment
79 name=" "
80 chain
81 first nter tail + * end
82 last cter head - * end
83 coordinates @$pdb.in.file.1
84 end
85 end
86 segment
87 name="AXIS"
88 chain
89 coordinates @$pdb.in.file.2
90 end
91 end
92 !write struc output=test.psf end
93
{*==========================================================================*}
94 {*================== SET VALUES FOR NONBONDED PARAMETERS
===================*}
95
{*==========================================================================*}
96 parameter
97 nbonds
98 nbxmod=5 atom cdiel shift
99 cutnb=9.5 ctofnb=8.5 ctonnb=6.5 eps=1.0 e14fac=0.4 inhibit 0.25
100 wmin=0.5
101 tolerance 0.5
102 end
103 end
104
{*==========================================================================*}
105 {*============== READ THE COORDINATES AND COPY TO REFERENCE SET
============*}
106
{*==========================================================================*}
107 coor @$pdb.in.file.1
108 coor @$pdb.in.file.2
109 !write coordinates from=main output=test.pdb end
110 vector do (refx = x) (all)
111 vector do (refy = y) (all)
112 vector do (refz = z) (all)
113
{*==========================================================================*}
114 {*========================= GENERATE THE WATER LAYER
=======================*}
115
{*==========================================================================*}
116 vector do (segid = "PROT") (segid " ")
117 @SOLVENT:generate_water.inp
118 vector do (segid = " ") (segid "PROT")
119
{*==========================================================================*}
120 {*========================= READ THE EXPERIMENTAL DATA
=====================*}
121
{*==========================================================================*}
122 noe reset
123 nrestraints = 10000
124 ceiling = 100
125 @$noe.file.1
126 averaging * $noe.ave
127 potential * soft
128 scale * 50
129 sqconstant * 1.0
130 sqexponent * 2
131 soexponent * 1
132 rswitch * 1.0
133 sqoffset * 0.0
134 asymptote * 2.0
135 end
136 restraints dihedral reset
137 @$cdih.file.1
138 scale = 200
139 end
140 {* collapse added *}
141 collapse
142 assign (resid 7:114) 100.0 12.03
143 scale 1.0
144 print
145 end
146 !distance/angle hbond term
147 hbda
148 nres 800
149 class back
150 @$hbda.file.1
151 force 200.0
152 end
153 rama
154 nres=10000
155 evaluate ($krama = 1.0)
156 @QUARTS:2D_quarts_new.tbl
157 @QUARTS:3D_quarts_new.tbl
158 @QUARTS:forces_torsion_prot_quarts_intra.tbl
159 end
160 @QUARTS:setup_quarts_torsions_intra_2D3D.tbl
161 ! did not tweak the coup
162 couplings
163 nres 400
164 potential harmonic
165 class phi
166 degen 1
167 force 1.0
168 coefficients 6.98 -1.38 1.72 -60.0
169 @TABLES:hnha3.tbl
170 end
171 evaluate ($ksani = 1)
172 evaluate ($ksani_CH = 1.0*$ksani)
173 evaluate ($ksani_CACO = 0.035*$ksani)
174 evaluate ($ksani_NCO = 0.050*$ksani)
175 evaluate ($ksani_HNC = 0.108*$ksani)
176 sani
177 nres=400
178 class t_nh
179 force $ksani
180 potential harmonic
181 coeff 0.0 -16.5 0.65
182 @TABLES:page_nh.tbl
183 end
184 sani class t_nh force $ksani end
185
{*==========================================================================*}
186 {*============================ SET FLAGS
===================================*}
187
{*==========================================================================*}
188 {* flags added hbda coll rama 20060419 *}
189 flags exclude *
190 include bond angle dihe impr vdw elec
191 noe cdih coup oneb carb ncs dani
192 vean sani prot harm hbda coll rama
193 end
194
{*==========================================================================*}
195 {*========================= START THE REFINEMENT
===========================*}
196
{*==========================================================================*}
197 dynamics internal
198 reset
199 !print=verbose=loopinfo
200 !print=verbose=loopdebug
201 ! print verbose=nodedef
202 !
203 ! dipolar axis
204 !
205 group (resid 500 )
206 hinge rotate (resid 500)
207 evaluate ($res = 1)
208 while ($res le $maxResid) loop group
209 {* group aromatic ring atoms. *}
210 group (resid $res and resname PHE and
211 (name CG or name CD1 or name CD2 or name CE1 or
name CE2 or name CZ))
212 group (resid $res and resname HIS and
213 (name CG or name ND1 or name CD2 or name CE1 or name NE2))
214 group (resid $res and resname TYR and
215 (name CG or name CD1 or name CD2 or name CE1 or
name CE2 or name CZ))
216 group (resid $res and resname TRP and
217 (name CG or name CD1 or name CD2 or name NE1 or
218 name CE2 or name CE3 or name CZ2 or name CZ3 or name CH2))
219 group (resid $res and resname PRO and
220 (name N or name CA or name CB or name CG or name CD))
221 group (resid $res and resname ASP and
222 (name CG or name OD1 or name OD2))
223 group (resid $res and resname ASN and
224 (name CG or name OD1 or name NE2))
225 group (resid $res and resname GLU and
226 (name CD or name OE1 or name OD2))
227 group (resid $res and resname GLN and
228 (name CD and name OE1 and name NE2))
229 group (resid $res and resname ARG and
230 (name NE or name CZ or name NH1 or name NH2))
231 evaluate ($res = $res + 1)
232 end loop group
233 set message on echo on end
234 cloop=false
235 auto torsion
236 maxe 10000
237 end
238 set seed $seed end
239 ! We loop until we have an accepted structure, maximum trials=3
240 evaluate ($end_count = 3)
241 evaluate ($count = 0)
242 while ($count < $end_count ) loop main
243 evaluate ($accept = 0)
244 ! since we do not use SHAKe, increase the water bond angle
energy constant
245 parameter
246 angle (resn tip3) (resn tip3) (resn tip3) 500 TOKEN
247 end
248 ! reduce improper and angle force constant for some atoms
249 evaluate ($kbonds = 1000)
250 evaluate ($kangle = 50)
251 evaluate ($kimpro = 5)
252 evaluate ($kchira = 5)
253 evaluate ($komega = 5)
254 parameter
255 angle (not resn tip3)(not resn tip3)(not resn tip3)
$kangle TOKEN
256 improper (all)(all)(all)(all) $kimpro TOKEN TOKEN
257 bond (not resn tip3 and not name H*) (not resn tip3 and
not name H*) $kbonds TOKEN
258 end
259 ! fix the protein for initial minimization
260 constraints fix (not resn tip3) end
261 dynamics internal
262 !print=verbose=loopinfo
263 itype=powell
264 nstep=200
265 maxcalls=200
266 nprint=10
267 etol=1e-6
268 gtol=0.1
269 depred=0.01
270 end
271 ! release protein and restrain harmonically
272 constraints fix (not all) end
273 vector do (refx=x) (all)
274 vector do (refy=y) (all)
275 vector do (refz=z) (all)
276 restraints harmonic
277 exponent = 2
278 end
279 vector do (harmonic = 0) (all)
280 vector do (harmonic = 10) (not name h*)
281 vector do (harmonic = 20.0)(resname ANI and name OO)
282 vector do (harmonic = 0.0) (resname ANI and name Z )
283 vector do (harmonic = 0.0) (resname ANI and name X )
284 vector do (harmonic = 0.0) (resname ANI and name Y )
285 constraints
286 interaction (not resname ANI) (not resname ANI)
287 interaction ( resname ANI) ( resname ANI)
288 end
289 dynamics internal
290 itype=powell
291 nstep=200
292 maxcalls=200
293 nprint=10
294 etol=1e-6
295 gtol=0.1
296 depred=0.01
297 end
298 vector do (refx=x) (not resname ANI)
299 vector do (refy=y) (not resname ANI)
300 vector do (refz=z) (not resname ANI)
301 dynamics internal
302 itype=powell
303 nstep=200
304 maxcalls=200
305 nprint=10
306 etol=1e-6
307 gtol=0.1
308 depred=0.01
309 end
310 vector do (refx=x) (not resname ANI)
311 vector do (refy=y) (not resname ANI)
312 vector do (refz=z) (not resname ANI)
313 vector do (mass =50) (all)
314 vector do (mass=1000) (resname ani)
315 vector do (fbeta = 0) (all)
316 vector do (fbeta = 20. {1/ps} ) (not resn ani)
317 evaluate ($kharm = 50)
318 ! heat to 500 K
319 for $bath in (100 200 300 400 500) loop heat
320 vector do (harm = $kharm) (not name h* and not resname ANI)
321 vector do (vx=maxwell($bath)) (all)
322 vector do (vy=maxwell($bath)) (all)
323 vector do (vz=maxwell($bath)) (all)
324 dynamics internal
325 itype=pc6
326 nstep=$mdsteps.heat
327 timestep=0.002 !{ps}
328 tbath=$bath
329 response= 20
330 response= 5
331 nprint=50
332 end
333 evaluate ($kharm = max(0, $kharm - 4))
334 vector do (refx=x) (not resname ANI)
335 vector do (refy=y) (not resname ANI)
336 vector do (refz=z) (not resname ANI)
337 end loop heat
338 ! refinement at high T:
339 constraints
340 interaction (not resname ANI) (not resname ANI) weights * 1
dihed 3 end
341 interaction ( resname ANI) ( resname ANI) weights * 1 end
342 end
343 vector do (harm = 0) (not resname ANI)
344 vector do (vx=maxwell($bath)) (all)
345 vector do (vy=maxwell($bath)) (all)
346 vector do (vz=maxwell($bath)) (all)
347 dynamics internal
348 itype=pc6
349 nstep=$mdsteps.hot
350 timest=0.002{ps}
351 tbath=$bath
352 response= 20
353 response= 5
354 nprint=50
355 end
356 constraints
357 interaction (not resname ANI) (not resname ANI) weights *
1 dihed 5 end
358 interaction ( resname ANI) ( resname ANI) weights * 1 end
359 end
360 ! cool
361 evaluate ($bath = 500)
362 while ($bath >= 25) loop cool
363 evaluate ($kbonds = max(1000,$kbonds / 1.1))
364 evaluate ($kangle = min(1000,$kangle * 1.1))
365 evaluate ($kimpro = min(300,$kimpro * 1.4))
366 evaluate ($kchira = min(1000,$kchira * 1.4))
367 evaluate ($komega = min(100,$komega * 1.4))
368 parameter
369 bond (not resn tip3 and not name H*)(not resn tip3
and not name H*) $kbonds TOKEN
370 angle (not resn tip3 and not name H*)(not resn tip3
and not name H*)(not resn tip3 and not name H*) $kangle TOKEN
371 improper (all)(all)(all)(all) $kimpro TOKEN TOKEN
372 !VAL: stereo CB
373 improper (name HB and resn VAL)(name CA and resn
VAL)(name CG1 and resn VAL)(name CG2 and resn VAL) $kchira TOKEN TOKEN
374 !THR: stereo CB
375 improper (name HB and resn THR)(name CA and resn
THR)(name OG1 and resn THR)(name CG2 and resn THR) $kchira TOKEN TOKEN
376 !LEU: stereo CG
377 improper (name HG and resn LEU)(name CB and resn
LEU)(name CD1 and resn LEU)(name CD2 and resn LEU) $kchira TOKEN TOKEN
378 !ILE: chirality CB
379 improper (name HB and resn ILE)(name CA and resn
ILE)(name CG2 and resn ILE)(name CG1 and resn ILE) $kchira TOKEN TOKEN
380 !chirality CA
381 improper (name HA)(name N)(name C)(name CB) $kchira TOKEN TOKEN
382 improper (name O) (name C) (name N) (name CA) $komega
TOKEN TOKEN
383 improper (name HN) (name N) (name C) (name CA) $komega
TOKEN TOKEN
384 improper (name CA) (name C) (name N) (name CA) $komega
TOKEN TOKEN
385 improper (name CD) (name N) (name C) (name CA) $komega
TOKEN TOKEN
386 end
387 vector do (vx=maxwell($bath)) (all)
388 vector do (vy=maxwell($bath)) (all)
389 vector do (vz=maxwell($bath)) (all)
390 dynamics internal
391 itype=pc6
392 nstep=$mdsteps.cool
393 timestep=0.002{ps}
394 tbath=$bath
395 nprint=50
396 end
397 evaluate ($bath = $bath - 25)
398 end loop cool
399 !final minimization:
400 dynamics internal
401 itype=powell
402 nstep=2000
403 maxcalls=2000
404 nprint=50
405 etol=1e-6
406 gtol=0.1
407 depred=0.01
408 end
409 !mini powell nstep=500 nprint=50 end
410
{*==========================================================================*}
411 {*======================= CHECK RESTRAINT VIOLATIONS
=======================*}
412
{*==========================================================================*}
413 print threshold=$accept.noe noe
414 evaluate ($v_noe = $violations)
415 print threshold=0.5 noe
416 evaluate ($v_noe_0.5 = $violations)
417 print threshold=0.4 noe
418 evaluate ($v_noe_0.4 = $violations)
419 print threshold=0.3 noe
420 evaluate ($v_noe_0.3 = $violations)
421 print threshold=0.2 noe
422 evaluate ($v_noe_0.2 = $violations)
423 print threshold=0.1 noe
424 evaluate ($v_noe_0.1 = $violations)
425 evaluate ($rms_noe = $result)
426 print threshold=$accept.cdih cdih
427 evaluate ($rms_cdih=$result)
428 evaluate ($v_cdih = $violations)
429 coupl print thres = $accept.coup class phi end
430 evaluate ($rms_coup = $result)
431 evaluate ($v_coup = $violations)
432 sani print threshold = $accept.sani class t_nh end
433 evaluate ($rms_sani = $result)
434 evaluate ($v_sani = $violations)
435 vean print threshold = $accept.vean class vea1 end
436 evaluate( $rms_vean = $result)
437 evaluate( $v_vean = $violations)
438 print thres=0.05 bonds
439 evaluate ($rms_bonds=$result)
440 evaluate ($v_bonds = $violations)
441 print thres=5. angles
442 evaluate ($rms_angles=$result)
443 evaluate ($v_angles = $violations)
444 print thres=5. impropers
445 evaluate ($rms_impropers=$result)
446 evaluate ($v_impropers = $violations)
447 if ($v_noe > 0) then evaluate ( $accept=$accept + 1 ) end if
448 if ($v_cdih > 0) then evaluate ( $accept=$accept + 1 ) end if
449 if ($v_coup > 0) then evaluate ( $accept=$accept + 1 ) end if
450 if ($v_sani > 0) then evaluate ( $accept=$accept + 1 ) end if
451 if ($v_vean > 0) then evaluate ( $accept=$accept + 1 ) end if
452 if ($accept = 0 ) then
453 evaluate ( $label = "ACCEPTED" )
454 exit main
455 else
456 evaluate ( $label = "NOT ACCEPTED" )
457 evaluate ( $count = $count + 1 )
458 end if
459 end loop main
460
{*==========================================================================*}
461 {*================ COMPUTE RMS DIFFERENCE BETWEEN OBSERVED
=================*}
462 {*========================= AND MODEL DISTANCES.
===========================*}
463
{*==========================================================================*}
464 constraints interaction
465 (not resname TIP* and not resname ANI)
466 (not resname TIP* and not resname ANI)
467 end
468 energy end
469 remarks Structure $label
470 remarks E-overall: $ener
471 remarks E-NOE_restraints: $noe
472 remarks E-CDIH_restraints: $cdih
473 remarks E-COUP_restraints: $coup
474 remarks E-SANI_restraints: $sani
475 remarks E-VEAN_restraints: $vean
476 remarks RMS-NOE_restraints: $rms_noe
477 remarks RMS-CDIH_restraints: $rms_cdih
478 remarks RMS-COUP_restraints: $rms_coup
479 remarks RMS-SANI_restraints: $rms_sani
480 remarks RMS-VEAN_restraints: $rms_vean
481 remarks NOE Acceptance criterium: $accept.noe
482 remarks NOE violations > $accept.noe: $v_noe
483 remarks All NOE viol. (>0.5,0.4,0.3,0.2,0.1A): $v_noe_0.5
$v_noe_0.4 $v_noe_0.3 $v_noe_0.2 $v_noe_0.1
484 remarks CDIH Acceptance criterium: $accept.cdih
485 remarks CDIH violations > $accept.cdih: $v_cdih
486 remarks COUP Acceptance criterium: $accept.coup
487 remarks COUP violations > $accept.coup: $v_coup
488 remarks SANI Acceptance criterium: $accept.sani
489 remarks SANI violations > $accept.sani: $v_sani
490 remarks VEAN Acceptance criterium: $accept.vean
491 remarks VEAN violations > $accept.vean: $v_vean
492 write coordinates sele= (not resn TIP3) output =$pdb.out.file.1 end
493 stop
