! ! slow cooling protocol in torsion angle space for protein G. Uses only ! NOE and dihedral angle constraints. ! ! CDS 5/2/00 ! eval ($numStructs = 3) !total number of structures to calculate eval ($randomSeed = 785) !random seed eval ($use_graphics = false) ! VMD-xplor graphics are not used. eval ($maxResid = 56) !used for properly configuring TA dynamics ! ! get parallel info, which computer and how many processes. ! cpyth "from os import environ as env" cpyth "xplor.command('eval ($proc_num=%s)' % env['XPLOR_PROCESS'] )" cpyth "xplor.command('eval ($num_procs=%s)' % env['XPLOR_NUM_PROCESSES'])" eval ($num_procs=min($num_procs,$numStructs)) if ( $proc_num >= $num_procs ) then stop end if eval ($firstStruct = int($proc_num * $numStructs / $num_procs)) eval ($lastStruct = int(($proc_num+1) * $numStructs / $num_procs)) ! ! read in the PSF and initial PDB files ! parameter @TOPPAR:protein.par @TOPPAR:axes.par end structure @protG.psf @axis_600.psf end ! ! get parallel info ! !cpyth "from os import environ as env" !cpyth "xplor.command('eval ($proc_num=%s)' % env['XPLOR_PROCESS'] )" !cpyth "xplor.command('eval ($num_procs=%s)' % env['XPLOR_NUM_PROCESSES'])" noe {====>} nres=3000 {*Estimate greater than the actual number of NOEs.*} class all {====>} @protG_NOE.tbl {*Read NOE distance ranges.*} end couplings nres 400 potential harmonic class phi degen 1 force 1.0 coefficients 6.98 -1.38 1.72 -60.0 {====>} @jna_coup.tbl {* regular (non-Glycine) hnha couplings *} end restraints dihed reset nass = 5000 {====>} @protG_dih.tbl {*Read dihedral angle restraints.*} end ! ! starting coords ! coor @template_protG.pdb coor @axis_600.pdb ! Fix dipolar coupling tensor axis ! this applies for the powell minimization- but not for the torsion angle ! dynamics constraints fix (resname "ANI") end flags exclude * include bond angle impr end mini powell nstep 1000 nprint 100 end {* Minimize with only the covalent constraints. *} coor copy end ! ! annealing settings ! eval ($init_t = 3500.01) eval ($final_t = 100) eval ($cool_steps = 25000) eval ($tempstep = 25) eval ($ncycle_vdw = 100) eval ($ncycle = ($init_t-$final_t)/$tempstep) eval ($nstep = int($cool_steps/$ncycle)) eval ($ini_rad = 0.4) eval ($fin_rad = 0.80) eval ($radfact = ($fin_rad/$ini_rad)^(1/$ncycle)) eval ($ini_con = 0.004) eval ($fin_con = 4.0) eval ($k_vdwfact = ($fin_con/$ini_con)^(1/$ncycle)) eval ($ini_ang = 0.4) eval ($fin_ang = 1.0) eval ($ang_fac = ($fin_ang/$ini_ang)^(1/$ncycle)) eval ($ini_imp = 0.4) ! was 0.1 eval ($fin_imp = 1.0) eval ($imp_fac = ($fin_imp/$ini_imp)^(1/$ncycle)) eval ($hitemp_noe = 20.0) ! was 150 eval ($ini_noe = 20.0) ! was 2 eval ($fin_noe = 30.0) eval ($noe_fac = ($fin_noe/$ini_noe)^(1/$ncycle)) eval ($ini_sani = 0.01) eval ($fin_sani = 1.0) evaluate ($sani_fac = ($fin_sani/$ini_sani)^(1/$ncycle)) eval ($ini_timestep = 0.010) !reduced from 0.015 eval ($fin_timestep = 0.003) eval ($timestep_fac = ($fin_timestep/$ini_timestep)^(1/$ncycle)) !scaling factors for the vdw loop eval ($radfact_vdw = ($fin_rad/$ini_rad)^(1/$ncycle_vdw)) eval ($k_vdwfact_vdw = ($fin_con/$ini_con)^(1/$ncycle_vdw)) eval ($ang_fac_vdw = ($fin_ang/$ini_ang)^(1/$ncycle_vdw)) eval ($imp_fac_vdw = ($fin_imp/$ini_imp)^(1/$ncycle_vdw)) vector do (mass = 100.0) (all) vector do (fbeta = 10.0) (all) set echo on message on end ! ! set up the torsion angle dynamics topology stuff ! to eliminate degrees of freedom in the aromatics dynamics internal set echo off message off end evaluate ( $res = 1 ) while ($res le $maxResid) loop group {* group aromatic ring atoms. *} group (resid $res and resname PHE and (name CG or name CD1 or name CD2 or name CE1 or name CE2 or name CZ)) group (resid $res and resname PRO and (name N or name CA or name CB or name CG or name CD)) group (resid $res and resname HIS and (name CG or name ND1 or name CD2 or name CE1 or name NE2)) group (resid $res and resname TYR and (name CG or name CD1 or name CD2 or name CE1 or name CE2 or name CZ)) group (resid $res and resname TRP and (name CG or name CD1 or name CD2 or name NE1 or name CE2 or name CE3 or name CZ2 or name CZ3 or name CH2)) evaluate ( $res = $res + 1 ) end loop group !set up the dipolar coupling tensor axis group (resid 600) hinge rotate (resid 600) set echo on message on end cloop=false auto torsion maxe 10000 end constraints interaction (not resname ANI) (not resname ANI) weights * 1 angl $ini_ang impr $ini_imp {* Scale covalent constraints. *} end end parameter {*Parameters for the repulsive energy term.*} nbonds repel=0.5 {*Initial value for effective atom radius *} {*--modified later.*} rexp=2 irexp=2 rcon=1. nbxmod=-2 {*Initial value for nbxmod--modified later.*} wmin=0.01 cutnb=4.5 ctonnb=2.99 ctofnb=3. tolerance=0.5 end end restraints dihedral scale=5. {*Initial weight--modified later.*} end evaluate ($ksani = $ini_sani) sani nres=1600 class JNHp force $ksani potential harmonic {====>} coeff 0.0 -9.9 0.23 @bicelles_new_nh.tbl {* dipolar coupling restraints for protein amide NH. *} end eval ($count = $firstStruct) while ($count < $lastStruct) loop structure ! ! Be sure different structures start with different seeds ! eval ($seed = $randomSeed+$count) set seed $seed end if ( $use_graphics = true ) then ps defi $count bonds (not name h) end end end if ! ! reset the force constants and call the energy terms ! eval ($bath = $init_t) eval ($radius = $ini_rad) eval ($k_vdw = $ini_con) eval ($k_ang = $ini_ang) eval ($k_imp = $ini_imp) eval ($knoe = $ini_noe) eval ($ksani = $ini_sani) eval ($timestep = $ini_timestep) flags exclude * include bonds angl impr noe coup cdih sani end ! ! re-init the coords ! coor swap end coor copy end {* Save the first structure to copy to use as a reference *} ! ! set some high-temp force constants ! noe ceiling=1000 averaging * sum potential * soft sqconstant * 1.0 sqexponent * 2 sqoffset * 0.0 soexponent * 1 asymptote * 1.0 rswitch * 0.5 rswitch * 3 scale * $hitemp_noe end sani class JNHp force $ksani end parameters nbonds atom nbxmod 4 {* Can use 4 here, due to internal coordinate dynamics. *} wmin 0.01 cutnb 4.5 tolerance 0.5 rexp 2 irex 2 repel $radius rcon $k_vdw end end vector do (vx = maxwell(3500)) (all) vector do (vy = maxwell(3500)) (all) vector do (vz = maxwell(3500)) (all) {* Set initial velocities to fit a temperature of 3500K. *} {* High temperature to promote convergence. *} ! ! high temp dynamics ! evaluate ($tol = $bath/1000) dynamics internal nstep 0 endt 20 timestep $timestep tbath $bath response 20 nprint 100 etol $tol end if ( $use_graphics = true ) then ps append $count bonds (not name h) end end end if flags exclude * include bond angle impr vdw noe coup cdih sani end eval ($bath = $init_t) evaluate ($i_vdw = 0) while ($i_vdw < $ncycle_vdw) loop vdw {* iterate changes in van der *} {* Waals interaction. Scale and*} {* radius. *} evaluate ($i_vdw=$i_vdw+1) {* Min function is used to keep variables within ini_variable and *} {* fin_variable. *} eval ($k_vdw = min($fin_con,$k_vdw*$k_vdwfact_vdw)) eval ($radius = min($fin_rad,$radius*$radfact_vdw)) eval ($k_ang = min($fin_ang,$k_ang*$ang_fac_vdw)) eval ($k_imp = min ($fin_imp,$k_imp*$imp_fac_vdw)) parameter nbonds rcon $k_vdw repel $radius end end constraints interaction (all) (all) weights * 1 angl $k_ang impr $k_imp end end ! ! vdw dynamics {* dynamics with Van der Waals repulsion turned on. *} ! dynamics internal timestep $timestep endt 1 tbath $bath nprint $nstep ntrfr 10 end if ( $use_graphics = true ) then ps append $count bonds (not name h) end end end if end loop vdw ! ! cooling ! restraints dihedral scale=200. end {* increase dihedral term *} evaluate ($i_cool = 0) while ($i_cool < $ncycle) loop cool evaluate ($i_cool=$i_cool+1) eval ($bath = $bath - $tempstep) eval ($knoe = $knoe*$noe_fac) {* Scale during cooling. *} evaluate ($ksani = $ksani*$sani_fac) eval ($timestep = $timestep*$timestep_fac) eval ($k_vdw = min($fin_con,$k_vdw*$k_vdwfact)) eval ($radius = min($fin_rad,$radius*$radfact)) eval ($k_ang = min($fin_ang,$k_ang*$ang_fac)) eval ($k_imp = min($fin_imp,$k_imp*$imp_fac)) noe scale * $knoe end sani class JNHp force $ksani end parameter nbonds rcon $k_vdw repel $radius end end constraints interaction (not resname ANI) (not resname ANI) weights * 1 angl $k_ang impr $k_imp {* Scale impropers and angles *} end end ! ! cooling dynamics ! evaluate ($tol = $bath/1000) dynamics internal endt 2 tbath $bath nprint $nstep ntrfr 10 end if ( $use_graphics = true ) then ps append $count bonds (not name h) end end end if end loop cool ! ! final minimization - with fixed tensor axis ! constraints fix (resname "ANI") end mini powell nstep 500 nprint 50 end if ( $use_graphics = true ) then ps append $count bonds (not name h) end end end if ! ! recenter ! coor orient end ! ! analysis ! evaluates RMS deviations and violations print threshold 0.5 noe eval ($rms_noe = $result) eval ($violations_noe = $violations) sani print threshold=0.0 class JNHp end evaluate ($rms_sani_JNH_prot=$result) evaluate ($viol_sani_JNH_prot=$violations) evaluate ($R_JNH_prot=$result/12.77) print thres 0.05 bonds eval ($rms_bonds = $result) print thres 5. angles eval ($rms_angles = $result) print thres 5. impropers eval ($rms_impropers = $result) print thres 5. cdih eval ($rms_cdih = $result) eval ($violations_cdih=$violations) couplings print threshold 1.0 class phi end evaluate ($rms_coup = $result) evaluate ($viols_coup = $violations) remarks overall bonds angles improper vdw noe coup sani remarks energies: $ener $bond $angl $impr $vdw $noe $coup $sani remarks =============================================================== remarks bonds angles impropers noe Jcoupings sani remarks $rms_bonds $rms_angles $rms_impropers $rms_noe $rms_coup $rms_sani_JNH_prot remarks =============================================================== remarks noe dihedrals couplings remarks violations : $violations_noe $violations_cdih $viols_coup remarks =============================================================== eval ($suffix = ".pdb_nonconverged") {* Test for convergence. *} if ($violations_noe < 4) then if ($violations_cdih < 2) then eval ($suffix = ".pdb") {* Change the suffix if converged. *} end if end if evaluate ($file = "1gb1_anneal_" + encode($count) + $suffix) write coor output= $file end eval ($count = $count + 1) end loop structure stop