[ViewVC] Annotation of: group/trunk/OOPSE-4/src/applications/utilities/md-solvator

#!/usr/bin/env perl

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

# author    = "Dan Gezelter"
# version   = "$Revision: 1.1 $"
# date      = "$Date: 2007-01-11 21:46:14 $"
# copyright = "Copyright (c) 2007 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;
$nSolvent = 0;

# get our options
getopts('hd:v:u:o:');

# 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 [-h] [-d distance tolerance] [-v file name (solvent)]";
  print "  [-u file name (solute)] [-o file name (output)]\n\n";
  print "  -h : show this message\n\n";
  print "  -d real : overlap removal distance\n";
  print "                 (default: 2.75 ang)\n";
  print "  -u char : solute input file (oopse .md file)\n";
  print "  -v char : solvent input file (oopse .md file)\n";
  print "  -o char : carved solvent output file (oopse .md format)\n";
  print "                 (default: 'solvatedSystem.md')\n";
  print "Example:\n";
  die   "   md-solvator -v freshWater.md -u mySolute.md -o mySystem.md\n";
}

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

$soluteFileName = $opt_u if defined($opt_u);
$fileName = $opt_o if defined($opt_o);

if (!defined($opt_u)){
  die "Error: No solute file specifed\n       Please select a solute file with the -u flag\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";


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

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

  if ($startSnap == 1){
    # save the frame data
    push(@soluteFrameData, $_) if $startFrame == 1;
    if (/Hmat/){
      @line = split;
      $soluteHxx = $line[2];
      chop($soluteHxx);
      $soluteHyy = $line[8];
      chop($soluteHyy);
      $soluteHzz = $line[14];
    }
    if ($startStunts == 1){
      @line = split;
      $x_val = $line[2] - ($soluteHxx*round($line[2]/$soluteHxx));
      $y_val = $line[3] - ($soluteHyy*round($line[3]/$soluteHyy));
      $z_val = $line[4] - ($soluteHzz*round($line[4]/$soluteHzz));

      push(@solute_x, $x_val);
      push(@solute_y, $y_val);
      push(@solute_z, $z_val);
      $saveLine = "";
      for ($i = 0; $i <= $#line; $i++){
        $saveLine = join(' ', $saveLine, $line[$i]);
      }
      push(@soluteStuntDoubles, $saveLine);
      $soluteCount++;
    }
  }
  $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/;
}

$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
  for ($i=0; $i<=$#solute_x; $i++){
    print OUTFILE "$soluteStuntDoubles[$i]\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

  for ($i = 0; $i<=$#soluteMetaLines; $i++) {
    print OUTFILE "$soluteMetaLines[$i]";
  }
  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:	3114
Committed:	Thu Jan 11 21:46:14 2007 UTC (17 years, 5 months ago) by gezelter
File size:	6887 byte(s)
Log Message:	adding md-solvator
#	User	Rev	Content
1	gezelter	3114	#!/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.1 $"
8			# date = "$Date: 2007-01-11 21:46:14 $"
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 = "";
105			for ($i = 0; $i <= $#line; $i++){
106			$saveLine = join(' ', $saveLine, $line[$i]);
107			}
108			push(@soluteStuntDoubles, $saveLine);
109			$soluteCount++;
110			}
111			}
112			$startSnap = 1 if /Snapshot/;
113			$startFrame = 1 if /FrameData/;
114			$startStunts = 1 if /StuntDoubles/;
115			$startMeta = 1 if /MetaData/;
116			# check again since "/value" contains "value"
117			$startSnap = 0 if /\/Snapshot/;
118			$startFrame = 0 if /\/FrameData/;
119			$startStunts = 0 if /\/StuntDoubles/;
120			$startMeta = 0 if /\/MetaData/;
121			}
122
123			$solventCount = $soluteCount;
124			# okay, let's open the solvent box and write a carved file
125			open(SOLVENT, "$solventName") \|\| die "Error: Can't open file $solventName\n";
126
127
128			while (<SOLVENT>){
129			$startSnap = 0 if /\/Snapshot/;
130			$startFrame = 0 if /\/FrameData/;
131			$startStunts = 0 if /\/StuntDoubles/;
132			$startMeta = 0 if /\/MetaData/;
133
134			# save the MetaData lines
135			push(@metaLines, $_) if $startMeta == 1;
136
137			if ($startSnap == 1){
138			# save the frame data
139			push(@frameData, $_) if $startFrame == 1;
140			if (/Hmat/){
141			@line = split;
142			$hxx = $line[2];
143			chop($hxx);
144			$hyy = $line[8];
145			chop($hyy);
146			$hzz = $line[14];
147			}
148			if ($startStunts == 1){
149			@line = split;
150
151			# wrap positions back into the box
152			$x_val = $line[2] - ($hxx*round($line[2]/$hxx));
153			$y_val = $line[3] - ($hyy*round($line[3]/$hyy));
154			$z_val = $line[4] - ($hzz*round($line[4]/$hzz));
155
156			$saveFlag = 1;
157			for($i=0; $i<=$#solute_x; $i++){
158			$diff_x = $x_val - $solute_x[$i];
159			$diff_y = $y_val - $solute_y[$i];
160			$diff_z = $z_val - $solute_z[$i];
161			$dist2 = $diff_x$diff_x + $diff_y$diff_y + $diff_z*$diff_z;
162			if ($dist2 < $d2tolerance) {
163			$saveFlag = 0;
164			last;
165			}
166			}
167			if ($saveFlag == 1){
168			$saveLine = "$solventCount\t$line[1]\t$line[2]";
169			for ($i = 3; $i <= $#line; $i++){
170			$saveLine = join(' ', $saveLine, $line[$i]);
171			}
172			push(@goodSolventMolecules, $saveLine);
173			$solventCount++;
174			}
175			}
176			}
177			$startSnap = 1 if /Snapshot/;
178			$startFrame = 1 if /FrameData/;
179			$startStunts = 1 if /StuntDoubles/;
180			$startMeta = 1 if /MetaData/;
181			# check again since "/value" contains "value"
182			$startSnap = 0 if /\/Snapshot/;
183			$startFrame = 0 if /\/FrameData/;
184			$startStunts = 0 if /\/StuntDoubles/;
185			$startMeta = 0 if /\/MetaData/;
186			}
187
188			$nSolvent = $#goodSolventMolecules + 1;
189
190			# write out the final file
191			writeOutFile();
192			print "The solvated system \"$fileName\" was generated.\n";
193
194
195			sub writeOutFile {
196			# write out the header
197			print OUTFILE "<OOPSE version=4>\n";
198			printMetaData();
199			printFrameData();
200			print OUTFILE " <StuntDoubles>\n";
201			# now print out the atoms
202			for ($i=0; $i<=$#solute_x; $i++){
203			print OUTFILE "$soluteStuntDoubles[$i]\n";
204			}
205			for ($i=0; $i<=$#goodSolventMolecules; $i++) {
206			print OUTFILE "$goodSolventMolecules[$i]\n";
207			}
208			print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n";
209			}
210
211			sub printMetaData() {
212			print OUTFILE " <MetaData>\n";
213
214			writeSoluteDescription();
215
216			for ($i = 0; $i<=$#metaLines; $i++) {
217			# reset the number of solvent molecules
218			if ($metaLines[$i] =~ /nMol/) {
219			$metaLines[$i] = " nMol = $nSolvent;\n";
220			}
221
222			print OUTFILE "$metaLines[$i]";
223			}
224
225			print OUTFILE " </MetaData>\n";
226			}
227
228			sub writeSoluteDescription {
229			# include a solute model description in the meta data region
230
231			for ($i = 0; $i<=$#soluteMetaLines; $i++) {
232			print OUTFILE "$soluteMetaLines[$i]";
233			}
234			print "The solute model definition was included in '$fileName'.\n";
235			}
236
237			sub printFrameData {
238			print OUTFILE " <Snapshot>\n <FrameData>\n";
239
240			for($i = 0; $i<=$#frameData; $i++) {
241			print OUTFILE "$frameData[$i]";
242			}
243
244			print OUTFILE " </FrameData>\n";
245			}
246
247			sub round {
248			return int( $_[0] + 0.5 * ($_[0] <=> 0) );
249			}