Disclaimer: I've never used this method, and don't know anybody,
personally, that actually *has*. My advice is based on how I would do it
based on my knowledge of the method (from reading that paper) and my
knowledge of Amber and MM/PBSA.
On Thu, Apr 16, 2015 at 2:53 PM, Laura Tociu <ltociu.princeton.edu> wrote:
> Dear all,
>
> I am trying to calculate the free energy of binding of a ligand to a
> protein that has around 6000 atoms. Needless to say running NMode on the
> whole protein takes forever, the minimization gets stuck very frequently
> and it also doesn't seem like the drms is going any lower after a while.
>
> I am looking into truncating the trajectory using the method described
> here:
>
> http://www.teokem.lu.se/~ulf/Methods/mm_pbsa.html
>
> I would like to know if the code changes outlined there are still the code
> changes that would be needed in the amber14 code. I am also unsure as to
> why so many confusing steps seem to be needed for this to work.
>
> Wouldn't taking the following steps be enough:
>
> 1) Determine which residues are within 12 angstroms of the ligand (I would
> personally do this by using the VMD command same residue as within 12 of
> ligand) and which are within 8 of the ligand in the same way.
>
> 2) Truncate the trajectory in cpptraj:
>
> parm full_system.parm7
> trajin full_system.nc
> strip : "residues that are not within 12 of ligand" outprefix trunc
> trajout trunc.nc netcdf
>
This produces a "working" prmtop file trunc.full_system.parm7.
> 3) Change the mm_pbsa_createinput file so that it writes input files for
> nmode in which the residues between 8 and 12 Angstroms are restrained.
>
> 4) Remove all the lines in leaprc.ff14SB under addPdbResMap. *--- Is this
> really necessary? *Wouldn't we want the input files for NMode to consist of
> the ligand + the residues within 12 Angstroms of it, capped nicely at their
> loose ends? Or just the residues within 12 Angstroms of the ligand in the
> case of the receptor, but with their loose ends capped nicely? Or is it
> that by truncating the trajectory, the residues are not capped, in which
> case there would be a disagreement between the number of atoms in the
> topology file and that in the input files?
>
I don't think you want to do this. I would suggest just using
trunc.full_system.parm7 generated by cpptraj. Your concerns about atom
mismatching are justified (and even if it doesn't add terminating groups,
it may try to bond two residues that should not be bonded together).You can
then pass this through ante-MMPBSA.py if you needed to to get receptor and
ligand prmtops for your "truncated" system.
> 5) Make nmode able to run thermo analysis even if belly is used by:
>
> In nmode.f (around line 296):
> write(6,*) 'Thermo analysis not supported for belly calc.'
>
> ! Code added by Jacob Kongsted 2007 for new MM/PBSA entropy
> ! JK
> nvecs = 3*natsys - 6
> call thermo(natsys,nvecs,ilevel,x(mx),x(mamass),x(mcval), &
> x(mh),x(mh+ns3),x(mh+2*ns3), &
> x(mh+3*ns3),t,patm)
> ! End of new code
>
> In make_crd_hg.f:
> parameter (NMO=400)
>
>
> Lastly, I am not sure why in the procedure I was referring to, the authors
> were trying to include some water molecules in the buffer zone. Were those
> crystallographic waters? Otherwise I am not sure why there would be waters
> in the trajectory in the first place. Isn't the first thing that MM-PBSA
> does to strip the trajectories of water molecules?
>
PB and GB model bulk solvent. If there are solvent molecules that behave
very much NOT like bulk solvent (i.e., they are structurally important),
then GB and PB treat those waters poorly, and including them explicitly may
be better.
HTH,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Apr 16 2015 - 13:30:02 PDT
