From: Chaloin Laurent (laurent.chaloin_at_univ-montp1.fr)
Date: Wed Oct 07 2009 - 05:22:01 CDT
Hi all,
I am trying to minimize under distance constraints a protein dimer. Thus, I
used tcl scripting interface in order to set up three distances between the
two monomers and performed energy minimization. It seemed to work at the
beginning but when the minimum of total energy was reached, nothing evolves
afterwards even when I applied a new force (another distance constraint). Does
it mean that I have to carry out short dynamics runs followed by an energy
mini step or that my tcl script is not good at all ? (see below)
Thanks for your precious help.
Laurent
# Scripting for setting up several distance restraints during DYNAMICS
# switch on tclforces
tclforces on
# Call tcl force
tclForcesScript {
# Set atom index (index number .coor file not from VMD)
set at1 1242
set at2 4519
set at3 7796
set at4 17627
set at5 11073
set at6 14350
addatom $at1
addatom $at2
addatom $at3
addatom $at4
addatom $at5
addatom $at6
# Set constant force and eq. dist. V=-0.5*K(r-ro)^2
set k 1000.0
set r0 10
set tclfreq 500
set ts 0
print " TCLFORCE ACTION "
# Call procedure
proc calcforces {} {
# set change variable
global tclfreq
global ts
global at1 at2 k r0
loadcoords c
set r12 [vecsub $c($at2) $c($at1) ]
set ra [veclength $r12 ]
global at3 at4 k r0
loadcoords c
set r34 [vecsub $c($at4) $c($at3) ]
set rb [veclength $r34 ]
global at5 at6 k r0
loadcoords c
set r56 [vecsub $c($at6) $c($at5) ]
set rc [veclength $r56 ]
# set unit vector
set n0 [vecnorm $r12]
set n1 [vecnorm $r34]
set n2 [vecnorm $r56]
# Force: F=-dV/dr
set forcea [expr $k*($ra-$r0)]
set forceaa [expr -$k*($ra-$r0)]
set forceb [expr $k*($rb-$r0)]
set forcebb [expr -$k*($rb-$r0)]
set forcec [expr $k*($rc-$r0)]
set forcecc [expr -$k*($rc-$r0)]
# set vector force for atom 1 to 6
set f1 [vecscale $forcea $n0 ]
set f2 [vecscale $forceaa $n0 ]
set f3 [vecscale $forceb $n1 ]
set f4 [vecscale $forcebb $n1 ]
set f5 [vecscale $forcec $n2 ]
set f6 [vecscale $forcecc $n2 ]
# add force f1 and f2 to atoms 1 and 2 and also up to atom 6
addforce $at1 $f1
addforce $at2 $f2
addforce $at3 $f3
addforce $at4 $f4
addforce $at5 $f5
addforce $at6 $f6
if { [expr $ts % $tclfreq == 0]} {
#print forces
loadforces f
foreach {atomid force} [array get f] {
puts "FORCEext:$atomid $ts $force"
}
}
incr ts
return
set ts [expr $ts+2]
}
# There is a definition of subroutine veclenght and vecnorm
# Returns: the vector length
proc veclength {v} {
set retval 0
foreach term $v {
set retval [expr $retval + $term * $term ]
}
return [expr sqrt($retval)]
}
# Returns: the normal vector
proc vecnorm {v} {
set sum 0
foreach term $v {
set sum [expr $sum + $term*$term]
}
set sum [expr sqrt($sum)]
set retval {}
foreach term $v {
lappend retval [expr $term / $sum]
}
return $retval
}
}
-- _______________________________________________________________ L. Chaloin Biophysics and Bioinformatics Centre d'études d'agents Pathogènes et Biotechnologies pour la Santé UMR5236 - CNRS - University Montpellier Institut de Biologie - 4 Bd Henri IV F-34965 Montpellier cedex 2 Tel: 33(0) 467 600 231 Fax: 33(0) 467 604 420
This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:20 CST