Dear all,
I am using IVM and Pararestraints to refine one
structure. But I checked the output PDB file and the
tensors axes are not orthogonal. Can anybody help me
check my script to see where the problem comes from.
Thanks
Here is the script:
! CDS 5/2/00
!
eval ($numStructs = 20) !total number of
structures to calculate
eval ($randomSeed = 735) !random seed
eval ($use_graphics = false) ! VMD-xplor graphics are
not used.
eval ($maxResid = 500) !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:parallhdg.pro @par_axis_3.pro
@ion.param @carbohydrate.param end
structure @new.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 tensors
@tensors.tbl
class others
{====>}
@dist.tbl {*Read NOE
distance ranges.*}
@air.tbl {*Read ambiguous
distance restraints.*}
end
xpcs
nres=2000
class EU
force 50.0
coeff 2350 -1092
@pcs_protein.tbl
class ligand
force 500.0
coeff 2350 -1092
@pcshifts.tbl
end
xrdc
nres=2000
class EU
force 1
coeff 2.82 -1.31
@protein_adj.tbl
class ligand
force=10
coeff 2.82 -1.31
@lac_adj.tbl
end
restraints dihed
reset
nassign = 400
scale=50
{====>}
@dih.tbl {*Read dihedral
angle restraints.*}
end
!
! starting coords
!
coor @new.pdb
! Fix dipolar coupling tensor axis
! this applies for the powell minimization- but not
for the torsion angle
! dynamics
!constraints fix=(resid 114:250) end
flags exclude * include bond angle impr cdih 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 = 150.0) ! was 150
eval ($ini_noe = 2.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
reset
set echo off message off end
group (resid 1:155)
group (resid 250 or resid 251 and
(name C1 or name C2 or name C3 or name C4 or
name C5 or name O4))
!set up the dipolar coupling tensor axis
group (resid 500)
! hinge rotate (resid 251)
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)
!write coor output=test2.pdb 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 xpcs xrdc sani
end
!
! re-init the coords
!
! write coor output=test3.pdb end
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 MDY 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 xpcs xrdc 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))
eval ($k_imp = min
($fin_imp,$k_imp*$imp_fac_vdw))
parameter
nbonds
rcon $k_vdw
repel $radius
end
end
constraints
interaction (not resname ANI) (not resname
ANI)
! interaction (not resname ANI) (not
resname 1EU)
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 MDY force $ksani end
parameter
nbonds
rcon $k_vdw
repel $radius
end
end
constraints
interaction (not resname ANI) (not resname
ANI)
! interaction (not resname ANI) (not resname
1EU)
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
! write coor output=test4.pdb end
!
! final minimization - with fixed tensor axis
!
constraints fix (resname "ANI") end
! constraints fix (resname "GAL" or resname
"GLC") end
mini powell nstep 500 nprint 50 end
! xpcs
! foff
! son
! evaluate
($filename="1pse"+encode($count)+"."+encode($count))
! save $filename
! frun 1
! soff
! end
! xrdc
! foff
! son
! evaluate
($filename="1dip"+encode($count)+"."+encode($count))
! save $filename
! frun 1
! soff
! 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 1 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)
sani print threshold=1 class MDY end
evaluate($rms_sani=$result)
evaluate($viol_sani=$violations)
! sani print threshold=1 class gal3 end
! xpcs print threshold=0.1 all end
!evaluate($rms_pcs=$result)
!evaluate($violations_pcs=$violations)
! couplings print threshold 1.0 class phi end
! evaluate ($rms_coup = $result)
! evaluate ($viols_coup = $violations)
! evaluate ($viols_coup = $violations)
xpcs print threshold=0 all end
xrdc print threshold=0 class EU end
xrdc print threshold=0 class ligand end
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 $rms_pcs
remarks
===============================================================
remarks noe dihedrals
couplings
remarks violations : $violations_noe
$violations_cdih $viols_coup $viol_sani
$violations_pcs
remarks
===============================================================
eval ($suffix = ".pdb_nonconverged")
{* Test for convergence. *}
if ($violations_noe < 10) then
if ($violations_cdih < 2) then
eval ($suffix = ".pdb") {* Change the
suffix if converged. *}
end if
end if
evaluate ($file = "rdc_dock_nodist_" +
encode($count) + $suffix)
write coor output= $file end
eval ($count = $count + 1)
end loop structure
stop
___________________________________________________________
??????????????????
http://cn.mail.yahoo.com/promo/carnival07/
