Dear Tom,
Thank you very much for the helpful response and providing your useful
script. It does make better sense to have the same number of waters rather
than the same box size -- thanks for suggesting that, too.
All the best,
Lillian
On 8/4/07, Thomas Cheatham <tec3.utah.edu> wrote:
>
>
> [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
>
--
Lillian T. Chong
Assistant Professor
University of Pittsburgh
Department of Chemistry
219 Parkman Avenue
Pittsburgh, PA 15260
(412) 624-6026
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Aug 08 2007 - 06:07:06 PDT
