REMARK This protocol has REMARKS very slow cooling with increase of vdw evaluate ($seed=8294876) set seed $seed end !---------------------------------------------------------------------- ! read in the PSF file and initial structure param @TOPPAR:protein.par @TOPPAR:axes.par end structure @e1_his.psf @TOPPAR:axis_500.psf end coor @e1_lr_his_final_ave.min coor @axis_new_alternate.pdb vector do (mass=0.5) (resid 500) !---------------------------------------------------------------------- ! set the weights for the experimental energy terms evaluate ($knoe = 30.0) ! noes evaluate ($asym = 0.1) ! slope of NOE potential evaluate ($kcdi = 10.0) ! torsion angles evaluate ($krama = 0.002) !rama !---------------------------------------------------------------------- ! The next statement makes sure the experimental energies are used in the ! calculation, and switches off the unwanted energy terms. ! note that the NMR torsions are only switched on in the cooling stage ! we include the noncrystallographic symmetry right from the start !--------------------------------------------------------------------------- ! Read experimental restraints set echo on message on end noe reset nrestraints = 6000 ! allocate space for NOEs set echo off message off end class all @noe.tbl @hbonds.tbl @hb_his.tbl end set echo on message on end noe ceiling 1000 averaging all sum potential all square scale all $knoe sqconstant all 1.0 sqexponent all 2 end set echo on message on end couplings nres 300 potential harmonic class phi degen 1 force 1.0 set echo off message off end coefficients 6.98 -1.38 1.72 -60.0 @phij_nogly.tbl !regular hnha couplings class gly !for gly where you don't know stereoassignement degen 2 force 1.0 0.2 coefficients 6.98 -1.38 1.72 0.00 @phij_gly.tbl !gly hnha couplings set echo on message on end end carbon phistep=180 psistep=180 nres=300 class all force 0.5 potential harmonic @C13SHIFTS:rcoil_c13.tbl !rcoil shifts @C13SHIFTS:expected_edited_c13.tbl !13C shift database set echo off message off end @cacb_shifts.tbl !carbon shifts set echo on message on end end evaluate ($kdani = 0.01) dani nres=300 class relax force $kdani potential harmonic ! coef 13.0 2.0 0.4 600.13 60.82 !Tc Anis Rhombicity Wh Wn! coef 13.16 1.71 0.30 600.13 60.82 set echo off message off end @T1T2_rhomb.tbl set echo on message on end end restraints dihed scale $kcdi nass = 5000 set message off echo off end @dihed_180.tbl !dihedral angle restraints. @dihed_his.tbl set message on echo on end end vector do (refx=x) (all) vector do (refy=y) (all) vector do (refz=z) (all) flags exclude * include bonds angl impr vdw noe cdih coup carb dani end evaluate ($cool_steps = 3000) evaluate ($init_t = 3000.01) vector do (mass = 100.0) (all) ! uniform mass for all atoms vector do (fbeta = 10.0) (all) ! coupling to heat bath coor copy end evaluate ($rcon = 0.003) parameters nbonds atom nbxmod 4 wmin = 0.01 ! warning off cutnb = 4.5 ! nonbonded cutoff tolerance 0.5 repel= 0.9 ! scale factor for vdW radii = 1 ( L-J radii) rexp = 2 ! exponents in (r^irex - R0^irex)^rexp irex = 2 rcon=$rcon ! actually set the vdW weight end end {* Generate Structures 1 -> 50 *} evaluate ($count =0) while ($count < 1) loop structure evaluate ($count = $count + 1) {====>} {*Filename(s) for embedded coordinates.*} vector do (x=xcomp) (all) vector do (y=ycomp) (all) vector do (z=zcomp) (all) vector do (vx =maxwell($init_t)) (not resid ANI) vector do (vy =maxwell($init_t)) (not resid ANI) vector do (vz =maxwell($init_t)) (not resid ANI) evaluate ($ini_rad = 0.9) evaluate ($fin_rad = 0.80) evaluate ($ini_con= 0.004) evaluate ($fin_con= 4.0) evaluate ($ini_ang = 0.4) evaluate ($fin_ang = 1.0) evaluate ($ini_imp = 0.1) evaluate ($fin_imp = 1.0) evaluate ($ini_noe = 2.0) evaluate ($fin_noe = 30.0) evaluate ($knoe = $ini_noe) ! slope of NOE potential evaluate ($ini_rama = 0.002) evaluate ($fin_rama = 1.0) evaluate ($krama = $ini_rama) evaluate ($kcdi = 10.0) ! torsion angles evaluate ($ini_dani = 0.01) evaluate ($fin_dani = 1.0) evaluate ($kdani = $ini_dani) dani class relax force $kdani end noe averaging all sum potential all square scale all $knoe sqconstant all 1.0 sqexponent all 2 end restraints dihed scale $kcdi end evaluate ($rcon = 1.0) parameters nbonds atom nbxmod 4 wmin = 0.01 ! warning off cutnb = 100 ! nonbonded cutoff tolerance 45 repel= 1.2 ! scale factor for vdW radii = 1 ( L-J radii) ! Just for CA atoms rexp = 2 ! exponents in (r^irex - R0^irex)^rexp irex = 2 rcon=$rcon ! actually set the vdW weight end end constraints interaction (not name ca) (all) weights * 1 angl 0.4 impr 0.1 vdw 0 elec 0 end interaction (name ca) (name ca) weights * 1 angl 0.4 impr 0.1 vdw 1.0 end end vector do (store1 = decode(resid)) (name CA) vector show min (store1) (name CA) evaluate ($first_residue = $result) vector show max (store1) (name CA) evaluate ($last_residue = $result) {* This is just a clever way to count the amino acids *} {* but it doesn't need to be changed if you change the protein.*} dynamics internal reset group (resid 500 ) hinge rotate (resid 500) {* Treat the axis as a single unit. *} evaluate ($res = $first_residue) while ($res le $last_residue) 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 end {* Must exit the dynamics set-up in order to use vector show command to *} {* identify prolines so they can be broken up. *} evaluate ( $res = 1 ) while ($res le $last_residue) loop protein_proline vector show element (resname) (resid $res) eval ($variable =$result) if ($variable = PRO) then dynamics internal reuse=on {* Reuse to keep the assigned grouping of the aromatics. *} break (resid $res and name CD) (resid $res and name N) hinge bendtorsion (resid $res and name CA) (resid $res and name CB) (resid $res and name N) hinge bendtorsion (resid $res and name CB) (resid $res and name CG) (resid $res and name CA) hinge bendtorsion (resid $res and name CG) (resid $res and name CD) (resid $res and name CB) end end if evaluate ( $res = $res + 1 ) end loop protein_proline dynamics internal set message on echo on end cloop=false {* This is changed when the dynamics actually starts. *} maxe = 10000.0 auto torsion {* a ? prints out the parameters associated with the internal dynamics *} {* so that you can tell what you actually did in setting them. This *} {* is mostly useful in testing your scripts. *} end evaluate ($tol = $init_t/1000) dynamics internal {* High temperature dynamics with VDW between CA only.*} nstep 0 endt 20 timestep 0.001 tbath $init_t response 20 nprint 100 etol $tol end parameters nbonds atom nbxmod 4 wmin = 0.01 ! warning off cutnb = 4.5 ! nonbonded cutoff tolerance 0.5 repel= 0.9 ! scale factor for vdW radii = 1 ( L-J radii) rexp = 2 ! exponents in (r^irex - R0^irex)^rexp irex = 2 rcon =1.0 ! actually set the vdW weight end end evaluate ($kcdi = 200) restraints dihed scale $kcdi end evaluate ($final_t = 100) { K } evaluate ($tempstep = 25) { K } evaluate ($ncycle = ($init_t-$final_t)/$tempstep) evaluate ($nstep = int($cool_steps*4.0/$ncycle)) evaluate ($bath = $init_t) evaluate ($k_vdw = $ini_con) evaluate ($k_vdwfact = ($fin_con/$ini_con)^(1/$ncycle)) evaluate ($radius= $ini_rad) evaluate ($radfact = ($fin_rad/$ini_rad)^(1/$ncycle)) evaluate ($k_ang = $ini_ang) evaluate ($ang_fac = ($fin_ang/$ini_ang)^(1/$ncycle)) evaluate ($k_imp = $ini_imp) evaluate ($imp_fac = ($fin_imp/$ini_imp)^(1/$ncycle)) evaluate ($noe_fac = ($fin_noe/$ini_noe)^(1/$ncycle)) evaluate ($knoe = $ini_noe) evaluate ($dani_fac = ($fin_dani/$ini_dani)^(1/$ncycle)) evaluate ($kdani = $ini_dani) evaluate ($rama_fac = ($fin_rama/$ini_rama)^(1/$ncycle)) evaluate ($krama = $ini_rama) vector do (vx = maxwell($bath)) (all) vector do (vy = maxwell($bath)) (all) vector do (vz = maxwell($bath)) (all) evaluate ($i_cool = 0) while ($i_cool < $ncycle) loop cool evaluate ($i_cool=$i_cool+1) evaluate ($bath = $bath - $tempstep) evaluate ($k_vdw=min($fin_con,$k_vdw*$k_vdwfact)) evaluate ($radius=max($fin_rad,$radius*$radfact)) evaluate ($k_ang = $k_ang*$ang_fac) evaluate ($k_imp = $k_imp*$imp_fac) evaluate ($knoe = $knoe*$noe_fac) evaluate ($kdani = $kdani*$dani_fac) evaluate ($krama = $krama*$rama_fac) constraints interaction (all) (all) weights * 1 angles $k_ang improper $k_imp end end parameter nbonds cutnb=4.5 rcon=$k_vdw nbxmod=4 repel=$radius end end noe scale all $knoe end dani class relax force $kdani end dynamics internal timestep 0.002 endt 3 tbath $bath nprint $nstep ntrfr 1 end end loop cool mini powell nstep= 500 nprint= 50 end print threshold=0.5 noe evaluate ($rms_noe=$result) evaluate ($violations_noe=$violations) print threshold=5. cdih evaluate ($rms_cdih=$result) evaluate ($violations_cdih=$violations) print thres=0.05 bonds evaluate ($rms_bonds=$result) print thres=5. angles evaluate ($rms_angles=$result) print thres=5. impropers evaluate ($rms_impropers=$result) couplings print threshold 1.0 class phi end evaluate ($rms_coup = $result) evaluate ($end_viols = $violations) couplings print threshold 1.0 class gly end evaluate ($rms_coup_g = $result) evaluate ($end_viols_g = $violations) carbon print threshold = 1.0 end evaluate ($rms_ashift = $rmsca) evaluate ($rms_bshift = $rmscb) evaluate ($viol_shift = $violations) ! angledb print threshold=0.5 type torsion end ! evaluate ($rms_tordb = $rms) ! angledb print threshold=0.5 type angle end ! evaluate ($rms_angdb = $rms) dani print threshold=0.3 all end evaluate ($rms_dani=$result) evaluate ($violations_dani=$violations) remarks coeff: 13.16 1.71 0.30 remarks =============================================================== remarks overall,bonds,angles,improper,vdw,cdih,noe,coup, shift, dani remarks energies: $ener, $bond, $angl, $impr, $vdw, $cdih, $noe, $coup, $carb, remarks $dani remarks =============================================================== remarks bonds,angles,impropers,cdih,noe,coup,dani remarks bonds etc: $rms_bonds,$rms_angles,$rms_impropers,$rms_cdih,$rms_noe,$rms_coup, remarks dani etc: $rms_dani remarks shifts RMS a, b: $rms_ashift, $rms_bshift remarks =============================================================== remarks cdih end_coup end_coup_gly noe dani ! remarks violations : $violations_cdih $end_viols $violations_noe remarks violations : $violations_cdih $end_viols $end_viols_g $violations_noe $violations_dani remarks shifts: $viol_shift remarks =============================================================== remarks jcoup stats: end_rms end_rms_gly ! remarks rms-d: $rms_coup remarks rms-d: $rms_coup $rms_coup_g remarks =============================================================== {====>} {*Name(s) of the family of final structures.*} evaluate ($file = "e1_t1t2_rhomby_" + encode($count) + ".pdb") write coor output= $file end end loop structure stop