!$Revision: 2.6 $ !$Date: 2002/07/23 16:19:27 $ !$RCSfile: re_h2o.inp,v $ {* ========================================================================== *} remarks re_h2o.inp (derived from sa_l.inp) remarks Author: Michael Nilges, Jens Linge 26.8.98 remarks seperated from ARIA1.2 and modified by Jinwon Jung 2003.2.25 remarks modified for refinement in water of monomers, dimers, RDC, metals remarks Roberto Tejero, 2011 {* ========================================================================== *} {* ============ things you want to change ================= *} evaluate ($seed = 139411524 ) {* seed for random generator *} evaluate ($count = 20 ) {* starting structure number *} evaluate ($strucfile = "or36_h2o.mtf") {* molec. topology file *} evaluate ($Noeavg = sum) {* Your choice of center, sum *} evaluate ($xtimes = 1.0) {* apply user selected scale *} evaluate ($miparm = "PARAM19") {* choice of OPLSX|PARAM19|PARMALLH6|CONTACT *} {* === data for the dynamics: cycles and timesteps ==== *} evaluate ($heatingc = "400") {* heating cycles -- 200 (doubled if RDC) *} evaluate ($hotcycle = "2000") {* hot stage cycles -- 1000 (doubled if RDC) *} evaluate ($coolcycl = "400") {* cooling cycles -- 200 (doubled if RDC) *} evaluate ($timestepHeat = 0.004 ) {* timestep at heating -- 0.004 *} evaluate ($timestepHot = 0.004 ) {* timestep at hot stage -- 0.004 *} evaluate ($timestepCool = 0.001 ) {* timestep at cooling -- 0.003 *} evaluate ($WaterWanted = "yes") {* ============= Things you _may_ want to change ============ *} evaluate ($rootname = "cnsPDB/sa_cns_") evaluate ($outname = "refinedPDB/resa_") evaluate ($outwatr = "refinedPDB_w/resa_") evaluate ($finomega = 100) evaluate ($wpdb = "TOPOWAT:boxtyp20.pdb") {* For generate_water.cns *} evaluate ($HaveDisu = "no") {* -- Getting ready for S-S bridges -- *} evaluate ($HaveNoe = "yes") ! change to "no" if no NOE data evaluate ($HaveMetalNoe = "no") ! change to "no" if no metal NOE data evaluate ($HaveHbond = "yes") ! evaluate ($HaveDih = "yes") ! change to "no" if no DIH data evaluate ($HaveRDC1 = "yes") ! evaluate ($HaveRDC2 = "yes") ! evaluate ($noe_rstrs = "or36_noe.tbl") ! File holding the NOE data evaluate ($metalnoe_rstrs = "or36_metal.tbl") ! File holding NOEs from and to metal evaluate ($hbn_rstrs = "or36_hbond.tbl") ! File for hydrogen bond evaluate ($dih_rstrs = "or36_dihe.tbl") ! File holding DIHEDRAL data evaluate ($HaveDimer = "no") evaluate ($HaveSymmetry = "no") evaluate ($Elecs = "short") ! Electrostatics cutoffs, can be long or short {* Activate if dimer *} if ( $HaveSymmetry = "yes" ) then evaluate ($ini_ncs = WNCSI) evaluate ($fin_ncs = WNCSF) evaluate ($ncs_fac = ($fin_ncs/$ini_ncs)^(1/20)) evaluate ($k_ncs = $ini_ncs) end if {* RDC data only loaded if needed *} if ( $HaveRDC1 = "yes" ) then evaluate ($ksani1_f = 0.2 ) evaluate ($ksani1_i = 0.01) evaluate ($sani1_fctr = ($ksani1_f/$ksani1_i)^(1/20)) evaluate ($axispdb1 = "ORI1_or36_20.pdb") ! file with axis COORDS evaluate ($axismtf1 = "TOPOWAT:axis1.mtf") ! file with axis MTF evaluate ($file1sani = "or36_rdc1.tbl") ! file with RDC restraints media 1 evaluate ($da1 = 6.359 ) evaluate ($rhomb1 = 0.586 ) if ( $HaveRDC2 = "yes" ) then evaluate ($ksani2_f = 0.2 ) evaluate ($ksani2_i = 0.01) evaluate ($sani2_fctr = ($ksani2_f/$ksani2_i)^(1/20)) evaluate ($axispdb2 = "ORI2_or36_20.pdb") ! file with axis COORDS evaluate ($axismtf2 = "TOPOWAT:axis2.mtf") ! file with axis MTF evaluate ($file2sani = "or36_rdc2.tbl") ! file with RDC restraints media 2 evaluate ($da2 = -6.655 ) evaluate ($rhomb2 = 0.325 ) end if end if {******************* not yet developed evaluate ($HaveJcoup = "no") evaluate ($filejcoup = "or36_Jhnha.tbl") evaluate ($HaveCaCbShifts = "no") evaluate ($fileshifts_CaCb = "or36_shiftsCACB.tbl") **********************} {* =============== things yo don't want to change ================ *} {* -- from here down you need to change nothing unless you know w hat you are doing.. -- *} {* ==== protocol starts ==== *} set seed $seed end set abort normal end topology @@TOPOWAT:protein-allhdg5-4.top @@TOPOWAT:water-allhdg5-4.top if ( $HaveRDC1 = "yes" ) then @@TOPOWAT:axis.top end if @@TOPOWAT:ion.top end evaluate ($par_nonbonded = $miparm) !read the parameter files: PARAMETER @@TOPOWAT:protein-allhdg5-4.param @@TOPOWAT:water-allhdg5-4.param if ( $HaveRDC1 = "yes" ) then @@TOPOWAT:axis.param end if @@TOPOWAT:ion.param nbonds {* load a few defaults, main load is done below *} nbxmod=5 tolerance 0.5 wmin=1.5 end end {* -- adjust NBONds depending on param being used -- *} {* -- depending on the choice of paramaters to be used here the -- most advised selection for nonbonded parameters is done. *} if ( $miparm = "OPLSX" ) then parameter if ( $Elecs = "short" ) then nbonds cutnb=9.5 ctofnb=8.5 ctonnb=6.5 wmin=1.0 nbxmod=5 atom cdiel shift repel=0.0 tolerance 0.5 eps=1.0 e14fac=0.4 inhibit 0.25 end else {* ---- long cutoffs ---- *} nbonds cutnb=12 ctonnb=5.5 ctofnb=11.0 wmin=1.0 nbxmod=5 tolerance=0.5 repel=0.0 cdie shift eps=1.0 e14fac=0.4 inhibit=0.25 end end if end {-----------------PARAM19---------------------} elseif ( $miparm = "PARAM19" ) then parameter if ( $Elecs = "short" ) then nbonds cutnb=9.5 ctofnb=8.5 ctonnb=6.5 wmin=1.0 repel=0.0 tolerance 0.5 nbxmod=5 atom cdiel shift eps=1.0 e14fac=0.4 inhibit 0.25 ! test end else {* ---- long cutoffs ---- *} nbonds cutnb=12 ctonnb=5.5 ctofnb=11.0 wmin=1.0 nbxmod=5 tolerance=0.5 repel=0.0 cdie shift eps=1.0 e14fac=0.4 inhibit=0.25 end end if end {----------------PARMALLH6--------------------} elseif ( $miparm = "PARMALLH6" ) then parameter nbonds cutnb=9.5 ctofnb=8.5 ctonnb=6.5 wmin=1.3 nbxmod=5 repel=0.8 tolerance 0.5 rexponent=2 irexponent=2 rconst=5.0 end end {---------------- CONTACT -----------------} elseif ( $miparm = "CONTACT" ) then parameter nbonds cutnb=7.0 ctofnb=6.0 ctonnb=5.5 wmin=1.3 nbxmod=5 repel=1.0 tolerance 0.5 rexponent=4 irexponent=1 rconst=16.0 end end end if {* --- begins --- *} evaluate ($pdbfile = $rootname+encode($count)+".pdb") structure reset end structure @@$strucfile {* load RDC from 1 alignment media *} if ( $HaveRDC1 = "yes" ) then @@$axismtf1 end if {* load RDC from 2 alignment media *} if ( $HaveRDC2 = "yes" ) then @@$axismtf2 end if end coor init end coor @@$pdbfile if ( $HaveRDC1 = "yes" ) then coor @@$axispdb1 end if if ( $HaveRDC2 = "yes" ) then coor @@$axispdb2 end if {* ------------ coord translation -------------- *} show max (x) ( not resn ani) evaluate ($xmax = $result) show min (x) ( not resn ani) evaluate ($xmin = $result) evaluate ($x_centr = -0.5*($xmin + $xmax) ) show max (y) ( not resn ani) evaluate ($ymax = $result) show min (y) ( not resn ani) evaluate ($ymin = $result) evaluate ($y_centr = -0.5*($ymin + $ymax) ) show max (z) ( not resn ani) evaluate ($zmax = $result) show min (z) ( not resn ani) evaluate ($zmin = $result) evaluate ($z_centr = -0.5*($zmin + $zmax) ) coord translate selection = ( not resn ani) vector = ( $x_centr $y_centr $z_centr ) end do (refx = x) (all) do (refy = y) (all) do (refz = z) (all) {* ------- generate water layer ------ *} if ( $WaterWanted = "yes" ) then do (segid = "PROT") (segid " ") @TOPOWAT:generate_water.cns do (segid = " ") (segid "PROT") end if if ( $HaveSymmetry = "yes" ) then flags include bonds angle vdw impro cdih dihe noe elec sani ncs end else flags include bonds angle vdw impro cdih dihe noe elec sani end end if {* -- write the coordinates with water to see what's going on -- *} if ( $WaterWanted = "yes" ) then if ( $count < 10 ) then evaluate ($waternam0 = $outwatr+"0"+encode($count)+"_waterIni.pdb") else evaluate ($waternam0 = $outwatr+encode($count)+"_waterIni.pdb") end if write coordinates output = $waternam0 end end if noe {* -- NOE section -- *} reset nrestraints 15000 ceiling 1000 end if ( $HaveNoe = "yes" ) then {* -- load only if present -- *} noe class dist @@$noe_rstrs end end if if ( $HaveMetalNoe = "yes" ) then {* -- load only if present -- *} noe class metal @@$metalnoe_rstrs end end if if ( $HaveHbond = "yes" ) then {* -- load only if present -- *} noe class hbond @@$hbn_rstrs end end if noe averaging * $Noeavg potential * soft scale * $xtimes sqconstant * 1.0 sqexponent * 2 soexponent * 1 rswitch * 1.0 sqoffset * 0.0 asymptote * 2.0 msoexponent * 1 masymptote * -0.1 mrswitch * 1.0 avexpo hbond 20 end if ( $HaveDisu = "yes" ) then !SSBondsAsNoe {* -- Add SSbonds as noe restraints --* } end if if ( $HaveDih = "yes" ) then restraints dihedral {* -- Read dihedral constraints -- *} reset nassign = 2000 @@$dih_rstrs end end if if ( $HaveRDC1 = "yes" ) then {* -- Read RDCs constraints -- *} sani reset nrestraints=5000 class=rdc1 potential=harmonic coeff 0.0 $da1 $rhomb1 @@$file1sani if ( $HaveRDC2 = "yes" ) then class=rdc2 potential=harmonic coeff 0.0 $da2 $rhomb2 @@$file2sani end if end end if !!!Values taken from run.cns evaluate ($k_unamb= 50 ) evaluate ($k_amb= 50 ) evaluate ($k_hbond= 50 ) {* -- evaluate scales values -- *} evaluate ($scalambi = 50 * $xtimes) evaluate ($scaldist = 50 * $xtimes) evaluate ($scalhbnd = 100 * $xtimes) evaluate ($scaldihd = 200 * $xtimes) if ( $HaveNoe = "yes" ) then noe {* -- scale and configure NOE term -- *} rswitch ambi 0.5 rswitch dist 0.5 rswitch metal 0.5 mrswitch ambi 0.5 mrswitch dist 0.5 mrswitch metal 0.5 asym ambi 0.1 asym dist 0.1 asym metal 0.1 masym ambi -0.1 masym dist -0.1 masym metal -0.1 scale ambi $scalambi scale dist $scaldist scale metal $scaldist end end if if ( $HaveHbond = "yes" ) then {* -- scale and configure for HBONS -- *} noe rswitch hbon 0.5 mrswitch hbon 0.5 asym hbon 0.1 masym hbon -0.1 scale hbond $scalhbnd end end if if ( $HaveDih = "yes" ) then restraints dihedral {* -- scale dihedral term -- *} scale=$scaldihd end end if if ( $HaveRDC1 = "yes" ) then {* -- scale RDCs term --- *} evaluate ($ksani1 = $ksani1_i) if ( $HaveRDC2 = "yes" ) then evaluate ($ksani2 = $ksani2_i) end if sani class=rdc1 force=$ksani1 if ( $HaveRDC2 = "yes" ) then class=rdc2 force=$ksani2 end if end end if {* since we do not use SHAKe, increase the water bond angle energy constant *} if ( $WaterWanted = "yes" ) then parameter angle (resn tip3) (resn tip3) (resn tip3) 500 TOKEN end end if {* reduce improper and angle force constant for some atoms *} evaluate ($kangle = 50) evaluate ($kimpro = 5) evaluate ($komega = 5) parameter angle ((not resn tip3) and not resn ANI) ((not resn tip3) and not resn ANI) ((not resn tip3) and not resn ANI) $kangle TOKEN improper (all)(all)(all)(all) $kimpro TOKEN TOKEN end {* fix the protein for initial minimization *} if ( $WaterWanted = "yes" ) then fix sele = (not resn tip3) end minimize powell nstep=40 drop=100 end end if {* release protein and restrain harmonically *} fix sele = (not all) end do (refx=x) (all) do (refy=y) (all) do (refz=z) (all) restraints harmonic exponent = 2 end do (harm = 0) (all) do (harm = 10) (not name h*) igroup interaction (all) (all) end minimize powell nstep=40 drop=10 end do (refx=x) (all) do (refy=y) (all) do (refz=z) (all) minimize powell nstep=40 drop=10 end do (refx=x) (all) do (refy=y) (all) do (refz=z) (all) do (mass = 100) (all) do (fbeta = 0) (all) do (fbeta = 20. {1/ps} ) (all) evaluate ($kharm = 50) {* Activate if dimer present *} if ( $HaveSymmetry = "yes" ) then ncs restraints initialize group equi ( _SymmSel1_ ) equi ( _SymmSel2_ ) weight = $k_ncs end ? end end if {* ------------------------ heating to 500 K ------------------------------ *} for $bath in (100 200 300 400 500) loop heat do (harm = $kharm) (not name h* ) do (vx=maxwell($bath)) (all) do (vy=maxwell($bath)) (all) do (vz=maxwell($bath)) (all) dynamics cartesian nstep=$heatingc timest=$timestepheat {ps} temperature=$bath tcoupling = true nprint=50 end evaluate ($kharm = max(0, $kharm - 4)) {* -- abort condition -- *} evaluate ($critical=$temp/$bath) if ($critical > 5. ) then display ** rerun job with smaller timestep in heating, currently is: $timestepheat stop end if do (refx=x) (all) do (refy=y) (all) do (refz=z) (all) end loop heat {* ------------------------ refinement at high T: --------------------------- *} igroup interaction (all) (all) weights * 1 dihed 2 end end do (harm = 0) (all) do (vx=maxwell($bath)) (all) do (vy=maxwell($bath)) (all) do (vz=maxwell($bath)) (all) dynamics cartesian nstep=$hotcycle timest=$timestephot {ps} temperature=$bath tcoupling = true nprint=50 end igroup interaction (all) (all) weights * 1 dihed 3 end end {* ------------------------------- cooling --------------------------------- *} evaluate ($bath = 500) while ($bath ge 25) loop cool evaluate ($kangle = min(500,$kangle * 1.2)) evaluate ($kimpro = min(500,$kimpro * 1.5)) evaluate ($komega = min($finomega,$komega * 1.5)) if ( $HaveSymmetry = "yes" ) then evaluate ($k_ncs = min($fin_ncs, $k_ncs*$ncs_fac)) end if if ( $HaveRDC1 = "yes" ) then evaluate ($ksani1 = min($ksani1_f, $ksani1 * $sani1_fctr)) sani class=rdc1 force=$ksani1 if ( $HaveRDC2 = "yes" ) then evaluate ($ksani2 = min($ksani2_f, $ksani2 * $sani2_fctr)) class=rdc2 force=$ksani2 end if end end if {* Activate if dimer present *} if ( $HaveSymmetry = "yes" ) then ncs restraints initialize group equi ( _SymmSel1_ ) equi ( _SymmSel2_ ) weight = $k_ncs end ? end end if parameter angle((not resn tip3) and not resn ANI) ((not resn tip3) and not resn ANI) ((not resn tip3) and not resn ANI) $kangle TOKEN improper (all)(all)(all)(all) $kimpro TOKEN TOKEN improper (name CA) (name C) (name N) (name CA) $komega TOKEN TOKEN improper (name CA) (name CCIS) (name N) (name CA) $komega TOKEN TOKEN end do (vx=maxwell($bath)) (all) do (vy=maxwell($bath)) (all) do (vz=maxwell($bath)) (all) dynamics cartesian nstep=$coolcycl timest=$timestepcool {ps} temperature=$bath tcoupling = true nprint=50 end {* -- abort condition -- *} evaluate ($critical=$temp/$bath) if ($critical > 5. ) then display ** rerun job with smaller timestep in cooling, currently is: $timestepcool stop end if evaluate ($bath = $bath - 25) end loop cool {* ---------------------------- end of cooling ------------------------------ *} fix selection=(resn ANI) end mini powell nstep 200 end igroup interaction (all) (all) weights * 1 vdw 5 dihed 4 end end {* -- final minimization: -- *} mini powell nstep 200 end igroup interaction (not resname TIP* ) (not resname TIP* ) end {* -- prepare names for output files -- *} if ( $count < 10 ) then evaluate ($filename= $outname+"0"+encode($count)+".pdb") evaluate ($waternam= $outwatr+"0"+encode($count)+"_waterEnd.pdb") evaluate ($viofile= $outname+"0"+encode($count)+".vio") else evaluate ($filename= $outname+encode($count)+".pdb") evaluate ($waternam= $outwatr+encode($count)+"_waterEnd.pdb") evaluate ($viofile= $outname+encode($count)+".vio") end if set display_file=$viofile end set print_file=$viofile end energy end @TOPOWAT:print_coorheader.cns display overall,bonds,angles,improper,dihe,vdw,elec,noe,cdih,sani display energies: $ener, $bond, $angl, $impr, $dihe, $vdw, $elec, $noe, $cdih, $sani do (q=1) (all) write coordinates sele= ((not resn TIP3) and not resn ANI) output =$filename end {* -- write the coordinates with water to see what's going on -- *} if ( $WaterWanted = "yes " ) then write coordinates output = $waternam end end if !! delete sele = (resn TIP3) end stop