Re: [AMBER] MMPBSA.py error

From: Jason Swails <jason.swails.gmail.com>
Date: Tue, 26 Apr 2011 19:46:48 -0700

Hi Gabriel,

You should continue sending these emails to the mailing list so that others
may be able to assist and they're archived so others can benefit from the
answers.


On Tue, Apr 26, 2011 at 2:56 PM, Gabriel Urquiza <urquizagabes.gmail.com>wrote:

> Dear Jason
>
> I have tried as you told me and found out that there were a couple of
> crystallization waters along with the
> receptor, so basically you were right and the problem was with the strip
> mask. However, when I run the
> MMPBSA.py to make GBSA calculations using several trajectories, it makes
> averages for a single frame
> for each trajectory (a total of 10 frames for 10 trajectories).
>

The way MMPBSA.py works is that it takes the frames that correspond to
"startframe" "endframe" and "interval" from each trajectory -- it doesn't
combine them first (this will probably be changed in future versions, but
that's the way it is now). Thus, if you set startframe=50,
endframe=1000000, interval=50, and you specify 100 trajectories each with 90
frames, then it will only grab 1 snapshot (snapshot #50) from each
trajectory file.

HTH,
Jason


>
> Do you have any clue as to why is this so?
>
> Regards
>
> Gabriel
>
> P.S. I'm using quasi-harmonic approximation to calculate entropic
> contributions to free energy.
>
>
> 2011/4/20 Jason Swails <jason.swails.gmail.com>
>
>> Hi Gabriel,
>>
>> You should continue this thread on the mailing list so it's archived and
>> people can benefit from them in the future.
>>
>> On Wed, Apr 20, 2011 at 10:22 AM, Gabriel Urquiza <urquizagabes.gmail.com
>> > wrote:
>>
>>> Jason,
>>>
>>> I have done as you said and found out the topologies were kind of
>>> problematic. However, in tracking down the
>>> error I found out that leap is miswriting the topologies even though the
>>> files from which it writes are right. I can't
>>> find out what is going on. On the one hand, for instance, the ligand mol2
>>> is perfect. Then, after the topology has been
>>> written down, I check the results on VMD and it is hideous.
>>>
>>
>> Be careful here if you're using ASCII trajectories. There are two options
>> for loading MDCRD files in VMD, and choosing the wrong one can give you this
>> effect.
>>
>> Make sure that you're loading properly either "with periodic box" or not.
>> Note that if you run a solvated trajectory, you'll have to choose the one
>> "with box info". MMPBSA.py automatically strips the box information from
>> the trajectories, so in that case you DON'T want to load that kind of
>> trajectory.
>>
>> An alternative to all of this is to just use netcdf trajectories (ioutfm=1
>> in your sander/pmemd input file), which VMD reads the same way whether or
>> not there is box information there.
>>
>>
>>>
>>> There is also one more puzzling event. When I use ambpdb to generate a
>>> PDB from the wrecked topolgies, the pdb comes
>>> out perfect, as if it had been written from good topologies. How can the
>>> very same pair topology-coordinates generate junk on
>>> vmd but also a perfect pdb? Is it common? Does the topology scheme has
>>> changed in Amber11? Perhaps VMD isn't showing
>>> data accurately? Is there any other visualizer for Amber11 topologies and
>>> coordnates?
>>>
>>
>> I've never had problems with VMD/pmemd. My guess is that your topology is
>> not really "wrecked". Keep in mind that if you dump a PDB from the first
>> frame using a topology file and a trajectory file that has coordinates for
>> more atoms than you have in your prmtop file, your first frame will look
>> *fine*, since all of the atoms are present and in the correct place. After
>> the first frame, however, it will look awful, as the first atom of the
>> second frame will take the coordinates of the first atom after the last atom
>> present in the topology file (from the FIRST frame).
>>
>>
>>>
>>> I am more puzzled now than when I first started, but at least the
>>> checkCoordinates() problem was clarified =P. Do you have any clue
>>> on what's happening now? And by the way, is it even possible to run ten
>>> to twenty molecular dynamics trajectories with wrecked topologies?
>>> How come pmemd didn't crashed even before I generated the trajectories
>>> used in the mmgbsa calculations? Sorry to bomb you with all these
>>> doubts.
>>>
>>
>> The only reason it wouldn't have crashed is because the topology file was
>> fine. If you used the topology file and coordinate file printed out by the
>> same Leap session, then your topology file is guaranteed* to be consistent
>> with your inpcrd file (*barring bugs, but at least the atom count between
>> them will match).
>>
>> My guess is that you're having a problem with ions in your MM/PBSA
>> calculation. You adjusted the strip_mask so that it only gets rid of WAT
>> and Na+, yet are there other ions in your prmtop file? If there are, for
>> instance, Cl- ions in your solvated topology file that you took out when you
>> created your complex/receptor/ligand topology files, then you will run into
>> this problem. Alternatively, if you kept Na+ ions in your
>> complex/receptor/ligand topology files, then MMPBSA will have stripped them
>> out, your prmtops/trajectories won't match, and you'll also see this
>> problem.
>>
>> HTH,
>> Jason
>>
>>
>>>
>>> Regards
>>>
>>> Gabriel
>>>
>>> 2011/4/20 Jason Swails <jason.swails.gmail.com>
>>>
>>>> On Tue, Apr 19, 2011 at 7:34 PM, Gabriel Urquiza <
>>>> urquizagabes.gmail.com>wrote:
>>>>
>>>> > Dear amberists,
>>>> >
>>>> > I am having some problemas when I try to run MMPBSA.py on several
>>>> short
>>>> > dynamics. My system is an explicitly solvated protein with a
>>>> > ligand, sodium ions and a flavin residue (which has been parametrized
>>>> using
>>>> > antechamber and parmchk).
>>>> >
>>>> > I run it with this input:
>>>> >
>>>> > &general
>>>> > startframe=1,
>>>> > endframe=1000000,
>>>> > interval=50,
>>>> > verbose=2,
>>>> > keep_files=1,
>>>> > strip_mask=":WAT:Na+",
>>>> >
>>>>
>>>> You should leave this at its default value, unless you want to keep any
>>>> ions
>>>> in there.
>>>>
>>>>
>>>> > ligand_mask=":O01",
>>>> >
>>>>
>>>> If you don't specify both receptor_mask and ligand_mask, it will ignore
>>>> this
>>>> and just use the default.
>>>>
>>>> entropy=1,
>>>> > /
>>>> >
>>>> > &gb
>>>> > igb=5,
>>>> > saltcon=0.150,
>>>> > /
>>>> >
>>>> > &pb
>>>> > istrng=0.15, fillratio=4.0
>>>> > /
>>>> >
>>>> > The residue name is O01 because it is the first of a series of related
>>>> > compounds. I run MMPBSA.py and it gives me this:
>>>> >
>>>> > checkCoordinates(): Could not predict number of frames for AMBER
>>>> trajectory
>>>> > file: _MMPBSA_complex.mdcrd
>>>> > If this is not a compressed file then there is a problem
>>>> >
>>>>
>>>> This suggests that your prmtop files don't match the trajectory files
>>>> that
>>>> are generated by MMPBSA.py.
>>>>
>>>>
>>>> > checkCoordinates(): Could not predict number of frames for AMBER
>>>> trajectory
>>>> > file: _MMPBSA_complex.mdcrd
>>>> > If this is not a compressed file then there is a problem
>>>> > checkCoordinates(): Could not predict number of frames for AMBER
>>>> trajectory
>>>> > file: _MMPBSA_complex.mdcrd
>>>> > If this is not a compressed file then there is a problem
>>>> > checkCoordinates(): Could not predict number of frames for AMBER
>>>> trajectory
>>>> > file: _MMPBSA_complex.mdcrd
>>>> > If this is not a compressed file then there is a problem
>>>> > Error: Sander output is missing values!
>>>> > VDWAALS = ************* EEL = -22412.2843 EGB =
>>>> > -3345.9563
>>>> >
>>>>
>>>> This all but confirms it. Try visualizing _MMPBSA_complex.mdcrd with
>>>> your
>>>> complex prmtop, _MMPBSA_receptor.mdcrd with your receptor prmtop, and
>>>> _MMPBSA_ligand.mdcrd with your ligand prmtop. If the structures look
>>>> warped, then something is messed up and you have to fix your prmtops (or
>>>> just use the default strip_mask).
>>>>
>>>> HTH,
>>>> Jason
>>>>
>>>>
>>>> >
>>>> > ptraj found! Using /opt/amber11/exe/ptraj
>>>> > sander found! Using /opt/amber11/exe/sander (serial only!)
>>>> >
>>>> > Preparing trajectories with ptraj...
>>>> > Beginning quasi-harmonic entropy calculation with ptraj...
>>>> > 70 frames were read in and processed by ptraj for use in calculation.
>>>> >
>>>> > Starting sander calls
>>>> >
>>>> > Starting gb calculation...
>>>> >
>>>> > Starting pb calculation...
>>>> >
>>>> > NOTE: All files have been retained for debugging purposes. Type
>>>> MMPBSA.py
>>>> > --clean to erase these files.
>>>> >
>>>> > There are no values of GB or PB in the output. Something probably went
>>>> > wrong
>>>> > (no kidding?).
>>>> >
>>>> > I made PDB files for the minimization step and for several snapshots
>>>> of the
>>>> > dynamics using ambpdb. I also checked THE dynamics
>>>> > themselves. All of that using VMD. The snapshots seemed fine as well
>>>> as the
>>>> > trajectories. I also checked the topologies used by the
>>>> > MMPBSA.py and all of them were at the right place.
>>>> >
>>>> > I looked around the list for problems such as these but i didn't find
>>>> one
>>>> > that matched the context I am facing right now. Also, my trajectories
>>>> > aren't compressed.
>>>> >
>>>> > Looking forward to hear your ideas.
>>>> >
>>>> > Gabriel
>>>> >
>>>> > -----
>>>> > Graduation student at Federal University of Paraíba, Brazil.
>>>> > LQQC - Laboratory of Computational Quantum Chemistry
>>>> > _______________________________________________
>>>> > AMBER mailing list
>>>> > AMBER.ambermd.org
>>>> > http://lists.ambermd.org/mailman/listinfo/amber
>>>> >
>>>>
>>>>
>>>>
>>>> --
>>>> Jason M. Swails
>>>> Quantum Theory Project,
>>>> University of Florida
>>>> Ph.D. Candidate
>>>> 352-392-4032
>>>> _______________________________________________
>>>> AMBER mailing list
>>>> AMBER.ambermd.org
>>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>>
>>>
>>>
>>
>>
>> --
>> Jason M. Swails
>> Quantum Theory Project,
>> University of Florida
>> Ph.D. Candidate
>> 352-392-4032
>>
>
>


-- 
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Apr 26 2011 - 20:00:03 PDT
Custom Search