Hi Jason,
I followed your suggestion and set radiopt=0 and inp=1 in your &pb section and try again… but there are warnings when running MMPBSA.py regarding the inp=1 which I don’t understand… (please see below). The energies are much more similar to the tutorial.
Are then those parameters you suggested the correct ones despite teh warnings?
All the details below.
Thanks lot,
Fabian
Here is the input file modified as you suggested
Input file for running PB and GB
&general
verbose=1,
entropy=0,
/
&gb
igb=2, saltcon=0.100
/
&pb
istrng=0.100,
radiopt=0, inp=1
/
And here is the output and the warnings I get when I run it with my own mdcrd file (10 ns in this case):
==> MMPBSA.py -O -i mmpbsa_2b.in -o FINAL_RESULTS_MMPBSA.dat -sp 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y *.mdcrd
Loading and checking parameter files for compatibility...
mmpbsa_py_energy found! Using /Users/admin/Data/TOOLS/amber14/bin/mmpbsa_py_energy
cpptraj found! Using /Users/admin/Data/TOOLS/amber14/bin/cpptraj
Preparing trajectories for simulation...
50 frames were processed by cpptraj for use in calculation.
Running calculations on normal system...
Beginning GB calculations with /Users/admin/Data/TOOLS/amber14/bin/mmpbsa_py_energy
calculating complex contribution...
calculating receptor contribution...
calculating ligand contribution...
Beginning PB calculations with /Users/admin/Data/TOOLS/amber14/bin/mmpbsa_py_energy
calculating complex contribution...
Warning: inp=1 was old default
Warning: cavity_surften=0.0378 not recommended for inp=1, switching to inp=1 default value: 0.0050
Warning: cavity_offset=-0.5692 not recommended for inp=1, switching to inp=1 default value: 0.000
Warning: sprob=.557 not recommended for inp=1, switching to inp=1 default value: 1.400
calculating receptor contribution...
Warning: inp=1 was old default
Warning: cavity_surften=0.0378 not recommended for inp=1, switching to inp=1 default value: 0.0050
Warning: cavity_offset=-0.5692 not recommended for inp=1, switching to inp=1 default value: 0.000
Warning: sprob=.557 not recommended for inp=1, switching to inp=1 default value: 1.400
calculating ligand contribution...
Warning: inp=1 was old default
Warning: cavity_surften=0.0378 not recommended for inp=1, switching to inp=1 default value: 0.0050
Warning: cavity_offset=-0.5692 not recommended for inp=1, switching to inp=1 default value: 0.000
Warning: sprob=.557 not recommended for inp=1, switching to inp=1 default value: 1.400
Timing:
Total setup time: 0.035 min.
Creating trajectories with cpptraj: 0.293 min.
Total calculation time: 16.376 min.
Total GB calculation time: 2.464 min.
Total PB calculation time: 13.912 min.
Statistics calculation & output writing: 0.000 min.
Total time taken: 16.706 min.
MMPBSA.py Finished! Thank you for using. Please cite us if you publish this work with this paper:
Miller III, B. R., McGee Jr., T. D., Swails, J. M. Homeyer, N. Gohlke, H. and Roitberg, A. E.
J. Chem. Theory Comput., 2012, 8 (9) pp 3314–3321
But the energies are correct this time.
Here is the FINAL_RESULTS_MMPBSA.dat file
| Run on Tue Oct 20 13:11:30 2015
|
|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB in serial
|&general
| endframe=50, keep_files=2,
|/
|&gb
| igb=2, saltcon=0.100,
|/
|&pb
| istrng=0.100,
| radiopt=0, inp=1
|/
|--------------------------------------------------------------
|MMPBSA.py Version=14.0
|Solvated complex topology file: 1err.solvated.prmtop
|Complex topology file: complex.prmtop
|Receptor topology file: receptor.prmtop
|Ligand topology file: ligand.prmtop
|Initial mdcrd(s): 1err_prod.mdcrd
|
|Receptor mask: ":1-240"
|Ligand mask: ":241"
|Ligand residue name is "RAL"
|
|Calculations performed using 50 complex frames.
|Poisson Boltzmann calculations performed using internal PBSA solver in mmpbsa_py_energy
|
|Generalized Born ESURF calculated using 'LCPO' surface areas
|
|All units are reported in kcal/mole.
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
GENERALIZED BORN:
Complex:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -2013.3861 20.1108 2.8441
EEL -16938.6434 84.8707 12.0025
EGB -3506.9860 67.1003 9.4894
ESURF 97.4744 1.3170 0.1863
G gas -18952.0295 81.3542 11.5052
G solv -3409.5116 66.5951 9.4180
TOTAL -22361.5411 36.3100 5.1350
Receptor:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1955.2346 19.0529 2.6945
EEL -16895.0333 84.6812 11.9757
EGB -3528.7058 67.6664 9.5695
ESURF 101.2621 1.2926 0.1828
G gas -18850.2679 80.7202 11.4156
G solv -3427.4437 67.2065 9.5044
TOTAL -22277.7116 34.9512 4.9428
Ligand:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1.8595 2.0320 0.2874
EEL -5.5800 2.0133 0.2847
EGB -28.4852 0.5979 0.0846
ESURF 4.4805 0.0425 0.0060
G gas -7.4395 3.0299 0.4285
G solv -24.0047 0.6040 0.0854
TOTAL -31.4442 3.2513 0.4598
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -56.2921 2.8956 0.4095
EEL -38.0301 3.1822 0.4500
EGB 50.2050 2.5622 0.3624
ESURF -8.2682 0.2700 0.0382
DELTA G gas -94.3222 3.3476 0.4734
DELTA G solv 41.9369 2.4569 0.3475
DELTA TOTAL -52.3853 2.4400 0.3451
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
POISSON BOLTZMANN:
Complex:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -2013.3861 20.1108 2.8441
EEL -16938.6434 84.8707 12.0025
EPB -3314.8738 66.2632 9.3710
ENPOLAR 63.9049 0.4744 0.0671
G gas -18952.0295 81.3542 11.5052
G solv -3250.9688 66.1312 9.3524
TOTAL -22202.9984 38.5993 5.4588
Receptor:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1955.2346 19.0529 2.6945
EEL -16895.0333 84.6812 11.9757
EPB -3340.9508 66.5434 9.4107
ENPOLAR 65.6147 0.4833 0.0684
G gas -18850.2679 80.7202 11.4156
G solv -3275.3361 66.4327 9.3950
TOTAL -22125.6040 37.0760 5.2433
Ligand:
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -1.8595 2.0320 0.2874
EEL -5.5800 2.0133 0.2847
EPB -30.7649 0.6756 0.0955
ENPOLAR 3.8717 0.0261 0.0037
G gas -7.4395 3.0299 0.4285
G solv -26.8932 0.6809 0.0963
TOTAL -34.3327 3.2575 0.4607
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS -56.2921 2.8956 0.4095
EEL -38.0301 3.1822 0.4500
EPB 56.8420 2.9272 0.4140
ENPOLAR -5.5814 0.1540 0.0218
EDISPER 0.0000 0.0000 0.0000
DELTA G gas -94.3222 3.3476 0.4734
DELTA G solv 51.2605 2.9326 0.4147
DELTA TOTAL -43.0617 2.8642 0.4051
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
____________________
Dr. Fabian Glaser
Head of the Structural Bioinformatics section
Bioinformatics Knowledge Unit - BKU
The Lorry I. Lokey Interdisciplinary Center for Life Sciences and Engineering
Technion - Israel Institute of Technology, Haifa 32000, ISRAEL
fglaser at technion dot ac dot il
Tel: +972 4 8293701
http://bku.technion.ac.il
> On 19 Oct 2015, at 6:24 PM, Jason Swails <jason.swails.gmail.com> wrote:
>
> That's because in the time between when the tutorial was written and the
> current version of the code, the default nonpolar solvation model changed,
> unfortunately. If you compare the *polar* solvation free energy component
> (EPB), you *should* find that they are very similar.
>
> Indeed, if you compare them, the tutorial lists 57.6402 kcal/mol as the
> polar solvation free energy of binding, while your output gives 56.7534
> kcal/mol. I don't know why the discrepancy is as large as 1.1 kcal/mol,
> but it's certainly a lot smaller than 35 kcal/mol.
>
> If you want to recover the original tutorial results (or values close to
> it), set radiopt=0 and inp=1 in your &pb section and try again. I really
> should try to find some time to update the tutorial...
>
> Hope this helps,
> Jason
>
> On Mon, Oct 19, 2015 at 10:17 AM, Fabian gmail <fabian.glaser.gmail.com <mailto:fabian.glaser.gmail.com>>
> wrote:
>
>> Dear Amber colleagues,
>>
>> I am doing all the MMPBSA tutorials, and I found consistently that while
>> the GENERALISED BORN results are almost identical to the tutorial (using
>> the tutorial .mdrcd and .prmtop files) the POISSON BOLTZMANN results are
>> completely different. As I said I only did the MMPBSA.py calculation, using
>> given files. For example using the following command for tutorial files A
>> 3.2 (
>> http://ambermd.org/tutorials/advanced/tutorial3/py_script/section2.htm <http://ambermd.org/tutorials/advanced/tutorial3/py_script/section2.htm> <
>> http://ambermd.org/tutorials/advanced/tutorial3/py_script/section2.htm <http://ambermd.org/tutorials/advanced/tutorial3/py_script/section2.htm>>)
>> I got the following FINAL_RESULTS_MMPBSA.dat file, which as you see for
>> the POISSON BOLTZMANN yields a much lower DELTA TOTAL -6.6887 compared to
>> the tutorial DELTA G binding = -41.7256 !!
>>
>> Why is that?????
>>
>> As explained above the GENERALISED BORN results are much more similar, for
>> my calculation is DELTA TOTAL -52.3853 and the tutorial is
>> DELTA G binding = -52.2672. That is great.
>>
>> I used the command and input files provided on the tutorial:
>>
>> MMPBSA.py -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp
>> 1err.solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp
>> ligand.prmtop -y *.mdcrd
>>
>> I want to test this method on several systems, so it’s important for me to
>> understand the differences…. any hint will really help me. The same happens
>> on 3.1 so I suspect there is something different or wrong with the
>> calculation.
>>
>> Thanks a lot in advance,
>>
>> Fabian
>>
>> ____________________
>> Dr. Fabian Glaser
>>
>> Head of the Structural Bioinformatics section
>>
>> Bioinformatics Knowledge Unit - BKU
>> The Lorry I. Lokey Interdisciplinary Center for Life Sciences and
>> Engineering
>> Technion - Israel Institute of Technology, Haifa 32000, ISRAEL
>>
>> fglaser at technion dot ac dot il
>> Tel: +972 4 8293701
>> http://bku.technion.ac.il
>>
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org <mailto:AMBER.ambermd.org>
>> http://lists.ambermd.org/mailman/listinfo/amber <http://lists.ambermd.org/mailman/listinfo/amber>
>>
>
>
>
> --
> Jason M. Swails
> BioMaPS,
> Rutgers University
> Postdoctoral Researcher
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org <mailto:AMBER.ambermd.org>
> http://lists.ambermd.org/mailman/listinfo/amber <http://lists.ambermd.org/mailman/listinfo/amber>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Oct 20 2015 - 03:30:05 PDT
