[re: automatically building prmtops to have a set number of waters]
> I would like to solvate two different peptide conformations with truncated
> octahedral boxes of water that are of the same dimensions. Can anyone
> suggest a good way to specify the dimensions of an isotropic truncated
> octahedron?
I think what would be better is to have the same number of waters (rather
than the same box size) since then you could compare the energies (and
having the same box size initially does not mean post constant pressure
equilibration that the box sizes will be equivalent anyways). If I really
want the box sizes equivalent, what I would do is equilibrate both runs
and then manually edit the box size of the smaller to match the larger,
but again it makes much more sense to have an equivalent number of waters
rather than box size.
solvateOct only takes a single parameter for the buffer as all the box
lengths are equivalent. To get LEaP to produce the same number of waters
(or the same box size) requires trial and error with the buffer...
Below is a perl script I hacked together to get the same number of waters,
in this case 4000 waters. If you need a larger initial buffer to get it
to fly, increase the "size". What this does is sequentially run LEaP
increasing or decreasing the buffer size (and eventually the closeness if
necessary) until a topology with 4000 waters is built. This could be
modified to do the same for box size (or you can do it by hand by fooling
with the "closeness").
*
solvateoct aa TIP3PWAT 10.0 iso 1.0
The script below reads in an initial pdb named "initial.pdb", solvates,
adds net-neutralizing Na+ and then 15 extra Na+ and Cl-. I built this
when trying to build a lot of inputs and got tired of doing it by hand...
--tom
#!/usr/bin/perl
$totalwaters = 4000;
$finished = 0;
$size = 15.0;
$waters = 0;
$increment = 1.0;
$closeness = 1.0;
$iteration = 0;
$prev= 0;
while ( ! $finished || $iteration > 250 ) {
open(OUT, "> buildit");
print OUT "\naa = loadpdb initial.pdb\n";
printf(OUT "\nsolvateoct aa TIP3PBOX %f iso %f\n\n", $size, $closeness);
print OUT "addions aa Na+ 0\n\n";
for ($i=0; $i < 15; $i++) {
print OUT "addions aa Na+ 1 Cl- 1\n";
}
print OUT "\nsaveamberparm aa solvated.topo solvated.coords\n";
print OUT "\nquit\n";
close(OUT);
open(INP, "tleap -f leaprc.ff99sb -f buildit |");
while ( <INP> ) {
# print $_;
$input = $_;
if ( /WAT/ ) {
# print "INPUT: ", $_;
($junk, $waters) = split(' ', $_);
printf("*** Size is %f, solvated with %i waters; increment is %f\n", $size, $waters, $increment);
}
}
$iteration++;
$prev = $cur;
$cur = $waters;
if ($prev != 0 && $cur < $totalwaters && $prev > $totalwaters) {
$increment = $increment/3.0;
$prevsize = $size;
$size = $size + $increment;
} elsif ($prev != 0 && $cur > $totalwaters && $prev < $totalwaters ) {
$increment = $increment/3.0;
$prevsize = $size;
$size = $prevsize - $increment;
} elsif ($cur < $totalwaters) {
$prevsize = $size;
$size = $size + $increment;
} elsif ($cur > $totalwaters) {
$prevsize = $size;
$size = $size - $increment;
}
if ($increment < 0.000001) {
$increment = 0.04;
$closeness = $closeness * 1.02;
printf(" Bumping closeness to %f\n", $closeness);
}
if ($cur == $totalwaters) {
$finished = 1;
}
}
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Aug 05 2007 - 06:08:13 PDT
