| 1 |
#!/usr/bin/env perl |
| 2 |
|
| 3 |
# program that reads in a water box and solute (md), merges the two systems, |
| 4 |
# and deletes overlapping molecules |
| 5 |
|
| 6 |
# author = "Dan Gezelter" |
| 7 |
# version = "$Revision: 1.3 $" |
| 8 |
# date = "$Date: 2008-01-23 21:22:18 $" |
| 9 |
# copyright = "Copyright (c) 2007 by the University of Notre Dame" |
| 10 |
# license = "OOPSE" |
| 11 |
|
| 12 |
|
| 13 |
use Getopt::Std; |
| 14 |
|
| 15 |
$d2tolerance = 7.5625; #distance to start cutting |
| 16 |
$fileName = 'solvatedSystem.md'; |
| 17 |
$startSnap = 0; |
| 18 |
$startFrame = 0; |
| 19 |
$startStunts = 0; |
| 20 |
$startMeta = 0; |
| 21 |
$nSolvent = 0; |
| 22 |
|
| 23 |
# get our options |
| 24 |
getopts('hd:v:u:o:'); |
| 25 |
|
| 26 |
# if we don't have a filename, drop to -h |
| 27 |
$opt_h = 'true' if $#ARGV != -1; |
| 28 |
|
| 29 |
# our option output |
| 30 |
if ($opt_h){ |
| 31 |
print "solvator: carves a solute void in an oopse water box\n\n"; |
| 32 |
print "usage: solvator [-h] [-d distance tolerance] [-v file name (solvent)]"; |
| 33 |
print " [-u file name (solute)] [-o file name (output)]\n\n"; |
| 34 |
print " -h : show this message\n\n"; |
| 35 |
print " -d real : overlap removal distance\n"; |
| 36 |
print " (default: 2.75 ang)\n"; |
| 37 |
print " -u char : solute input file (oopse .md file)\n"; |
| 38 |
print " -v char : solvent input file (oopse .md file)\n"; |
| 39 |
print " -o char : carved solvent output file (oopse .md format)\n"; |
| 40 |
print " (default: 'solvatedSystem.md')\n"; |
| 41 |
print "Example:\n"; |
| 42 |
die " md-solvator -v freshWater.md -u mySolute.md -o mySystem.md\n"; |
| 43 |
} |
| 44 |
|
| 45 |
# set some variables to be used in the code |
| 46 |
if (defined($opt_v)){ |
| 47 |
$solventName = $opt_v; |
| 48 |
} else { |
| 49 |
die "Error: No solvent box specified\n Please select a solvent box using the -v flag\n"; |
| 50 |
} |
| 51 |
|
| 52 |
$soluteFileName = $opt_u if defined($opt_u); |
| 53 |
$fileName = $opt_o if defined($opt_o); |
| 54 |
|
| 55 |
if (!defined($opt_u)){ |
| 56 |
die "Error: No solute file specifed\n Please select a solute file with the -u flag\n"; |
| 57 |
} |
| 58 |
|
| 59 |
if (defined($opt_d)){ |
| 60 |
if ($opt_d =~ /^[0-9]/) { |
| 61 |
$dval = $opt_d; |
| 62 |
$d2tolerance = $dval*$dval; |
| 63 |
} else { |
| 64 |
die "Error: The '-d' value ($opt_d) is not a valid number\n Please choose a positive real # value\n"; |
| 65 |
} |
| 66 |
} |
| 67 |
|
| 68 |
# open the file writer |
| 69 |
open(OUTFILE, ">./$fileName") || die "Error: can't open file $fileName\n"; |
| 70 |
|
| 71 |
# open and read the pdb or xyz file and get solute coordinates |
| 72 |
open(SOLUTEFILE, "$soluteFileName") || die "Error: Can't open file $soluteFileName\n"; |
| 73 |
|
| 74 |
|
| 75 |
while (<SOLUTEFILE>){ |
| 76 |
$startSnap = 0 if /\/Snapshot/; |
| 77 |
$startFrame = 0 if /\/FrameData/; |
| 78 |
$startStunts = 0 if /\/StuntDoubles/; |
| 79 |
$startMeta = 0 if /\/MetaData/; |
| 80 |
|
| 81 |
# save the MetaData lines |
| 82 |
push(@soluteMetaLines, $_) if $startMeta == 1; |
| 83 |
|
| 84 |
if ($startSnap == 1){ |
| 85 |
# save the frame data |
| 86 |
push(@soluteFrameData, $_) if $startFrame == 1; |
| 87 |
if (/Hmat/){ |
| 88 |
@line = split; |
| 89 |
$soluteHxx = $line[2]; |
| 90 |
chop($soluteHxx); |
| 91 |
$soluteHyy = $line[8]; |
| 92 |
chop($soluteHyy); |
| 93 |
$soluteHzz = $line[14]; |
| 94 |
} |
| 95 |
if ($startStunts == 1){ |
| 96 |
@line = split; |
| 97 |
$x_val = $line[2] - ($soluteHxx*round($line[2]/$soluteHxx)); |
| 98 |
$y_val = $line[3] - ($soluteHyy*round($line[3]/$soluteHyy)); |
| 99 |
$z_val = $line[4] - ($soluteHzz*round($line[4]/$soluteHzz)); |
| 100 |
|
| 101 |
push(@solute_x, $x_val); |
| 102 |
push(@solute_y, $y_val); |
| 103 |
push(@solute_z, $z_val); |
| 104 |
$saveLine = "$line[0]\t$line[1]\t$x_val\t$y_val\t$z_val\t"; |
| 105 |
for ($i = 5; $i <= $#line; $i++){ |
| 106 |
$saveLine = join(' ', $saveLine, $line[$i]); |
| 107 |
} |
| 108 |
push(@soluteStuntDoubles, $saveLine); |
| 109 |
$soluteCount++; |
| 110 |
} |
| 111 |
} |
| 112 |
|
| 113 |
$startSnap = 1 if /Snapshot/; |
| 114 |
$startFrame = 1 if /FrameData/; |
| 115 |
$startStunts = 1 if /StuntDoubles/; |
| 116 |
$startMeta = 1 if /MetaData/; |
| 117 |
# check again since "/value" contains "value" |
| 118 |
$startSnap = 0 if /\/Snapshot/; |
| 119 |
$startFrame = 0 if /\/FrameData/; |
| 120 |
$startStunts = 0 if /\/StuntDoubles/; |
| 121 |
$startMeta = 0 if /\/MetaData/; |
| 122 |
} |
| 123 |
|
| 124 |
$solventCount = $soluteCount; |
| 125 |
# okay, let's open the solvent box and write a carved file |
| 126 |
open(SOLVENT, "$solventName") || die "Error: Can't open file $solventName\n"; |
| 127 |
|
| 128 |
|
| 129 |
while (<SOLVENT>){ |
| 130 |
$startSnap = 0 if /\/Snapshot/; |
| 131 |
$startFrame = 0 if /\/FrameData/; |
| 132 |
$startStunts = 0 if /\/StuntDoubles/; |
| 133 |
$startMeta = 0 if /\/MetaData/; |
| 134 |
|
| 135 |
# save the MetaData lines |
| 136 |
push(@metaLines, $_) if $startMeta == 1; |
| 137 |
|
| 138 |
if ($startSnap == 1){ |
| 139 |
# save the frame data |
| 140 |
push(@frameData, $_) if $startFrame == 1; |
| 141 |
if (/Hmat/){ |
| 142 |
@line = split; |
| 143 |
$hxx = $line[2]; |
| 144 |
chop($hxx); |
| 145 |
$hyy = $line[8]; |
| 146 |
chop($hyy); |
| 147 |
$hzz = $line[14]; |
| 148 |
} |
| 149 |
if ($startStunts == 1){ |
| 150 |
@line = split; |
| 151 |
|
| 152 |
# wrap positions back into the box |
| 153 |
$x_val = $line[2] - ($hxx*round($line[2]/$hxx)); |
| 154 |
$y_val = $line[3] - ($hyy*round($line[3]/$hyy)); |
| 155 |
$z_val = $line[4] - ($hzz*round($line[4]/$hzz)); |
| 156 |
|
| 157 |
$saveFlag = 1; |
| 158 |
for($i=0; $i<=$#solute_x; $i++){ |
| 159 |
$diff_x = $x_val - $solute_x[$i]; |
| 160 |
$diff_y = $y_val - $solute_y[$i]; |
| 161 |
$diff_z = $z_val - $solute_z[$i]; |
| 162 |
$dist2 = $diff_x*$diff_x + $diff_y*$diff_y + $diff_z*$diff_z; |
| 163 |
if ($dist2 < $d2tolerance) { |
| 164 |
$saveFlag = 0; |
| 165 |
last; |
| 166 |
} |
| 167 |
} |
| 168 |
if ($saveFlag == 1){ |
| 169 |
$saveLine = "$solventCount\t$line[1]\t$x_val\t$y_val\t$z_val\t"; |
| 170 |
for ($i = 5; $i <= $#line; $i++){ |
| 171 |
$saveLine = join(' ', $saveLine, $line[$i]); |
| 172 |
} |
| 173 |
push(@goodSolventMolecules, $saveLine); |
| 174 |
$solventCount++; |
| 175 |
} |
| 176 |
} |
| 177 |
} |
| 178 |
$startSnap = 1 if /Snapshot/; |
| 179 |
$startFrame = 1 if /FrameData/; |
| 180 |
$startStunts = 1 if /StuntDoubles/; |
| 181 |
$startMeta = 1 if /MetaData/; |
| 182 |
# check again since "/value" contains "value" |
| 183 |
$startSnap = 0 if /\/Snapshot/; |
| 184 |
$startFrame = 0 if /\/FrameData/; |
| 185 |
$startStunts = 0 if /\/StuntDoubles/; |
| 186 |
$startMeta = 0 if /\/MetaData/; |
| 187 |
} |
| 188 |
|
| 189 |
$nSolvent = $#goodSolventMolecules + 1; |
| 190 |
|
| 191 |
# write out the final file |
| 192 |
writeOutFile(); |
| 193 |
print "The solvated system \"$fileName\" was generated.\n"; |
| 194 |
|
| 195 |
|
| 196 |
sub writeOutFile { |
| 197 |
# write out the header |
| 198 |
print OUTFILE "<OOPSE version=4>\n"; |
| 199 |
printMetaData(); |
| 200 |
printFrameData(); |
| 201 |
print OUTFILE " <StuntDoubles>\n"; |
| 202 |
# now print out the atoms |
| 203 |
for ($i=0; $i<=$#solute_x; $i++){ |
| 204 |
print OUTFILE "$soluteStuntDoubles[$i]\n"; |
| 205 |
} |
| 206 |
for ($i=0; $i<=$#goodSolventMolecules; $i++) { |
| 207 |
print OUTFILE "$goodSolventMolecules[$i]\n"; |
| 208 |
} |
| 209 |
print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n"; |
| 210 |
} |
| 211 |
|
| 212 |
sub printMetaData() { |
| 213 |
print OUTFILE " <MetaData>\n"; |
| 214 |
|
| 215 |
writeSoluteDescription(); |
| 216 |
|
| 217 |
for ($i = 0; $i<=$#metaLines; $i++) { |
| 218 |
# reset the number of solvent molecules |
| 219 |
if ($metaLines[$i] =~ /nMol/) { |
| 220 |
$metaLines[$i] = " nMol = $nSolvent;\n"; |
| 221 |
} |
| 222 |
|
| 223 |
print OUTFILE "$metaLines[$i]"; |
| 224 |
} |
| 225 |
|
| 226 |
print OUTFILE " </MetaData>\n"; |
| 227 |
} |
| 228 |
|
| 229 |
sub writeSoluteDescription { |
| 230 |
# include a solute model description in the meta data region |
| 231 |
|
| 232 |
for ($i = 0; $i<=$#soluteMetaLines; $i++) { |
| 233 |
print OUTFILE "$soluteMetaLines[$i]"; |
| 234 |
} |
| 235 |
print "The solute model definition was included in '$fileName'.\n"; |
| 236 |
} |
| 237 |
|
| 238 |
sub printFrameData { |
| 239 |
print OUTFILE " <Snapshot>\n <FrameData>\n"; |
| 240 |
|
| 241 |
for($i = 0; $i<=$#frameData; $i++) { |
| 242 |
print OUTFILE "$frameData[$i]"; |
| 243 |
} |
| 244 |
|
| 245 |
print OUTFILE " </FrameData>\n"; |
| 246 |
} |
| 247 |
|
| 248 |
sub round { |
| 249 |
return int( $_[0] + 0.5 * ($_[0] <=> 0) ); |
| 250 |
} |
