VMD-L Mailing List

From: Mike McCallum (mmccallum_at_PACIFIC.EDU)
Date: Tue Jul 03 2018 - 15:25:42 CDT

So, I think I solved my problem, thanks to Joshua V and the other, older posters. It turned out I used some of the information Joshua sent me, but I could not get his suggested solution to work for me 100%. I used what I learned from Joshua’s code suggestion and applied it to the older scripts. One of the things I did not understand was how to select and pass the index properly from one molid to another.

The nut of the problem is: given a spherical solvated droplet, which you are allowing to boil away, calculate the (velocity-based) temperature both inside and outside the droplet, you choose the radius. This is basically in order to observe evaporative cooling. I have a lot of trouble with TCL as I am not an expert in it, and I have trouble knowing when to free up variable names and so on.

This is not doing 100% what I want, but I think it’s correct as far as it goes. I really need to send the output to a file, clean it up, etc, but that is easy.

Cheers
Mike

—
########### TCL script to calculate temperature from velocity
########### trajectorySTART #########
# rigidbonds none=0; water=1; all=2
# to begin, to get this working, let's assume water rigid bonds
#
proc getT { rigidbonds fstart num_frames radius } {
    # so the main selection should be here... droplet not vapor.
    # make sure veldcd is in mol 0, and dcd is in mol 1.
    puts "## frame, drop T, vapor T, nDof, esum"
    set fspan [expr $fstart + $num_frames]
    set Tsum 0.0;
    set notTsum 0.0;
    for { set frame $fstart } { $frame < $fspan } { incr frame } {
set alldrop [atomselect 1 all frame $frame]
set notdrop [atomselect 1 "sqr(x)+sqr(y)+sqr(z) > sqr($radius)" frame $frame]
set drop [atomselect 1 "sqr(x)+sqr(y)+sqr(z) <= sqr($radius)" frame $frame]
set all [atomselect 0 all frame $frame]
#
# nothing changes if I put "frame $frame" at the end of $velsel. ??
#
# INSIDE drop
set dropvel [atomselect 0 "index [$drop get index]" frame $frame ]
set esum 0.0;
set nDoF 0;
set c1 [[ atomselect 0 "((type 'HW' or hydrogen) or water)" ] \
   get index ]
set c2 [[ atomselect 0 "water" ] get index ]
set PDBVELFACTOR 20.45482706
foreach m [$dropvel get mass] v [$dropvel get {x y z}] id [$dropvel get index] {
   set $v [vecscale $v $PDBVELFACTOR]
   set e [expr {0.5 * $m * [vecdot $v $v]}]
   set esum [expr {$esum + $e}]
   set DoF 3
   #SHAKE?
   if { ($rigidbonds == 2) && ($id in $c1) } {
set DoF 2
   } elseif { ($rigidbonds == 1) && ($id in $c2) } {
set DoF 2
   }
   set nDoF [expr {$nDoF + $DoF}]
}
#######################
# OUTSIDE drop
set notvel [atomselect 0 "index [$notdrop get index]" frame $frame ]
set notesum 0.0;
set notnDoF 0;
set notc1 [[ atomselect 0 "((type 'HW' or hydrogen) or water)" ] \
   get index ]
set notc2 [[ atomselect 0 "water" ] get index ]
set PDBVELFACTOR 20.45482706
foreach m [$notvel get mass] v [$notvel get {x y z}] id [$notvel get index] {
   set $v [vecscale $v $PDBVELFACTOR]
   set e [expr {0.5 * $m * [vecdot $v $v]}]
   set notesum [expr {$notesum + $e}]
   set DoF 3
   #SHAKE?
   if { ($rigidbonds == 2) && ($id in $c1) } {
set DoF 2
   } elseif { ($rigidbonds == 1) && ($id in $c2) } {
set DoF 2
   }
   set notnDoF [expr {$notnDoF + $DoF}]
}
####################################
$dropvel delete
$notvel delete
set kB 0.00198657 ;# kcal/(mol.K)
set T [expr { $esum * 2.0 / ($nDoF * $kB) }]
set notT [expr { $notesum * 2.0 / ($notnDoF * $kB) }]
set Tsum [expr {$Tsum + $T}]
set notTsum [expr {$notTsum + $notT}]
puts "$frame $T $notT $nDoF $esum"
##### frame loop
    }
    set Tsum [expr {$Tsum / $num_frames}]
    set notTsum [expr {$notTsum / $num_frames}]
    puts "Average inside T: $Tsum, Average outside T: $notTsum"
    #### main loop
}
########### TCL script to calculate temperature from velocity trajectory
########### END #########

—



On Jun 30, 2018, at 10:00 AM, Mike McCallum <mmccallum_at_PACIFIC.EDU<mailto:mmccallum_at_PACIFIC.EDU>> wrote:

I made a small error typing this in the “foreach m..” The “$sel” should be “$all” (that is the way it is in my script):
foreach m [$all get mass] v [$all get {x y z}] id [$all get index] {


--
C. Michael McCallum
Professor
Department of Chemistry, UOP
mmccallum .at. pacific .dot. edu (209) 946-2636 v / (209) 946-2607 fax
On Jun 30, 2018, at 10:00 AM, Mike McCallum <mmccallum_at_PACIFIC.EDU<mailto:mmccallum_at_PACIFIC.EDU>> wrote: