AMBER Archive (2009)

Subject: [AMBER] problem in calculation of binding energy in final step (MM-PBSA)

From: Maryam Hamzehee (maryam_h_7860_at_yahoo.com)
Date: Tue Apr 07 2009 - 08:59:03 CDT

Dear Amber users,
I am trying to find the binding energy for my ligand and receptor complex. I had done all the procedures according to the tutorial with minor changes desribed below:
 
I used the following file and script for running the production simulation and obtaining an ensemble of snapshots:
 
 
prod.in file :
 
heat egfr-egf
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=500000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=2, ntp=1, taup=2.0,
  ntpr=5000, ntwx=5000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0,
 /
 
 
Content of run.x script:
 
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod1.out \
-p egf-egfr_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod2.out \
-p egf-egfr_solvated.prmtop -c prod1.rst -r prod2.rst -x prod2.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod3.out \
-p egf-egfr_solvated.prmtop -c prod2.rst -r prod3.rst -x prod3.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod4.out \
-p egf-egfr_solvated.prmtop -c prod3.rst -r prod4.rst -x prod4.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod5.out \
-p egf-egfr_solvated.prmtop -c prod4.rst -r prod5.rst -x prod5.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod6.out \
-p egf-egfr_solvated.prmtop -c prod5.rst -r prod6.rst -x prod6.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod7.out \
-p egf-egfr_solvated.prmtop -c prod6.rst -r prod7.rst -x prod7.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod8.out \
-p egf-egfr_solvated.prmtop -c prod7.rst -r prod8.rst -x prod8.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod9.out \
-p egf-egfr_solvated.prmtop -c prod8.rst -r prod9.rst -x prod9.mdcrd
mpirun -np 32 /opt/pmemd -O -i prod.in -o prod10.out \
-p egf-egfr_solvated.prmtop -c prod9.rst -r prod10.rst -x prod10.mdcrd
gzip -9 prod*.mdcrd
 
 
Subsequently, I extracted snapshots from the production runs for use in the MM-PBSA calculation using extract_coords.mmpbsa file as input:
 
extract_coords.mmpbsa:
 
@GENERAL
PREFIX                snapshot
PATH                        ./
COMPLEX                1
RECEPTOR               1
LIGAND                    1
COMPT                 ./egf-egfr.prmtop
RECPT                  ./egfr.prmtop
LIGPT                   ./egf.prmtop
GC                              1
AS                              0
DC                              0
MM                            0
GB                              0
PB                               0
MS                              0
NM                             0
@MAKECRD
BOX                       YES
NTOTAL                 82201
NSTART                 1
NSTOP                  1000
NFREQ                    1
NUMBER_LIG_GROUPS     1
LSTART                              7852
LSTOP                                 8590
NUMBER_REC_GROUPS     1
RSTART                                     1
RSTOP                                   7851
@TRAJECTORY
TRAJECTORY            ./prod1.mdcrd
TRAJECTORY            ./prod2.mdcrd
TRAJECTORY            ./prod3.mdcrd
TRAJECTORY            ./prod4.mdcrd
TRAJECTORY            ./prod5.mdcrd
TRAJECTORY            ./prod6.mdcrd
TRAJECTORY            ./prod7.mdcrd
TRAJECTORY            ./prod8.mdcrd
TRAJECTORY            ./prod9.mdcrd
TRAJECTORY            ./prod10.mdcrd
@PROGRAMS
 
 
 
In this step I could extract the snapshots successfully using this command:
 
$AMBERHOME/exe/mm_pbsa.pl extract_coords.mmpbsa > extract_coords.log
 
In the final step (i. e. calculation of binding free energy and analyzing the results) I applied this command and input file:
 
$AMBERHOME/exe/mm_pbsa.pl binding_energy.mmpbsa > binding_energy.log
 
binding_energy.mmpbsa content:
@GENERAL
PREFIX                         snapshot
PATH                                ./
COMPLEX                        1
RECEPTOR                       1
LIGAND                            1
COMPT                        ./egf-egfr.prmtop
RECPT                          ./egfr.prmtop
LIGPT                           ./egf.prmtop
GC                                      0
AS                                      0
DC                                      0
MM                                    1
GB                                      1
PB                                       1
MS                                      1
NM                                      0
@PB
PROC                                   2
REFE                                   0
INDI                                  1.0
EXDI                                  80.0
SCALE                                2
LINIT                                  1000
PRBRAD                            1.4
ISTRNG                              0.0
RADIOPT                            0
NPOPT                                1
CAVITY_SURFTEN        0.0072
CAVITY_OFFSET            0.00
SURFTEN                         0.0072
SURFOFF                          0.00
@MM
DIELC                                1.0
@GB
IGB                                      2
GBSA                                   1
SALTCON                          0.00
EXTDIEL                            80.0
INTDIEL                             1.0
SURFTEN                           0.0072
SURFOFF                           0.00
@MS
PROBE                                0.0
@PROGRAMS
 
Unfortunately, calculation procedure stopped suddenly and  pbsa_com.162.out produced and in the binding_energy.log file there are not any information related to the binding energies.
 
Here is the content of this file (i.e. last sentences):
pbsa_com.162.out
total system charges (+/-) for PB        0.0000     1086.9642    -1086.9642
 cavity_surften =        0.0072 cavity_offset =        0.0000
 
  SAS Surface: surface dots generated:    366
 
--------------------------------------------------------------------------------
   3.  ATOMIC COORDINATES AND VELOCITIES
--------------------------------------------------------------------------------
 
 
 begin time read from input coords =     0.000 ps
 
 Number of triangulated 3-point waters found:        0
 
--------------------------------------------------------------------------------
   4.  RESULTS
--------------------------------------------------------------------------------
 
  NB-update: residue-based nb list  4064355
  NB-update: atom-based nb list   883408
 
 
 ======== Setting up Grid Parameters ========
 Using bounding box for grid setup
 Bounding Box Center:      50.500    67.000    18.000
 Xmin, Xmax, Xmax-Xmin:     7.029    94.424    87.395
 Ymin, Ymax, Ymax-Ymin:    13.286   120.412   107.126
 Zmin, Zmax, Zmax-Zmin:   -33.527    70.003   103.530
   beginning box center at level      1     50.500    67.000    18.000
   beginning box center at level      2     50.500    67.000    18.000
 Grid dimension at level     1    45   55   53
 Grid origin corrected at level     1    -41.500   -45.000   -90.000
 Grid dimension at level     2   189  229  221
 Grid origin corrected at level     2      3.000     9.500   -37.500
  SA surface: setting up working radii
  SA surface: found nonzero radii        8590
 Number of SA srf points exposed      110472
 SA Bomb in sa_arc(): Allocation aborted           0           0          41
0                      0
 
 
And the following is related to tail of binding_energy.log
 
Generate PDB
        Center PDB
        Calc PBSA
        Generate PQR
        Calc MS
    Calc contrib for ./snapshot_com.crd.162
        Calc MM/GB/SAS
        Generate PDB
        Center PDB
        Calc PBSA
 
I do not know what happened in this step. I will appreciate it if anyone could kindly help me in this regard.
 
Cheers,
Maryam

      
_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber