Here we describe in detail the operation of the script siliconNitridePsf.tcl, which is used to generate the structure file for our synthetic subsystems in section 1.3.
Looking below, you'll see that the script siliconNitridePsf.tcl
has a number of parameters. Besides having the usual input and output
files, we have three flags which determine whether the script should
search for angles (findAngles), whether bonds should be formed
across hexagonal periodic boundaries in the 
-plane (hexPeriodic), and whether bonds should be formed between the
hexagonal faces at the top and bottom (zPeriodic). We want to
determine the angles and bond across the periodic boundaries; however,
we wish to have water above and below the membrane, so we do not bond
the top of the membrane to its bottom.
The next important parameter to note is bondDistance. It is the threshold distance between atoms below which bonds are created. The remaining parameters define the properties of the silicon and nitrogen atoms--such as their masses and charges. NAMD matches the atom type to values in the parameter files that give the bond and non-bonded force constants between atoms. We take a look at one of these parameter files in section 1.4.
siliconNitridePsf.tcl
# Make a psf file for Si3N4.
set fileNamePrefix membrane
# Input:
set pdbFile ${fileNamePrefix}.pdb
set boundaryFile ${fileNamePrefix}.bound
# Output:
set psfFile ${fileNamePrefix}.psf
set surfPdb surf.pdb
# Parameters:
# Should angles be calculated in addition to bonds?
set findAngles 1
set hexPeriodic 1
set zPeriodic 0
# "bondDistance" is used to determine whether a bond exists between atoms.
set bondDistance 2.0
# Si parameters
set nameSi "SI.*"
set massSi 28.085500
set chargeSi 0.7710
set typePrefixSi SI_
set numBondsSi 4
# N parameters
set nameN "N.*"
set massN 14.00700
# chargeN is determined by neutrality.
set chargeN 0.
set typePrefixN N_
set numBondsN 3
The main procedure is the driver for the script. Note that it determines the nitrogen charge to enforce charge neutrality in the system. See the box ``Charge neutrality'' in section 1.3 for more information.
proc main {} {
global pdbFile boundaryFile psfFile surfPdb
global findAngles hexPeriodic zPeriodic
global bondDistance
global nameSi massSi chargeSi typePrefixSi numBondsSi
global nameN massN chargeN typePrefixN numBondsN
set selTextSi "name \"${nameSi}\""
set selTextN "name \"${nameN}\""
# Load the pdb.
mol load pdb $pdbFile
set nAtoms [molinfo top get numatoms]
# Get the number of nitrogen and silicon atoms.
set silicon [atomselect top $selTextSi]
set numSilicon [$silicon num]
$silicon delete
set nitrogen [atomselect top $selTextN]
set numNitrogen [$nitrogen num]
$nitrogen delete
# Determine the nitrogen charge.
set chargeN [expr -$chargeSi*$numSilicon/$numNitrogen]
puts "Charge on nitrogen atoms: $chargeN"
The procedure first calls bondAtoms to find bonds between internal atoms and then finds bonds across the periodic boundaries by bonding to periodic images with bondPeriodic (Fig. 3). The location of the periodic images are obtained by extracting information from the boundary file with readRadius and readLz. To save time, we do not attempt to bond all atoms to the periodic images, only those that did not receive a complete set of bonds (defined by numBondSi and numBondsN) during the first bonding step. To accomplish this, a temporary PDB file, surf.pdb, is created in which all the atoms that are incompletely bonded are marked 0.0 in the B column of the PDB. The procedure bondPeriodic is used to search for the bonds. If both hexPeriodic and zPeriodic are not true, then some atoms will never have a complete set of bonds. These are the true surface atoms--those that will be in contact with water molecules in the simulations. Next we put the bond lists in a more convenient form. We then call findAngles and finally write the PSF file with manifestPsf.
# Find the internal bonds.
puts "Bonding internal atoms..."
set bond [bondAtoms all $bondDistance]
puts "Internal bonds: [expr [llength $bond]/4]"
# Bond to periodic images.
if {$hexPeriodic || $zPeriodic} {
# Create the surface atom pdb.
set all [atomselect top all]
$all set beta 1.0
puts "Searching for surface atoms..."
set nSurfSi [markSurface $bond $selTextSi $numBondsSi]
set nSurfN [markSurface $bond $selTextN $numBondsN]
puts "Number of surface silicons: $nSurfSi"
puts "Number of surface nitrogens: $nSurfN"
$all writepdb $surfPdb
$all delete
# Load it up.
mol delete top
mol load pdb $surfPdb
if {$hexPeriodic} {
puts "The system has hexagonal periodic boundary conditions."
set radius [readRadius $boundaryFile]
puts "Hexagon radius: $radius"
# Determine the centers of the image hexagons.
set pi [expr 4.0*atan(1.0)]
set hexCen {}
set d [expr sqrt(3.)*$radius]
for {set i 0} {$i < 6} {incr i} {
set theta [expr $pi/6.*(2*$i-1)]
lappend hexCen [list [expr $d*cos($theta)] \
[expr $d*sin($theta)] 0.]
}
puts "Periodic image displacements: $hexCen"
# Find the bonds on the periodic boundaries.
puts "Bonding to the periodic image..."
foreach r $hexCen {
set bond [concat $bond \
[bondPeriodic all $bondDistance $r]]
}
}
if {$zPeriodic} {
puts "The system is periodic along the z-axis."
set lz [readLz $boundaryFile]
puts "Period in z: $lz"
set zCen [list [list 0 0 -${lz}] [list 0 0 $lz]]
# Find the bonds on the periodic boundaries.
puts "Bonding to the periodic image..."
foreach r $zCen {
set bond [concat $bond \
[bondPeriodic all $bondDistance $r]]
}
}
}
mol delete top
puts "Counting bonds on each atom..."
countBonds count $bond $nAtoms
puts "Reorganizing bond lists..."
set bond [reorganizeBonds $bond]
puts "Removing redundancy..."
set bond [removeRedundantBonds $bond]
set totalBonds [llength $bond]
puts "Number of bonds: $totalBonds"
set angle {}
if {$findAngles} {
puts "Determining the angles..."
set angle [findAngles $bond]
set totalAngles [llength $angle]
puts "Number of angles: $totalAngles"
}
puts "Writing psf file..."
manifestPsf $psfFile $pdbFile $nAtoms bond angle count
puts "The file $psfFile was written successfully."
}
The procedure bondAtoms uses VMD's atom selection interface to find other atoms within bondDistance of each atom. Because the procedure searches for neighbors of each atom, the resulting list contains each bond twice, since a bond between atom 1 and atom 2 is the same as a bond between atom 2 and atom 1. Note that the algorithm has a quadratic growth rate in the number of atoms. For the small systems in this tutorial, the method used here should be fast enough. However, by effecting a spatial decomposition of the system, we could reduce the growth rate to linear in the number of atoms.
# Find bonds between internal atoms and return them.
proc bondAtoms {selText bondDistance} {
set sel [atomselect top $selText]
set pos [$sel get {x y z}]
set index [$sel get index]
$sel delete
set bondDistance2 [expr $bondDistance*$bondDistance]
set bond {}
foreach r $pos ind $index {
# Select neighboring atoms.
foreach {x y z} $r { break }
set nearText "($x-x)^2+($y-y)^2+($z-z)^2 < $bondDistance2"
set near [atomselect top \
"$selText and $nearText and not index $ind"]
set nearNum [$near num]
set nearIndex [$near get index]
$near delete
# Add them to the bond list.
foreach i $nearIndex {lappend bond $ind $i}
}
return $bond
}
The following two procedures extract information about the system's geometry from the boundary file for use in bonding across periodic boundaries.
# Get the radius from the boundary file.
proc readRadius {boundaryFile} {
set in [open $boundaryFile r]
foreach line [split [read $in] \n] {
if {[string match "radius *" $line]} {
set radius [lindex $line 1]
break
}
}
close $in
return $radius
}
# Get the cellBasisVector3_z from the boundary file.
proc readLz {boundaryFile} {
set in [open $boundaryFile r]
foreach line [split [read $in] \n] {
if {[string match "cellBasisVector3 *" $line]} {
set lz [lindex $line 3]
break
}
}
close $in
return $lz
}
The procedure bondPeriodic acts much like bondAtoms except that the entire system is shifted to its periodic image using VMD's moveby command (Fig. 3).
# Try to bond surface atoms to the periodic image.
proc bondPeriodic {selText bondDistance periodicDisp} {
set selText "$selText and beta == 0.0"
set sel [atomselect top $selText]
set pos [$sel get {x y z}]
set index [$sel get index]
# Shift all of the atoms into this periodic image.
$sel moveby $periodicDisp
set bondDistance2 [expr $bondDistance*$bondDistance]
set bond {}
foreach r $pos ind $index {
# Select neighboring atoms.
foreach {x y z} $r { break }
set nearText "($x-x)^2+($y-y)^2+($z-z)^2 < $bondDistance2"
set near [atomselect top \
"$selText and $nearText and not index $ind"]
set nearNum [$near num]
set nearIndex [$near get index]
$near delete
# Add them to the bond list.
foreach i $nearIndex {lappend bond $ind $i}
}
# Return all atoms to their original position.
$sel set {x y z} $pos
$sel delete
return $bond
}
The following procedure sets the beta value to 0.0 for all atoms that do not have a full set of bonds. This includes both atoms that will later be bonded to the periodic image and those that are truly on the surface. Marking these atoms allows most of the atoms to be skipped when bonding to the periodic images.
# Find the atoms that have fewer than "numBonds" bonds.
# Mark surface atoms by beta = 0.0.
# Warning! The bond list is assumed to be flat and redundant.
proc markSurface {bond selText numBonds} {
set sel [atomselect top $selText]
set index [$sel get index]
set nSurfAtoms 0
foreach i $index {
# Find the number of bonds for each atom.
set n [llength [lsearch -all $bond $i]]
# Assume each bond is in the list twice.
set n [expr $n/2]
# Set the beta value to 0.0 if the atom is on the surface.
if {$n < $numBonds} {
set s [atomselect top "index $i"]
$s set beta 0.0
incr nSurfAtoms
$s delete
}
}
$sel delete
return $nSurfAtoms
}
The procedure countBonds creates a Tcl array linking each atom to the total number of bonds that it has. We need this number because the string in the PSF type column is determined by the number of bonds. For example, the type N_2 refers to a nitrogen atom with two bonds. Consequently, we can define different bond constants in the parameter files depending on the coordination of the atom. None of the parameter files we will use here discriminate in this way, however.
# Count the number of bonds on each atom and return an array (zero-based).
# The result is placed in a variable name countVar.
# Warning! The bond list is assumed to be flat and redundant.
proc countBonds {countVar bond nAtoms} {
upvar $countVar count
set num {}
for {set i 0} {$i < $nAtoms} {incr i} {
set n [llength [lsearch -all $bond $i]]
set n [expr $n/2]
lappend num $i $n
}
array set count $num
}
The following two procedures reformat the bond lists. The first converts the flat list of bonds into nested lists containing pairs of atom indices. It also adds 1 to all of the indices since the first PSF index is 1, but VMD atom indices are 0-based. The second of the two procedures removes bonds that are permutations of one another, as mentioned earlier.
# Put the bonds into sublists.
# Reindex to a 1-based index.
proc reorganizeBonds {bond} {
set ret {}
foreach {b0 b1} $bond {
incr b0
incr b1
lappend ret [list $b0 $b1]
}
return $ret
}
# We should now have all of the bonds twice.
# Find the unique bonds.
proc removeRedundantBonds {bond} {
set ret {}
foreach b $bond {
set bPerm [list [lindex $b 1] [lindex $b 0]]
set match [lsearch $ret $bPerm]
# Add the bond to "ret" only if it is unique.
if {$match == -1} {lappend ret $b}
}
return $ret
}
The findAngles procedure searches through all unique pairs of
bonds and finds triplets of atoms such that atom A is bonded
to atom B and atom B is bonded to atom
C. Since each atom not on the surface has a fixed number of
bonds, the number of bonds is proportional to number of atoms. Thus,
the algorithm is in 
where 
is the number of atoms. We
could reduce the asymptotic complexity by writing the algorithm more
cleverly; however, for our purposes here this method is fast enough.
Because of the quadratic complexity and the fact that it is written
entirely in Tcl--it uses no fast built-in VMD commands--this
procedure can take a long to time run for large systems.
# Find the angles.
proc findAngles {bond} {
set totalBonds [llength $bond]
set totalBonds1 [expr $totalBonds - 1]
# Find bonds that share atoms.
set angle {}
for {set i 0} {$i < $totalBonds1} {incr i} {
for {set j [expr $i+1]} {$j < $totalBonds} {incr j} {
foreach {a0 a1} [lindex $bond $i] {break}
foreach {b0 b1} [lindex $bond $j] {break}
if {$a0 == $b0} {
lappend angle [list $a1 $a0 $b1]
} elseif {$a0 == $b1} {
lappend angle [list $a1 $a0 $b0]
} elseif {$a1 == $b0} {
lappend angle [list $a0 $a1 $b1]
} elseif {$a1 == $b1} {
lappend angle [list $a0 $a1 $b0]
}
}
}
return $angle
}
The final procedure writes all of the information we have determined
thus far to a PSF file. There are three sections of the PSF format
important for our
systems. The first is the atom record section
which replicates much of the identifying information contained in the
PDB as well the atom's type, mass, and charge, which are essential for
simulations. The next section of the PSF contains the bonds. The bonds
are stored in eight columns of indices, with each pair of columns in a
row representing a single bond between two atoms. Hence, each line of
the bond section of the PSF describes four bonds (except the last, which
may not be full). The final section of import to us is the angles
section which contains nine columns of indices, which as groups of
three define three bonds in each row.
# Write the psf file.
proc manifestPsf {psfFile pdbFile nAtoms bondVar angleVar countVar} {
global nameSi massSi chargeSi typePrefixSi numBondsSi
global nameN massN chargeN typePrefixN numBondsN
# Import the big pass-by-reference stuff.
upvar $bondVar bond
upvar $angleVar angle
upvar $countVar count
set dummy " 0"
set totalBonds [llength $bond]
set totalAngles [llength $angle]
set out [open $psfFile w]
##### HEADER
puts $out "PSF"
puts $out ""
puts $out " 1 !NTITLE"
puts $out " REMARKS original generated structure x-plor psf file"
##### ATOMS
puts "Writing atom records..."
puts $out ""
puts $out "[format %8i $nAtoms] !NATOM"
# Open the pdb to extract the atom records.
set inStream [open $pdbFile r]
set atom 1
foreach line [split [read $inStream] \n] {
set string0 [string range $line 0 3]
if {![string match $string0 "ATOM"]} {continue}
# Extract each pdb field.
set record [string range $line 0 5]
set serial [string range $line 6 10]
set name [string range $line 12 15]
set altLoc [string range $line 16 16]
set resName [string range $line 17 19]
set chainId [string range $line 21 21]
set resId [string range $line 22 25]
set iCode [string range $line 26 26]
set x [string range $line 30 37]
set y [string range $line 38 45]
set z [string range $line 46 53]
set occupancy [string range $line 54 59]
set beta [string range $line 60 65]
set segName [string range $line 72 75]
set element [string range $line 76 77]
set charge [string range $line 78 79]
# Determine the type names.
set numBonds $count([expr $atom-1])
set typeSi ${typePrefixSi}${numBonds}
set typeN ${typePrefixN}${numBonds}
# Write the atom record.
puts -nonewline $out [format "%8i " $atom]
puts -nonewline $out [format "%-4s " $segName]
puts -nonewline $out [format "%-4i " $resId]
puts -nonewline $out [format "%-3s " $resName]
puts -nonewline $out [format "%-4s " $name]
if {[regexp $nameSi $name]} {
puts -nonewline $out [format "%-4s " $typeSi]
puts -nonewline $out [format "% 5.6f " $chargeSi]
puts -nonewline $out [format "%6.4f " $massSi]
} else {
puts -nonewline $out [format "%-4s " $typeN]
puts -nonewline $out [format "% 5.6f " $chargeN]
puts -nonewline $out [format "%6.4f " $massN]
}
puts $out $dummy
incr atom
}
close $inStream
puts $out ""
##### BONDS
# Write the bonds.
set total [format %8i $totalBonds]
puts $out "$total !NBOND: bonds"
set num 0
foreach b $bond {
puts -nonewline $out [format "%8i%8i" [lindex $b 0] [lindex $b 1]]
incr num
if {$num == 4} {
puts $out ""
set num 0
}
}
puts $out ""
##### ANGLES
# Write the angles.
puts $out "[format %8i $totalAngles] !NTHETA: angles"
set num 0
foreach a $angle {
puts -nonewline $out \
[format "%8i%8i%8i" [lindex $a 0] [lindex $a 1] [lindex $a 2]]
incr num
if {$num == 3} {
puts $out ""
set num 0
}
}
puts $out ""
# Write everything else.
##### DIHEDRALS
set nDihedrals 0
puts $out ""
puts $out "[format %8i $nDihedrals] !NPHI: dihedrals"
puts $out ""
##### IMPROPERS
set nImpropers 0
puts $out ""
puts $out "[format %8i $nImpropers] !NIMPHI: impropers"
puts $out ""
##### DONORS
set nDonors 0
puts $out ""
puts $out "[format %8i $nDonors] !NDON: donors"
puts $out ""
##### ACCEPTORS
set nAcceptors 0
puts $out ""
puts $out "[format %8i $nAcceptors] !NACC: acceptors"
puts $out ""
##### NON-BONDED
set nNB 0
puts $out ""
puts $out "[format %8i $nNB] !NNB"
puts $out ""
set tmp [expr int($nAtoms/8)]
set tmp2 [expr $nAtoms -$tmp*8]
for {set i 0} {$i <$tmp} {incr i} {
puts $out " 0 0 0 0 0 0 0 0"
}
set lastString ""
for {set i 0} {$i <$tmp2} {incr i} {
set lastString "${lastString} 0"
}
puts $out $lastString
####### GROUPS
puts $out ""
puts $out " 1 0 !NGRP"
puts $out " 0 0 0"
puts $out ""
puts $out ""
close $out
}
main