Building solvent around a protein

The following script illustrates how psfgen and VMD can be used together to add water around a protein structure. It assumes you already have a psf and pdb file for your protein, as well as a box of water which is large enough to contain the protein. For more information on how atomselections can be used within VMD scripts, see the VMD User's Guide.

proc addwater { psffile pdbfile watpsf watpdb } {
	# Create psf/pdb files that contain both our protein as well as
	# a box of equilibrated water.  The water box should be large enough
	# to easily contain our protein.
	resetpsf
	readpsf $psffile pdb $pdbfile
	readpsf $watpsf pdb $watpdb

	# Load the combined structure into VMD   
	writepsf combine.psf
	writepdb combine.pdb
	mol load psf combine.psf pdb combine.pdb

	# Assume that the segid of the water in watpsf is QQQ
	# We want to delete waters outside of a box ten Angstroms
	# bigger than the extent of the protein. 
	set protein [atomselect top "not segid QQQ"]
	set minmax [measure minmax $protein]
	foreach {min max} $minmax { break }
	foreach {xmin ymin zmin} $min { break }
	foreach {xmax ymax zmax} $max { break }
    set xmin [expr $xmin - 10]
    set ymin [expr $ymin - 10]
    set zmin [expr $zmin - 10]
    set xmax [expr $xmax + 10]
    set ymax [expr $ymax + 10]
    set zmax [expr $zmax + 10]

	# Center the water on the protein.  Also update the coordinates held
	# by psfgen.
	set wat [atomselect top "segid QQQ"]
	$wat moveby [vecsub [measure center $protein] [measure center $wat]]
	foreach atom [$wat get {segid resid name x y z}] {
		foreach {segid resid name x y z} $atom { break }
		coord $segid $resid $name [list $x $y $z]
	}

	# Select waters that we don't want in the final structure.
	set outsidebox [atomselect top "segid QQQ and (x <= $xmin or y <= $ymin \
		or z <= $zmin or x >= $xmax or y >= $ymax or z >= $xmax)"]
	set overlap [atomselect top "segid QQQ and within 2.4 of (not segid QQQ)"]

	# Get a list of all the residues that are in the two selections, and delete
	# those residues from the structure.
	set reslist [concat [$outsidebox get resid] [$overlap get resid]]
	set reslist [lsort -unique -integer $reslist]

	foreach resid $reslist {
		delatom QQQ $resid
	}

	# That should do it - write out the new psf and pdb file. 
	writepsf solvate.psf 
	writepdb solvate.pdb

	# Delete the combined water/protein molecule and load the system that
	# has excess water removed.
	mol delete top
	mol load psf solvate.psf pdb solvate.pdb

	# Return the size of the water box
	return [list [list $xmin $ymin $zmin] [list $xmax $ymax $zmax]]
}