• ## Outreach

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: