[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.3 $"
# date      = "$Date: 2008-01-23 21:22:18 $"
# 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 = "$line[0]\t$line[1]\t$x_val\t$y_val\t$z_val\t";
      for ($i = 5; $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$x_val\t$y_val\t$z_val\t";
        for ($i = 5; $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:	3323
Committed:	Wed Jan 23 21:22:18 2008 UTC (17 years, 9 months ago) by xsun
File size:	6948 byte(s)
Log Message:	added stat2visco
#	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	xsun	3323	# version = "$Revision: 1.3 $"
8			# date = "$Date: 2008-01-23 21:22:18 $"
9	gezelter	3114	# 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	xsun	3323	$saveLine = "$line[0]\t$line[1]\t$x_val\t$y_val\t$z_val\t";
105	gezelter	3117	for ($i = 5; $i <= $#line; $i++){
106	gezelter	3114	$saveLine = join(' ', $saveLine, $line[$i]);
107			}
108			push(@soluteStuntDoubles, $saveLine);
109			$soluteCount++;
110			}
111			}
112	xsun	3323
113	gezelter	3114	$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	xsun	3323	$saveLine = "$solventCount\t$line[1]\t$x_val\t$y_val\t$z_val\t";
170			for ($i = 5; $i <= $#line; $i++){
171	gezelter	3114	$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			}