Thank you for your reply.
Because of my workflow, I find myself needing to input fully solvated pdbs (with octahedral boxes) back into LEaP for further simulations and modifications. This step involves manually setting the box with the "set PDB box NUM" command (and changing the angles in the prmtop and inpcrd from square to octahedral). Up to now, I've taken NUM from the restart file I created the pdb from, but this was clunky and always required me to check the restart file. If you tell me there is no way to calculate the box size from the pdb, I will go back to it, but I thought it would be a relatively simple coordinate calculation to output the correct value.
As for the charge calculation, I know pdbs have no charge, but they do have the number of ASP, GLU, LYS and ARG, as well as counterions, so I thought it would be trivial to sum accordingly and give the overall charge. If there is no existing program, I guess I'll just write my own.
Kind regards
Pietro
-----Original Message-----
From: Jason Swails [mailto:jason.swails.gmail.com]
Sent: 22 July 2015 13:23
To: AMBER Mailing List
Subject: Re: [AMBER] Periodic Box Size and Charge Determination
On Wed, Jul 22, 2015 at 7:31 AM, Aronica, Pietro < pietro.aronica07.imperial.ac.uk> wrote:
> Hello,
>
> I have a solvated restart file and its appropriate prmtop. On the last
> line of the restart, there are the three dimensions of the box (all
> three equal, because I used solvateoct) and the three angles. I have
> created a pdb file from the restart and the prmtop using cpptraj
> (center the protein, then autoimage) and then used this<
> http://ambermd.org/tutorials/advanced/tutorial16/#Estimate_Periodic_Bo
> x_Size>
> file to calculate the box dimensions.
After glancing at that script, it seems to me that it is *only* for rectangular boxes. Since you used a truncated octahedron, it does not apply. Basically what that script does is take the waters and draw a rectangular box around it that encompasses all atomic centers. Which runs the implicit assumption that all particles are part of the same unit cell.
This will necessarily give you a periodic box that is at least large enough, but almost always *too* large. Amber images centers of mass of individual molecules, which means that some of the atoms -- like a particular hydrogen -- may actually be protruding from its unit cell and the image that is in the unit cell is on the *other* side of the box. But this would lead to breaking bonds over periodic boundaries and would be a visualization and analysis headache, so it is what it is.
tleap goes a step further when creating its box dimensions. It not only encloses all atom centers in the target box, but adds a buffer region roughly equivalent to the vdW radius of each atom near the boundary. So the resulting boxes are usually *way* too big (relative to the "true" box size that you get after running NpT or the approximation that uses a bounding box like the script above).
The output is very different from the dimensions given at the bottom of the
> restart file. My question is, then, how can I calculate the "correct"
> dimensions from the pdb?
Others may have a better answer than me, but the approach I usually take
is: don't calculate it if you don't have to. The box was already generated for you by tleap -- just use that one. Especially for a non-orthgonal cell that uses familiar (rather than triclinic) imaging, it's much easier just to take the box dimensions that tleap gave you and use NpT to equilibrate the density.
What are the necessary calculations? Also, is there a way that does not
> involve vmd?
All that VMD script is doing is taking the minimum x, y, and z dimensions and subtracting it from the maximum x, y, and z dimensions of all water atoms in a system. You can do this just as easily with your own script.
> Somewhat unrelated, is there an amber functionality that outputs a
> pdb's overall charge? I know there's a leap command, but I tried
> looking for a command-line version and couldn't find one.
>
PDBs don't have charges in the sense that topology files do. You can use ParmEd to calculate the net charge of a prmtop file using the following script (let's call it parmed.in):
parm <your.prmtop>
netCharge :*
Then run it with parmed.py -i parmed.in
You can actually *also* use ParmEd to calculate the same box that the VMD script does, using something like this:
parm <your.prmtop>
netCharge :*
!!
xmin = xmax = ymin = ymax = zmin = zmax = None for atom in amber_prmtop.parm.atoms
if atom.residue.name != 'WAT': continue
if xmin is None:
xmin = xmax = atom.xx
ymin = ymax = atom.xy
zmin = zmax = atom.xz
else:
xmin = min(xmin, atom.xx)
xmax = max(xmax, atom.xx)
ymin = min(ymin, atom.xy)
ymax = max(ymax, atom.xy)
zmin = min(zmin, atom.xz)
zmax = max(zmax, atom.xz)
print('Box dimensions: %s x %s x %s' % (xmax-xmin, ymax-ymin, zmax-zmin)) !!
This time you need the -e flag to have it run the Python code:
parmed.py -i parmed.in -e
But again, this doesn't really apply to a truncated octahedron generated by tleap.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jul 22 2015 - 06:00:05 PDT
