# hello emacs, this is -*- tcl -*- code.
# 
# $Id: nacl-follow.vmd,v 1.1 2003/07/22 09:46:55 wwwadmin Exp $
# Time-stamp: <akohlmey 20.07.2003 17:41:42 timburton.bochum>
#
# Copyright (c) 2003 by <Axel.Kohlmeyer@theochem.ruhr-uni-bochum.de>
#
############################################################
# display settings
color Display {Background} silver
display projection   Orthographic
display nearclip set 0.000000
display farclip  set 10.000000

# load trajectory:
mol new {800h2o-16nacl.xyz} type xyz waitfor all
mol addfile {800h2o-16nacl.dcd} type dcd waitfor all
# delete default representation 
mol delrep 0 top

# selection macros this is the first part of the magic
atomselect macro hyd {name H}
atomselect macro oxy {name O}
# first and second solvation shell of Na+ and oxygens in transition
atomselect macro naox1 {oxy within 3.3 of index 2405}
atomselect macro naox2 {(oxy within 5.5 of index 2405) and not naox1}
atomselect macro naotr {(oxy within 3.39 of index 2405) 
    and not (oxy within 3.21 of index 2405)}
atomselect macro nahy1 {hyd within 1.31 of naox1}
atomselect macro nahy2 {hyd within 1.31 of naox2}
# first and second solvation shell of Na+ and oxygens in transition
atomselect macro clox1 {oxy within 3.9 of index 2431}
atomselect macro clox2 {(oxy within 6.2 of index 2431) and not clox1}
atomselect macro clotr {(oxy within 3.95 of index 2431)
    and not (oxy within 3.85 of index 2431)}
atomselect macro clhy1 {hyd within 1.31 of clox1}
atomselect macro clhy2 {hyd within 1.31 of clox2}
# all oxys and hydrogens in the two shells around the two ions
atomselect macro oxtot {(oxy within 5.5 of index 2405) or 
    (oxy within 6.2 of index 2431)}
atomselect macro hytot {hyd within 1.31 of oxtot}


# the ions are the simple: we only want to show a single NaCl pair.
# index 2405 the Na+ and index 2431 is the Cl- 
mol representation VDW 0.600000 21.000000
mol color Name
mol selection {index 2405 or index 2431}
mol material Opaque
mol addrep top

# oxygen atoms in the first hydration shell of the Na+
mol representation VDW 0.300000 20.000000
mol color Name
mol material Opaque
mol selection {naox1}
mol addrep top
mol representation VDW 0.250000 20.000000
# now the hydrogen atoms 
mol selection {nahy1}
mol addrep top
# now the bonds between hydrogen atoms and oxygens
mol selection {naox1 or nahy1}
mol representation DynamicBonds 1.310000 0.150000 15.000000
mol addrep top

# oxygen atoms in the first hydration shell of the Cl-
mol representation VDW 0.300000 20.000000
mol color Name
mol material Opaque
mol selection {clox1}
mol addrep top
mol representation VDW 0.250000 20.000000
# now the hydrogen atoms 
mol selection {clhy1}
mol addrep top
# now the bonds between hydrogen atoms and oxygens
mol selection {clox1 or clhy1}
mol representation DynamicBonds 1.310000 0.150000 15.000000
mol addrep top

# the second hydration shell is made transparent and with different colors
mol material Transparent
# again we start with the Na+ shell
mol color ColorID 10
mol selection {naox2}
mol representation VDW 0.300000 20.000000
mol addrep top
mol color Name
mol selection {nahy2}
mol representation VDW 0.250000 20.000000
mol addrep top
mol color ColorID 6
mol selection {naox2 or nahy2}
mol representation DynamicBonds 1.310000 0.150000 15.000000
mol addrep top
# ... and continue with the Cl- shell
mol color ColorID 12
mol selection {clox2}
mol representation VDW 0.300000 20.000000
mol addrep top
mol color Name
mol selection {clhy2}
mol representation VDW 0.250000 20.000000
mol addrep top
mol color ColorID 6
mol selection {clox2 or clhy2}
mol representation DynamicBonds 1.310000 0.150000 15.000000
mol addrep top

# oxygen atoms in transition between the first and the second hydration 
# shell get a special highlight
mol representation VDW 0.350000 20.000000
mol material Opaque
mol color ColorID 11
mol selection {naotr}
mol addrep top
mol color ColorID 4
mol selection {clotr}
mol addrep top

# finally hydrogen bonds
# mol material Opaque
# mol color ColorID 16
# mol selection {oxtot or hytot}
# mol representation HBonds 3.000000 30.000000 4.000000
# mol addrep top
# mol selection {clox1 or index 2431}
# mol representation HBonds 4.000000 30.000000 4.000000
# mol addrep top
mol rename top {NaCl Ion Pair with Water}

# let all selections be recalculated for each frame
# and smooth the trajectory a little bit for all representations
# that's part two of the magic.
set molid 0
set n [molinfo $molid get numreps]
for {set i 0} {$i < $n} {incr i} {
    mol selupdate $i $molid on
    mol smoothrep $molid $i 2
}

# realign the ion pair to be horizontal along the 
# x-axis and keep it in the center of the screen
proc do_realign {args} {
    global chist dhist hcount hoffs molid

    # this is the axis to align to
    set avec [vecnorm {1.0 0.0 0.0}]
    # this is the sliding window size
    set asize 20

    # initialize the cache counters
    if ([info exists hcount]) { } else {
        set hcount 0
        set hoffs  -1
    }

    # calculate center and direction of the ion pair
    set sel [atomselect $molid "index 2405 or index 2431"]
    lassign [$sel get {x y z}] na cl
    set cent [vecscale [vecadd $na $cl] 0.5]
    set dir [vecsub $na $cl]

    # store data in cache for sliding window averaring
    if {$hcount < $asize} then { incr hcount }
    incr hoffs
    if {$hoffs >= $hcount} then { set hoffs 0 }
    set chist($hoffs) $cent
    set dhist($hoffs) $dir

    # calculate averages
    set csum [veczero]
    set dsum [veczero]
    for {set i 0} {$i < $hcount} {incr i} {
        set csum [vecadd $csum $chist($i)]
        set dsum [vecadd $dsum $dhist($i)]
    }
    set csum [vecscale [expr 1.0/[expr $hcount * 1.0]] $csum]
    set dsum [vecnorm $dsum]

    # get rotation axis
    set rvec [vecnorm [veccross $dsum $avec]]

    # set origin and rotation
    molinfo $molid set \
        { center_matrix rotate_matrix} \
        [list [trans origin $csum] \
             [trans axis $rvec [expr acos([vecdot $dsum $avec])] rad]]
    # clean up selections
    $sel delete

#########3
# debug code, draw current and smoothed alignment vectors
#    graphics $molid delete all
#    graphics $molid color iceblue
#    vmd_draw_vector $molid $csum $dsum 4
#    graphics $molid color white
#    vmd_draw_vector $molid $csum $rvec 4
#    graphics $molid color purple
#    vmd_draw_vector $molid $cent $dir 0.5
}

# 
scale to 0.15
trace variable vmd_frame(0) w do_realign


# go back to the start of the trajectory.
animate style rock
animate goto start

############################################################
# Local Variables:
# mode: tcl
# mode: turn-on-auto-revert
# mode: auto-fill
# time-stamp-format: "%u %02d.%02m.%y %02H:%02M:%02S %s"
# End:
############################################################