### TcB Tulane university 3/2004
###
###  SCRIPTNAME:  DashOfSalt.vmd
###
############################################
### DashOfSalt.vmd  
###   replaces random bulk water molecules 
###   with  Na-Cl ion pairs
###   NOTE: this does NOT change the charge 
###		of the system, but could
###             be easily adapted to do this
#############################################
##   Overview:
###   select bulk water (defined as water greater than 5A from not water)
##  LOOP for number of ion pairs to adD
###   a) pick a random water in  bulk
###   b) replace it with a Na+ ion
###   c) redo the selection
###   d) pick a random water in bulk
###   e) replace it with a Cl- ion
### END LOOP
############################################
### USAGE:  
############################################
###	start vmd 
###   in the command window type these commands
###     set inputpdb "yoursolvated.pdb"
###     set outputpdb "yoursolvatedSalted.pdb"
###     source DashOfSalt.vmd
##
############################################
##  AlT.USAGE: 
##     change $inputpdb to "yoursolvated.pdb" 
##     change $outputpdb to "SaltySolvated.pdb" 
###        as indicated in the comments below then
##     vmd -e  "DashOfSalt.vmd"
############################################
##  INPUTS:  define $inputpdb as the input pdb
##  OUTPUTS:  writes $outputpdb as the output pdb
##     note: there  will be some loss of accuracy in the atom positions in last digit
##     because of rounding
##

############################################
## CAUTIONS
############################################
### BE VERY CAREFUL ABOUT SELECTION OF ATOMS by RESID since   the numbers as they appear in the 
#### pdb file are not necessarily unique especially for large systems OR  is different segnames
####  are used with each beginning at resid 1
####   the residue command seems to work as expected though

####################################################################
##  DEFINE INPUTS: define input/outputs here if  run as vmd -e DashOfSalt.vmd
####################################################################
#set inputpdb "solvated.pdb"
#set outputpdb "salty-solvated.pdb"


mol load pdb $inputpdb
puts " -------------     Adding a Dash of Salt ------------------------- "

## uncomment all lines containing fpna or fpcl and 
## vmd will create a file that contains information about the Na and Cl ions 
## that is created
#set fpna [ open "na.resnames" w ] 
#set fpcl [ open "cl.resnames" w ] 

## select the oxygen atom of the bulk waters
set bulk [atomselect top " water and not hydrogen and not within 5 of not water " ] 
set oxygenlist  [ $bulk get { index } ]
set numwaters [llength $oxygenlist ]
set randomindex [ expr round (rand() * $numwaters) ] 
puts " number of waters: $numwaters  random index selected: $randomindex "
set numpairs [ expr round ($numwaters * 0.15 /55.5) ]
puts "   Number of ion pairs to add $numpairs "

#  initialize the lists and data values
set replaced {  }
set nareslist { }
set clreslist { }
set all [atomselect top "all" ]
$all set  beta 0
$all set  occupancy  0

############################################
### MAIN LOOP ####
############################################
for {set i 0 } { $i < $numpairs } { incr i 1 } { 
	puts "placing bulk ion pair number $i of $numpairs"
#  get the index for all  bulk waters (tagged by oxygen)
        set oxygenlist  [ $bulk get { index } ]
	set numwaters [llength $oxygenlist ] 
	puts "   number of waters remaining in bulk $numwaters " 
#  pick a water
	set randomindex [ expr round (rand() * $numwaters) ] 
	set oxygenindex [ lindex $oxygenlist $randomindex ] 
	puts "   selected oxygen with index $oxygenindex "
	set selectoxygen [atomselect top " same residue as index $oxygenindex " ] 
	set selectoxygenatom [atomselect top " index $oxygenindex " ]
#  rename it to Sodium; and tag the  hydrogens for deletion
	$selectoxygen set beta  { 999 999 999 }  
	$selectoxygenatom set  {name resname beta }  {{ "Na+" "Na+" 1 }}  
	$selectoxygenatom set  name  "Na+"   
#  output what I did: comment/uncomment as you like
	set replaceoxygen [ $selectoxygenatom get {name  residue resname x y z } ] 
	puts "   NA replacement $replaceoxygen " 
#	puts $fpna " $replaceoxygen  "
# update the list of replaced atoms
	lappend replaced  $oxygenindex 
	lappend nareslist $oxygenindex
        set bulk [atomselect top " water and not hydrogen and not within 5 of (index $replaced or not water) " ] 

#  now replace a  water with a Chlorine
#  get the index for all  bulk waters (tagged by oxygen)
        set oxygenlist  [ $bulk get { index } ]
	set numwaters [llength $oxygenlist ] 
	puts "   number of waters remaining in bulk $numwaters " 
#  pick a water
	set randomindex [ expr round (rand() * $numwaters) ] 
	set oxygenindex [ lindex $oxygenlist $randomindex ] 
	puts "   selected oxygen with index $oxygenindex "
	set selectoxygen [atomselect top " same residue as index $oxygenindex " ] 
	set selectoxygenatom [atomselect top " index $oxygenindex " ]
#  rename it to Chlorine; and tag the  hydrogens for deletion
	$selectoxygen set beta  { 999 999 999 }  
	$selectoxygenatom set { name resname beta } {{ "Cl-" "Cl-" -1 } }
#  output what I did: comment/uncomment as you like
	set replaceoxygen [ $selectoxygenatom get {name  residue resname x y z } ] 
	puts "   CL replacement $replaceoxygen " 
#	puts $fpcl " $replaceoxygen  "
# update the list of replaced atoms
	lappend replaced  $oxygenindex 
	lappend clreslist $oxygenindex
        set bulk [atomselect top " water and not hydrogen and not within 5 of (index $replaced or not water) " ] 
}

############################################
### END MAIN LOOP ####
############################################

#close   $fpna
#close   $fpcl

set finalout [atomselect top "beta < 900" ]
$finalout writepdb  $outputpdb

puts " -------------     Salt Added.  Add pepper as desired. ------------------------- "

quit