Dave beat me to the response, but I'll send it anyway since it was already
done... :)
On Mon, May 13, 2013 at 6:55 AM, <dilrajl.bii.a-star.edu.sg> wrote:
> Dear Amber Users,
> I am trying to calculate the free energy of a
> protein-peptide system using mmpbsa.py script from Amber12.
>
> My input script file for the purpose is the following.
>
> sample input file for MM-GBSA calculation and DECOMPOSITION
> &general
> startframe=1, endframe=50000, interval=10,
> verbose=2, keep_files=2,
> /
> &gb
> igb=2, saltcon=0.100,
> /
> &decomp
> idecomp=2,
> dec_verbose=0,
> csv_format=0,
> /
>
> The calculation goes on fine. The final statistics is as follows.
>
> |MMPBSA.py Version=12.0
> |Complex topology file: ../../complex.parmtop
> |Receptor topology file: ../../protein.parmtop
> |Ligand topology file: ../../peptide.parmtop
> |Initial mdcrd(s):
> ../../../centering/prod_cen_strip_sol_ion.mdcrd
> |
> |Receptor mask: ":1-189"
> |Ligand mask: ":190-203"
> |
> |Calculations performed using 5000 complex frames.
> |
> |All units are reported in kcal/mole.
>
> -------------------------------------------------------------------------------
>
> -------------------------------------------------------------------------------
>
> GENERALIZED BORN:
>
> Complex:
> Energy Component Average Std. Dev. Std. Err. of
> Mean
>
> -------------------------------------------------------------------------------
> BOND 647.4911 21.6041
> 0.3055
> ANGLE 1682.8402 31.8009
> 0.4497
> DIHED 2132.8508 18.9976
> 0.2687
> VDWAALS 4539.8345 1530.7455
> 21.6480
> EEL -13699.1924 108.8936
> 1.5400
> 1-4 VDW 752.4023 11.6313
> 0.1645
> 1-4 EEL 7589.7396 52.8488
> 0.7474
> EGB -3163.9703 88.1210
> 1.2462
> ESURF 84.4256 1.7838
> 0.0252
>
> G gas 3645.9662 1532.2705
> 21.6696
> G solv -3079.5447 87.5957
> 1.2388
>
> TOTAL 566.4215 1530.7715
> 21.6484
>
>
> Receptor:
> Energy Component Average Std. Dev. Std. Err. of
> Mean
>
> -------------------------------------------------------------------------------
> BOND 599.0144 20.8510
> 0.2949
> ANGLE 1564.8283 30.4181
> 0.4302
> DIHED 1985.2735 18.0507
> 0.2553
> VDWAALS -1475.0583 18.4672
> 0.2612
> EEL -12804.7448 99.3452
> 1.4050
> 1-4 VDW 694.2489 11.2516
> 0.1591
> 1-4 EEL 7385.6220 52.5821
> 0.7436
> EGB -3188.6682 89.0650
> 1.2596
> ESURF 79.9245 1.4009
> 0.0198
>
> G gas -2050.8160 103.4176
> 1.4625
> G solv -3108.7437 88.4805
> 1.2513
>
> TOTAL -5159.5597 42.0311
> 0.5944
>
>
> Ligand:
> Energy Component Average Std. Dev. Std. Err. of
> Mean
>
> -------------------------------------------------------------------------------
> BOND 48.4767 5.9530
> 0.0842
> ANGLE 118.0119 8.4970
> 0.1202
> DIHED 147.5773 5.1734
> 0.0732
> VDWAALS 6069.3243 1530.2449
> 21.6409
> EEL -373.5075 16.7907
> 0.2375
> 1-4 VDW 58.1534 3.0488
> 0.0431
> 1-4 EEL 204.1176 10.3060
> 0.1457
> EGB -506.2401 12.3892
> 0.1752
> ESURF 14.2383 0.4659
> 0.0066
>
> G gas 6272.1537 1530.1827
> 21.6401
> G solv -492.0018 12.4146
> 0.1756
>
> TOTAL 5780.1519 1530.2265
> 21.6407
>
>
> Differences (Complex - Receptor - Ligand):
> Energy Component Average Std. Dev. Std. Err. of
> Mean
>
> -------------------------------------------------------------------------------
> BOND 0.0000 0.0000
> 0.0000
> ANGLE -0.0000 0.0000
> 0.0000
> DIHED -0.0000 0.0000
> 0.0000
> VDWAALS -54.4314 5.9129
> 0.0836
> EEL -520.9400 67.0862
> 0.9487
> 1-4 VDW 0.0000 0.0000
> 0.0000
> 1-4 EEL -0.0000 0.0001
> 0.0000
> EGB 530.9380 61.6686
> 0.8721
> ESURF -9.7372 0.8343
> 0.0118
>
> DELTA G gas -575.3715 67.4110
> 0.9533
> DELTA G solv 521.2008 61.2946
> 0.8668
>
> DELTA TOTAL -54.1707 8.7948
> 0.1244
>
>
>
> -------------------------------------------------------------------------------
>
> -------------------------------------------------------------------------------
>
>
> The "DELTA TOTAL" value for the complex and ligand doesn't look
> correct to me. The energy values are positive and the standard
> deviations are huge. The same values for the Receptors and finally for
> the Differences (Complex - Receptor - Ligand) looks fine.
>
> Is there anything wrong with the calculation or is there something I
> am missing??
>
This is where looking at the energy breakdown and understanding how to run
MM/PBSA calculations 'by hand' help (see the introduction to the MMPBSA.py
chapter in the AmberTools manual).
The key interaction that is giving rise to this behavior is the van der
Waals term. For large systems, this term is typically (but not
necessarily!) negative. This is clearly true for the complex and ligand
ensembles. The large standard deviation suggests to me that there are a
couple frames that you've analyzed that have partial atomic overlaps,
giving rise to such a large energy. You should visualize these
trajectories (_MMPBSA_ligand.mdcrd and _MMPBSA_complex.mdcrd) with the
complex and ligand topology files to see if you can find the offending
frames. The checkoverlaps command in cpptraj can also help you here (see
the manual for more instructions).
Also, the fact that the total values are OK suggest even more---the vdW
clashes are localized in the ligand. Because you are using a
single-trajectory approach for MM/PBSA binding, energy terms that contain
contributions just from the receptor or just from the ligand cancel
completely (since they are present equally in the bound and unbound
states). This includes all internal bonded terms as well as non-bonded
interactions from pairwise-decomposable terms (EEL and VDW, but none of the
solvation terms) between 2 receptor or 2 ligand atoms. So you can see why
my suspicion is that the problem is in the ligand only, and therefore
cancels in the difference.
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
Received on Mon May 13 2013 - 05:30:02 PDT
