[ViewVC] Annotation of: OpenMD/branches/development/src/applications/utilities/solvator

#!@PERL_EXECUTABLE@

# program that reads in a water box and solute xyz (or pdb), merges the two systems, and
# deletes overlapping molecules

# author    = "Chris Fennell
# version   = "$Revision$"
# date      = "$Date$"
# copyright = "Copyright (c) 2006 by the University of Notre Dame"
# license   = "OpenMD"


use Getopt::Std;

$d2tolerance = 7.5625; #distance to start cutting
$fileName = 'solvatedSystem.md';
$startSnap = 0;
$startFrame = 0;
$startStunts = 0;
$startMeta = 0;
$soluteName = 'SOLUTE';
$nSolvent = 0;

# get our options
getopts('fhd:i:n:o:p:x:');

# if we don't have a filename, drop to -h
$opt_h = 'true' if $#ARGV != -1;

# our option output 
if ($opt_h){
  print "solvator: carves a solute void in an OpenMD water box\n\n";
  print "usage: solvator [-fh] [-d distance tolerance] [-i file name (solvent)]";
  print " [-n solute name] [-px file name (solute)] [-o file name (output)]\n\n";
  print "  -f : include a flexible solute model description in the output file rather\n";
  print "       than the rigid body solute model description\n";
  print "  -h : show this message\n\n";
  print "  -d real : overlap removal distance\n";
  print "                 (default: 2.75 ang)\n";
  print "  -i char : solvent input file (OpenMD .md file)\n";
  print "  -n char : name for the solute\n";
  print "                 (default: SOLUTE)\n";
  print "  -o char : carved solvent output file (OpenMD .md format)\n";
  print "                 (default: 'solvatedSystem.md')\n";
  print "  -p char : solute input file (pdb)\n";
  print "  -x char : solute input file (xyz)\n\n";
  print "Example:\n";
  die   "   solvator -i myWater.md -p mySolute.pdb -o mySystem.md\n";
}

# set some variables to be used in the code
if (defined($opt_i)){
  $solventName = $opt_i;
} else {
  die "Error: No solvent box specified\n       Please select a solvent box using the -i flag\n";
}

$soluteFileName = $opt_p if defined($opt_p);
$soluteFileName = $opt_x if defined($opt_x);
$fileName = $opt_o if defined($opt_o);
$soluteName = $opt_n if defined($opt_n);

if (defined($opt_p) || defined($opt_x)){
} else {
  die "Error: No solute file specifed\n       Please select a solute file with the -p or -x flags (pdb or xyz respectively)\n";
}

if (defined($opt_d)){
  if ($opt_d =~ /^[0-9]/) {
    $dval = $opt_d;
    $d2tolerance = $dval*$dval;
  } else {
    die "Error: The '-d' value ($opt_d) is not a valid number\n       Please choose a positive real # value\n";
  }
}

# open the file writer
open(OUTFILE, ">./$fileName") || die "Error: can't open file $fileName\n";

# open and read the pdb or xyz file and get solute coordinates
open(SOLUTEFILE, "$soluteFileName") || die "Error: Can't open file $soluteFileName\n";

if (defined($opt_x)){
  $headerLines = 2;
  while (<SOLUTEFILE>){
    if ($headerLines > 0){
      $headerLines--;
    } else {
      chomp;
      @line = split;
      push(@solute_names, $line[0]);
      push(@solute_x, $line[1]);
      push(@solute_y, $line[2]);
      push(@solute_z, $line[3]);
    }
  }
}
if (defined($opt_p)){
  while (<SOLUTEFILE>){
    chomp;
    @line = split;
    if ($line[0] eq 'ATOM'){
      push(@solute_names, $line[2]);
      push(@solute_x, $line[5]);
      push(@solute_y, $line[6]);
      push(@solute_z, $line[7]);
    }
  }
}

# remap solute to the center of the box
$xSol = 0;
$ySol = 0;
$zSol = 0;
for ($i=0; $i<=$#solute_x; $i++){
  $xSol += $solute_x[$i];
  $ySol += $solute_y[$i];
  $zSol += $solute_z[$i];
}
$xSol /= $#solute_x + 1;
$ySol /= $#solute_y + 1;
$zSol /= $#solute_z + 1;
for ($i=0; $i<=$#solute_x; $i++){
  $solute_x[$i] -= $xSol;
  $solute_y[$i] -= $ySol;
  $solute_z[$i] -= $zSol;
}

if ($opt_f){
  $soluteCount = $#solute_x + 1;
} else {
  $soluteCount = 1;
}

$solventCount = $soluteCount;
# okay, let's open the solvent box and write a carved file
open(SOLVENT, "$solventName") || die "Error: Can't open file $solventName\n";


while (<SOLVENT>){
  $startSnap = 0 if /\/Snapshot/;
  $startFrame = 0 if /\/FrameData/;
  $startStunts = 0 if /\/StuntDoubles/;
  $startMeta = 0 if /\/MetaData/;

  # save the MetaData lines
  push(@metaLines, $_) if $startMeta == 1;

  if ($startSnap == 1){
    # save the frame data
    push(@frameData, $_) if $startFrame == 1;
    if (/Hmat/){
      @line = split;
      $hxx = $line[2];
      chop($hxx);
      $hyy = $line[8];
      chop($hyy);
      $hzz = $line[14];
    }
    if ($startStunts == 1){
      @line = split;

      # wrap positions back into the box
      $x_val = $line[2] - ($hxx*round($line[2]/$hxx));
      $y_val = $line[3] - ($hyy*round($line[3]/$hyy));
      $z_val = $line[4] - ($hzz*round($line[4]/$hzz));

      $saveFlag = 1;
      for($i=0; $i<=$#solute_x; $i++){
        $diff_x = $x_val - $solute_x[$i];
        $diff_y = $y_val - $solute_y[$i];
        $diff_z = $z_val - $solute_z[$i];
        $dist2 = $diff_x*$diff_x + $diff_y*$diff_y + $diff_z*$diff_z;
        if ($dist2 < $d2tolerance) {
          $saveFlag = 0;
          last;
        }
      }
      if ($saveFlag == 1){
        $saveLine = "$solventCount\t$line[1]\t$line[2]";
        for ($i = 3; $i <= $#line; $i++){
          $saveLine = join(' ', $saveLine, $line[$i]);
        }
        push(@goodSolventMolecules, $saveLine);
        $solventCount++;
      }
    }
  }
  $startSnap = 1 if /Snapshot/;
  $startFrame = 1 if /FrameData/;
  $startStunts = 1 if /StuntDoubles/;
  $startMeta = 1 if /MetaData/;
  # check again since "/value" contains "value"
  $startSnap = 0 if /\/Snapshot/;
  $startFrame = 0 if /\/FrameData/;
  $startStunts = 0 if /\/StuntDoubles/;
  $startMeta = 0 if /\/MetaData/;
}

$nSolvent = $#goodSolventMolecules + 1;

# write out the final file
writeOutFile();
print "The solvated system \"$fileName\" was generated.\n";


sub writeOutFile {
  # write out the header
  print OUTFILE "<OpenMD version=1>\n";
  printMetaData();
  printFrameData();
  print OUTFILE "    <StuntDoubles>\n";
  # now print out the atoms
  if (defined($opt_f)){
    for ($i=0; $i<=$#solute_x; $i++){
      print OUTFILE "$i\tp\t$solute_x[$i] $solute_y[$i] $solute_z[$i]\n";
    }
  } else {
    print OUTFILE "0\tpq\t0.0 0.0 0.0 1.0 0.0 0.0 0.0\n";
  }
  for ($i=0; $i<=$#goodSolventMolecules; $i++) {
    print OUTFILE "$goodSolventMolecules[$i]\n";
  }
  print OUTFILE "    </StuntDoubles>\n  </Snapshot>\n</OpenMD>\n";
}

sub printMetaData() {
  print OUTFILE "  <MetaData>\n";

  writeSoluteDescription();

  for ($i = 0; $i<=$#metaLines; $i++) {
    # reset the number of solvent molecules
    if ($metaLines[$i] =~ /nMol/) {
      $metaLines[$i] = "  nMol = $nSolvent;\n";
    }

    print OUTFILE "$metaLines[$i]";
  }

  print OUTFILE "  </MetaData>\n";
}

sub writeSoluteDescription {
  # include a solute model description in the meta data region

  print OUTFILE "\nmolecule{\n  name = \"$soluteName\";\n\n";

  if (defined($opt_f)) {
    # for flexible solutes
    for ($i=0; $i<=$#solute_x; $i++){
      print OUTFILE "  atom[$i]{\n    type = \"$solute_names[$i]\";\n  }\n";
    }
    print OUTFILE "}\n";
  } else {
    # the default is a rigid body solute
    for ($i=0; $i<=$#solute_x; $i++){
      print OUTFILE "  atom[$i]{\n    type = \"$solute_names[$i]\";\n    position($solute_x[$i], $solute_y[$i], $solute_z[$i]);\n  }\n";
    }
    print OUTFILE "\n  rigidBody[0]{\n    members(";
    for ($i=0; $i<$#solute_x; $i++){
      print OUTFILE "$i, ";
    }
    print OUTFILE "$#solute_x);\n  }\n}\n";
  }

  # now back to the metaData output
  print OUTFILE "\ncomponent{
  type = \"$soluteName\";
  nMol = 1;
}\n";

  print "The solute model definition was included in '$fileName'.\n";
}

sub printFrameData {
  print OUTFILE "  <Snapshot>\n    <FrameData>\n";

  for($i = 0; $i<=$#frameData; $i++) {
    print OUTFILE "$frameData[$i]";
  }

  print OUTFILE "    </FrameData>\n";
}

sub round {
  return int( $_[0] + 0.5 * ($_[0] <=> 0) );
}
Revision:	1646
Committed:	Mon Sep 26 13:30:00 2011 UTC (13 years, 7 months ago) by gezelter
File size:	7778 byte(s)
Log Message:	update for cmake build and install
#	User	Rev	Content
1	gezelter	1646	#!@PERL_EXECUTABLE@
2	chrisfen	1063
3			# program that reads in a water box and solute xyz (or pdb), merges the two systems, and
4			# deletes overlapping molecules
5
6			# author = "Chris Fennell
7	gezelter	1442	# version = "$Revision$"
8			# date = "$Date$"
9	chrisfen	1063	# copyright = "Copyright (c) 2006 by the University of Notre Dame"
10	gezelter	1390	# license = "OpenMD"
11	chrisfen	1063
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			$soluteName = 'SOLUTE';
22			$nSolvent = 0;
23
24			# get our options
25			getopts('fhd:i:n:o:p:x:');
26
27			# if we don't have a filename, drop to -h
28			$opt_h = 'true' if $#ARGV != -1;
29
30			# our option output
31			if ($opt_h){
32	gezelter	1390	print "solvator: carves a solute void in an OpenMD water box\n\n";
33	chrisfen	1063	print "usage: solvator [-fh] [-d distance tolerance] [-i file name (solvent)]";
34			print " [-n solute name] [-px file name (solute)] [-o file name (output)]\n\n";
35			print " -f : include a flexible solute model description in the output file rather\n";
36			print " than the rigid body solute model description\n";
37			print " -h : show this message\n\n";
38			print " -d real : overlap removal distance\n";
39			print " (default: 2.75 ang)\n";
40	gezelter	1390	print " -i char : solvent input file (OpenMD .md file)\n";
41	chrisfen	1063	print " -n char : name for the solute\n";
42			print " (default: SOLUTE)\n";
43	gezelter	1390	print " -o char : carved solvent output file (OpenMD .md format)\n";
44	chrisfen	1063	print " (default: 'solvatedSystem.md')\n";
45			print " -p char : solute input file (pdb)\n";
46			print " -x char : solute input file (xyz)\n\n";
47			print "Example:\n";
48			die " solvator -i myWater.md -p mySolute.pdb -o mySystem.md\n";
49			}
50
51			# set some variables to be used in the code
52			if (defined($opt_i)){
53			$solventName = $opt_i;
54			} else {
55			die "Error: No solvent box specified\n Please select a solvent box using the -i flag\n";
56			}
57
58			$soluteFileName = $opt_p if defined($opt_p);
59			$soluteFileName = $opt_x if defined($opt_x);
60			$fileName = $opt_o if defined($opt_o);
61	chrisfen	1068	$soluteName = $opt_n if defined($opt_n);
62	chrisfen	1063
63			if (defined($opt_p) \|\| defined($opt_x)){
64			} else {
65			die "Error: No solute file specifed\n Please select a solute file with the -p or -x flags (pdb or xyz respectively)\n";
66			}
67
68			if (defined($opt_d)){
69			if ($opt_d =~ /^[0-9]/) {
70			$dval = $opt_d;
71			$d2tolerance = $dval*$dval;
72			} else {
73			die "Error: The '-d' value ($opt_d) is not a valid number\n Please choose a positive real # value\n";
74			}
75			}
76
77			# open the file writer
78			open(OUTFILE, ">./$fileName") \|\| die "Error: can't open file $fileName\n";
79
80			# open and read the pdb or xyz file and get solute coordinates
81			open(SOLUTEFILE, "$soluteFileName") \|\| die "Error: Can't open file $soluteFileName\n";
82
83			if (defined($opt_x)){
84			$headerLines = 2;
85			while (<SOLUTEFILE>){
86			if ($headerLines > 0){
87			$headerLines--;
88			} else {
89			chomp;
90			@line = split;
91			push(@solute_names, $line[0]);
92			push(@solute_x, $line[1]);
93			push(@solute_y, $line[2]);
94			push(@solute_z, $line[3]);
95			}
96			}
97			}
98			if (defined($opt_p)){
99			while (<SOLUTEFILE>){
100			chomp;
101			@line = split;
102			if ($line[0] eq 'ATOM'){
103			push(@solute_names, $line[2]);
104			push(@solute_x, $line[5]);
105			push(@solute_y, $line[6]);
106			push(@solute_z, $line[7]);
107			}
108			}
109			}
110
111			# remap solute to the center of the box
112			$xSol = 0;
113			$ySol = 0;
114			$zSol = 0;
115			for ($i=0; $i<=$#solute_x; $i++){
116			$xSol += $solute_x[$i];
117			$ySol += $solute_y[$i];
118			$zSol += $solute_z[$i];
119			}
120			$xSol /= $#solute_x + 1;
121			$ySol /= $#solute_y + 1;
122			$zSol /= $#solute_z + 1;
123			for ($i=0; $i<=$#solute_x; $i++){
124			$solute_x[$i] -= $xSol;
125			$solute_y[$i] -= $ySol;
126			$solute_z[$i] -= $zSol;
127			}
128
129			if ($opt_f){
130			$soluteCount = $#solute_x + 1;
131			} else {
132			$soluteCount = 1;
133			}
134
135			$solventCount = $soluteCount;
136			# okay, let's open the solvent box and write a carved file
137			open(SOLVENT, "$solventName") \|\| die "Error: Can't open file $solventName\n";
138
139
140			while (<SOLVENT>){
141			$startSnap = 0 if /\/Snapshot/;
142			$startFrame = 0 if /\/FrameData/;
143			$startStunts = 0 if /\/StuntDoubles/;
144			$startMeta = 0 if /\/MetaData/;
145
146			# save the MetaData lines
147			push(@metaLines, $_) if $startMeta == 1;
148
149			if ($startSnap == 1){
150			# save the frame data
151			push(@frameData, $_) if $startFrame == 1;
152			if (/Hmat/){
153			@line = split;
154			$hxx = $line[2];
155			chop($hxx);
156			$hyy = $line[8];
157			chop($hyy);
158			$hzz = $line[14];
159			}
160			if ($startStunts == 1){
161			@line = split;
162
163			# wrap positions back into the box
164			$x_val = $line[2] - ($hxx*round($line[2]/$hxx));
165			$y_val = $line[3] - ($hyy*round($line[3]/$hyy));
166			$z_val = $line[4] - ($hzz*round($line[4]/$hzz));
167
168			$saveFlag = 1;
169			for($i=0; $i<=$#solute_x; $i++){
170			$diff_x = $x_val - $solute_x[$i];
171			$diff_y = $y_val - $solute_y[$i];
172			$diff_z = $z_val - $solute_z[$i];
173			$dist2 = $diff_x$diff_x + $diff_y$diff_y + $diff_z*$diff_z;
174			if ($dist2 < $d2tolerance) {
175			$saveFlag = 0;
176			last;
177			}
178			}
179			if ($saveFlag == 1){
180			$saveLine = "$solventCount\t$line[1]\t$line[2]";
181			for ($i = 3; $i <= $#line; $i++){
182			$saveLine = join(' ', $saveLine, $line[$i]);
183			}
184			push(@goodSolventMolecules, $saveLine);
185			$solventCount++;
186			}
187			}
188			}
189			$startSnap = 1 if /Snapshot/;
190			$startFrame = 1 if /FrameData/;
191			$startStunts = 1 if /StuntDoubles/;
192			$startMeta = 1 if /MetaData/;
193			# check again since "/value" contains "value"
194			$startSnap = 0 if /\/Snapshot/;
195			$startFrame = 0 if /\/FrameData/;
196			$startStunts = 0 if /\/StuntDoubles/;
197			$startMeta = 0 if /\/MetaData/;
198			}
199
200			$nSolvent = $#goodSolventMolecules + 1;
201
202			# write out the final file
203			writeOutFile();
204			print "The solvated system \"$fileName\" was generated.\n";
205
206
207			sub writeOutFile {
208			# write out the header
209	gezelter	1390	print OUTFILE "<OpenMD version=1>\n";
210	chrisfen	1063	printMetaData();
211			printFrameData();
212			print OUTFILE " <StuntDoubles>\n";
213			# now print out the atoms
214			if (defined($opt_f)){
215			for ($i=0; $i<=$#solute_x; $i++){
216			print OUTFILE "$i\tp\t$solute_x[$i] $solute_y[$i] $solute_z[$i]\n";
217			}
218			} else {
219			print OUTFILE "0\tpq\t0.0 0.0 0.0 1.0 0.0 0.0 0.0\n";
220			}
221			for ($i=0; $i<=$#goodSolventMolecules; $i++) {
222			print OUTFILE "$goodSolventMolecules[$i]\n";
223			}
224	gezelter	1390	print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OpenMD>\n";
225	chrisfen	1063	}
226
227			sub printMetaData() {
228			print OUTFILE " <MetaData>\n";
229
230			writeSoluteDescription();
231
232			for ($i = 0; $i<=$#metaLines; $i++) {
233			# reset the number of solvent molecules
234			if ($metaLines[$i] =~ /nMol/) {
235			$metaLines[$i] = " nMol = $nSolvent;\n";
236			}
237
238			print OUTFILE "$metaLines[$i]";
239			}
240
241			print OUTFILE " </MetaData>\n";
242			}
243
244			sub writeSoluteDescription {
245			# include a solute model description in the meta data region
246
247			print OUTFILE "\nmolecule{\n name = \"$soluteName\";\n\n";
248
249			if (defined($opt_f)) {
250			# for flexible solutes
251			for ($i=0; $i<=$#solute_x; $i++){
252			print OUTFILE " atom[$i]{\n type = \"$solute_names[$i]\";\n }\n";
253			}
254			print OUTFILE "}\n";
255			} else {
256			# the default is a rigid body solute
257			for ($i=0; $i<=$#solute_x; $i++){
258			print OUTFILE " atom[$i]{\n type = \"$solute_names[$i]\";\n position($solute_x[$i], $solute_y[$i], $solute_z[$i]);\n }\n";
259			}
260			print OUTFILE "\n rigidBody[0]{\n members(";
261			for ($i=0; $i<$#solute_x; $i++){
262			print OUTFILE "$i, ";
263			}
264			print OUTFILE "$#solute_x);\n }\n}\n";
265			}
266
267			# now back to the metaData output
268			print OUTFILE "\ncomponent{
269			type = \"$soluteName\";
270			nMol = 1;
271			}\n";
272
273			print "The solute model definition was included in '$fileName'.\n";
274			}
275
276			sub printFrameData {
277			print OUTFILE " <Snapshot>\n <FrameData>\n";
278
279			for($i = 0; $i<=$#frameData; $i++) {
280			print OUTFILE "$frameData[$i]";
281			}
282
283			print OUTFILE " </FrameData>\n";
284			}
285
286			sub round {
287			return int( $_[0] + 0.5 * ($_[0] <=> 0) );
288			}