applications/utilities/solvator.in

#!@PERLINTERP@ -w

# 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-03 22:38:50 $"
# 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:	3032
Committed:	Tue Oct 3 22:38:50 2006 UTC (17 years, 9 months ago) by chrisfen
File size:	7757 byte(s)
Log Message:	added solvator and adjusted waterBoxer
#	Content
1	#!@PERLINTERP@ -w
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-03 22:38:50 $"
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	}