Thanks Jason for the swift reply and obviously for parmed.py...
I done as told, and it kind of works out... I still have one problem, very high vdw and internal energies.
Let me sum up (for the benefit of others) what I did:
I have a trajectory done with explicit solvent, and now I want to extract single residue energies from
the trajectory. This doesn´t seem to work with the explicit solvent since I have to use GB.
So, thanks to Jason, I went this way of preparing my files:
1. create a water and ion stripped .mdcrd and .inpcrd with ptraj
2. create a equally stripped .prmtop using parmed.py (only in AmberTools 12 I figured out ;-)
This works out fine. Now I do the energy calculation using:
#decompose energy on a per-residue basis on the trajectory
&cntrl
imin = 5, igb = 5, gbsa = 2,
ntx = 1, maxcyc = 1,
ntc = 1, ntf = 3,
ntb = 0, ntp = 0,
ntwe = 0, ntpr = 500, ntwx=1,
cut = 1000.0, idecomp = 1,
/
Residues to print
RES 1 227
END
END
What I get is the following. The first frame minimization values look all reasonable:
resid |internal |vdw |eel |pol |sas
============================================================
TDC 1 97.818 -4.235 -250.591 72.307 180.515
TDC 2 71.769 -4.982 -64.649 -3.288 4.889
TDC 3 88.488 -7.834 -50.514 -91.865 60.060
TDC 4 91.894 -5.629 -45.914 -97.714 113.038
TDC 5 61.024 -6.938 -62.792 -4.081 70.127
But in the second frame and following, the internal and vdw jump up to unreasonable high values (the original trajectory itself was stable):
resid |internal |vdw |eel |pol |sas
============================================================
TDC 1 1619.028 1018.080 -82.695 9.640 292.780
TDC 2 1809.027 46.723 -47.603 -2.444 10.724
TDC 3 57146.677 5799.336 21.396 -87.840 36.338
TDC 4 26555.148 ********* 254.986 -146.093 96.970
TDC 5 ********* 44964.405 -18.199 -22.389 83.803
I tried around a little bit with the different igb types (i come to understand for 5 i would have to change PBradii), gbsa and the
cut/rgbmax values.... to no avail, I always get this hight numbers. Might it be the problem that my original calculation used ntb=2/ntp=1?
Or are there other strange things at work here?
Does somebody have suggestions?
Thanks!
André
On 08.10.2012, at 23:40, Jason Swails wrote:
> On Mon, Oct 8, 2012 at 6:52 PM, "André C. Stiel" <
> andre.stiel.tuebingen.mpg.de> wrote:
>
>> Dear Mailing List...
>> I have several problems getting single residue energies out of a system.
>> I tried around for a time with various suggestions from the web/mailing
>> list, to no avail.
>> I have two basic problem:
>> 1. I get no vdw and eel
>> 2. the RRES/ LRES does not work
>> 3. can I do something like GB (i come to understand I maybe need this ;-),
>> but with 3-4 internal waters
>>
>> [snip]
>
>
>
>> My output is:
>>
>> PRINT DECOMP - TOTAL ENERGIES
>>
>> resid |internal |vdw |eel |pol |sas
>> ============================================================
>> TDC 1 19.398 0.000 0.000 0.000 0.000
>> .....
>>
>>> From reading I figured out, that this idecomp does not work properly wit
>> ntb > 0 and I better should
>> strip my molecule of all waters and ions and do igb > 0. Correct?
>>
>
> I would say most likely you are correct. My understanding is that the ONLY
> parts of the code common to electrostatic energy calculations in PME and GB
> is the 1-4 EEL calculation (which is included as part of internal for
> idecomp=1). This suggests that the PME code is not set up for
> decomposition at all (which is not entirely surprising, given that PME
> electrostatic energies are even less pairwise decomposable than GB
> electrostatic and solvation free energies).
>
>
>> So I tried the approach from /test:
>>
>> #decompose energy on a per-residue basis on the trajectory
>> &cntrl
>> imin = 5, igb = 5, gbsa = 2,
>> ntx = 1, maxcyc = 1,
>> ntc = 1, ntf = 3,
>> ntb = 0, ntp = 0,
>> ntwe = 0, ntpr = 500, ntwx=1,
>> cut = 1000.0, idecomp = 1,
>> /
>> Residues considered as REC
>> RRES 1 227
>> END
>> Residues considered as LIG
>> LRES 1 227
>> END
>> Residues to print
>> RES 1 227
>> END
>> END
>>
>> With a stripped molecule and .prmtop I get:
>>
>> FATAL: NATOM mismatch in coord and topology files
>
>
>> Is there a trick to go round this problem. I tried to produce the new -
>> water free - .prmtop with tleap on the original .pdb
>> But obviously i doesnt work... but ok, I will figure out this.
>>
>
> Try using ParmEd to strip out the waters from your topology file. It
> should be able to write out a new topology file without ever having to go
> through LEaP. Make sure you use the *exact* same strip mask for ParmEd
> that you used for ptraj to generate the 'dry' trajectories.
>
>
>> Second and Third Problem can be easily explained on this basis. 2. Even
>> with the RRES/LRES/RES from the /test file, I got:
>>
>> LOADING THE DECOMP ATOMS AS GROUPS
>>
>> ----- READING GROUP 1; TITLE:
>> Residues considered as REC
>> Number of atoms in this group = 0
>> ----- READING GROUP 2; TITLE:
>> Residues considered as LIG
>> Number of atoms in this group = 0
>> ----- READING GROUP 3; TITLE:
>> Residues to print
>> GRP 3 RES 1 TO 227
>> Number of atoms in this group = 3541
>>
>> So LRES and RRES don´t work. Why is that.
>>
>
> I'm not sure. To be honest, I've never understood what this distinction
> meant in the decomp code, or what it was actually used for.
>
>
>> Question 3. Some internal waters are important I guess, can I combine igb
>> and a stripped molecule with some internal waters?
>>
>
> Yes. Of course, I must include the obligatory comments about how GB was
> never intended to deal with actual water, but given the fact that GB decomp
> energies are not *really* pairwise decomposable and the results are at best
> useful as a qualitative measure for your system, this doesn't seem to be a
> serious limitation.
>
> HTH,
> Jason
>
> --
> 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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Oct 09 2012 - 10:00:03 PDT
