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

#!/usr/bin/env perl

# 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: 1.1 $"
# date      = "$Date: 2006-10-10 14:08:10 $"
# copyright = "Copyright (c) 2006 by the University of Notre Dame"
# license   = "OOPSE"


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 oopse 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 (oopse .md file)\n";
  print "  -n char : name for the solute\n";
  print "                 (default: SOLUTE)\n";
  print "  -o char : carved solvent output file (oopse .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);

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 "<OOPSE version=4>\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</OOPSE>\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:	1063
Committed:	Tue Oct 10 14:08:10 2006 UTC (19 years ago) by chrisfen
Original Path:	*trunk/src/applications/utilities/solvator*
File size:	7759 byte(s)
Log Message:	who needs .in files for perl? using env
#	User	Rev	Content
1	chrisfen	1063	#!/usr/bin/env perl
2
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			# version = "$Revision: 1.1 $"
8			# date = "$Date: 2006-10-10 14:08:10 $"
9			# copyright = "Copyright (c) 2006 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			$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			print "solvator: carves a solute void in an oopse water box\n\n";
33			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			print " -i char : solvent input file (oopse .md file)\n";
41			print " -n char : name for the solute\n";
42			print " (default: SOLUTE)\n";
43			print " -o char : carved solvent output file (oopse .md format)\n";
44			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
62			if (defined($opt_p) \|\| defined($opt_x)){
63			} else {
64			die "Error: No solute file specifed\n Please select a solute file with the -p or -x flags (pdb or xyz respectively)\n";
65			}
66
67			if (defined($opt_d)){
68			if ($opt_d =~ /^[0-9]/) {
69			$dval = $opt_d;
70			$d2tolerance = $dval*$dval;
71			} else {
72			die "Error: The '-d' value ($opt_d) is not a valid number\n Please choose a positive real # value\n";
73			}
74			}
75
76			# open the file writer
77			open(OUTFILE, ">./$fileName") \|\| die "Error: can't open file $fileName\n";
78
79			# open and read the pdb or xyz file and get solute coordinates
80			open(SOLUTEFILE, "$soluteFileName") \|\| die "Error: Can't open file $soluteFileName\n";
81
82			if (defined($opt_x)){
83			$headerLines = 2;
84			while (<SOLUTEFILE>){
85			if ($headerLines > 0){
86			$headerLines--;
87			} else {
88			chomp;
89			@line = split;
90			push(@solute_names, $line[0]);
91			push(@solute_x, $line[1]);
92			push(@solute_y, $line[2]);
93			push(@solute_z, $line[3]);
94			}
95			}
96			}
97			if (defined($opt_p)){
98			while (<SOLUTEFILE>){
99			chomp;
100			@line = split;
101			if ($line[0] eq 'ATOM'){
102			push(@solute_names, $line[2]);
103			push(@solute_x, $line[5]);
104			push(@solute_y, $line[6]);
105			push(@solute_z, $line[7]);
106			}
107			}
108			}
109
110			# remap solute to the center of the box
111			$xSol = 0;
112			$ySol = 0;
113			$zSol = 0;
114			for ($i=0; $i<=$#solute_x; $i++){
115			$xSol += $solute_x[$i];
116			$ySol += $solute_y[$i];
117			$zSol += $solute_z[$i];
118			}
119			$xSol /= $#solute_x + 1;
120			$ySol /= $#solute_y + 1;
121			$zSol /= $#solute_z + 1;
122			for ($i=0; $i<=$#solute_x; $i++){
123			$solute_x[$i] -= $xSol;
124			$solute_y[$i] -= $ySol;
125			$solute_z[$i] -= $zSol;
126			}
127
128			if ($opt_f){
129			$soluteCount = $#solute_x + 1;
130			} else {
131			$soluteCount = 1;
132			}
133
134			$solventCount = $soluteCount;
135			# okay, let's open the solvent box and write a carved file
136			open(SOLVENT, "$solventName") \|\| die "Error: Can't open file $solventName\n";
137
138
139			while (<SOLVENT>){
140			$startSnap = 0 if /\/Snapshot/;
141			$startFrame = 0 if /\/FrameData/;
142			$startStunts = 0 if /\/StuntDoubles/;
143			$startMeta = 0 if /\/MetaData/;
144
145			# save the MetaData lines
146			push(@metaLines, $_) if $startMeta == 1;
147
148			if ($startSnap == 1){
149			# save the frame data
150			push(@frameData, $_) if $startFrame == 1;
151			if (/Hmat/){
152			@line = split;
153			$hxx = $line[2];
154			chop($hxx);
155			$hyy = $line[8];
156			chop($hyy);
157			$hzz = $line[14];
158			}
159			if ($startStunts == 1){
160			@line = split;
161
162			# wrap positions back into the box
163			$x_val = $line[2] - ($hxx*round($line[2]/$hxx));
164			$y_val = $line[3] - ($hyy*round($line[3]/$hyy));
165			$z_val = $line[4] - ($hzz*round($line[4]/$hzz));
166
167			$saveFlag = 1;
168			for($i=0; $i<=$#solute_x; $i++){
169			$diff_x = $x_val - $solute_x[$i];
170			$diff_y = $y_val - $solute_y[$i];
171			$diff_z = $z_val - $solute_z[$i];
172			$dist2 = $diff_x$diff_x + $diff_y$diff_y + $diff_z*$diff_z;
173			if ($dist2 < $d2tolerance) {
174			$saveFlag = 0;
175			last;
176			}
177			}
178			if ($saveFlag == 1){
179			$saveLine = "$solventCount\t$line[1]\t$line[2]";
180			for ($i = 3; $i <= $#line; $i++){
181			$saveLine = join(' ', $saveLine, $line[$i]);
182			}
183			push(@goodSolventMolecules, $saveLine);
184			$solventCount++;
185			}
186			}
187			}
188			$startSnap = 1 if /Snapshot/;
189			$startFrame = 1 if /FrameData/;
190			$startStunts = 1 if /StuntDoubles/;
191			$startMeta = 1 if /MetaData/;
192			# check again since "/value" contains "value"
193			$startSnap = 0 if /\/Snapshot/;
194			$startFrame = 0 if /\/FrameData/;
195			$startStunts = 0 if /\/StuntDoubles/;
196			$startMeta = 0 if /\/MetaData/;
197			}
198
199			$nSolvent = $#goodSolventMolecules + 1;
200
201			# write out the final file
202			writeOutFile();
203			print "The solvated system \"$fileName\" was generated.\n";
204
205
206			sub writeOutFile {
207			# write out the header
208			print OUTFILE "<OOPSE version=4>\n";
209			printMetaData();
210			printFrameData();
211			print OUTFILE " <StuntDoubles>\n";
212			# now print out the atoms
213			if (defined($opt_f)){
214			for ($i=0; $i<=$#solute_x; $i++){
215			print OUTFILE "$i\tp\t$solute_x[$i] $solute_y[$i] $solute_z[$i]\n";
216			}
217			} else {
218			print OUTFILE "0\tpq\t0.0 0.0 0.0 1.0 0.0 0.0 0.0\n";
219			}
220			for ($i=0; $i<=$#goodSolventMolecules; $i++) {
221			print OUTFILE "$goodSolventMolecules[$i]\n";
222			}
223			print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n";
224			}
225
226			sub printMetaData() {
227			print OUTFILE " <MetaData>\n";
228
229			writeSoluteDescription();
230
231			for ($i = 0; $i<=$#metaLines; $i++) {
232			# reset the number of solvent molecules
233			if ($metaLines[$i] =~ /nMol/) {
234			$metaLines[$i] = " nMol = $nSolvent;\n";
235			}
236
237			print OUTFILE "$metaLines[$i]";
238			}
239
240			print OUTFILE " </MetaData>\n";
241			}
242
243			sub writeSoluteDescription {
244			# include a solute model description in the meta data region
245
246			print OUTFILE "\nmolecule{\n name = \"$soluteName\";\n\n";
247
248			if (defined($opt_f)) {
249			# for flexible solutes
250			for ($i=0; $i<=$#solute_x; $i++){
251			print OUTFILE " atom[$i]{\n type = \"$solute_names[$i]\";\n }\n";
252			}
253			print OUTFILE "}\n";
254			} else {
255			# the default is a rigid body solute
256			for ($i=0; $i<=$#solute_x; $i++){
257			print OUTFILE " atom[$i]{\n type = \"$solute_names[$i]\";\n position($solute_x[$i], $solute_y[$i], $solute_z[$i]);\n }\n";
258			}
259			print OUTFILE "\n rigidBody[0]{\n members(";
260			for ($i=0; $i<$#solute_x; $i++){
261			print OUTFILE "$i, ";
262			}
263			print OUTFILE "$#solute_x);\n }\n}\n";
264			}
265
266			# now back to the metaData output
267			print OUTFILE "\ncomponent{
268			type = \"$soluteName\";
269			nMol = 1;
270			}\n";
271
272			print "The solute model definition was included in '$fileName'.\n";
273			}
274
275			sub printFrameData {
276			print OUTFILE " <Snapshot>\n <FrameData>\n";
277
278			for($i = 0; $i<=$#frameData; $i++) {
279			print OUTFILE "$frameData[$i]";
280			}
281
282			print OUTFILE " </FrameData>\n";
283			}
284
285			sub round {
286			return int( $_[0] + 0.5 * ($_[0] <=> 0) );
287			}