[ViewVC] Contents of: group/trunk/OOPSE-4/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.2 $"
# date      = "$Date: 2006-10-11 19:24:37 $"
# 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);
$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 "<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:	3043
Committed:	Wed Oct 11 19:24:37 2006 UTC (17 years, 9 months ago) by chrisfen
File size:	7800 byte(s)
Log Message:	now the -n command works
#	Content
1	#!/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.2 $"
8	# date = "$Date: 2006-10-11 19:24:37 $"
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	$soluteName = $opt_n if defined($opt_n);
62
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	print OUTFILE "<OOPSE version=4>\n";
210	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	print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n";
225	}
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	}