AMBER Archive (2002)

Subject: Energies/coordinate files from mm_pbsa problem

From: Peter Gannett (pgannett_at_hsc.wvu.edu)
Date: Tue Jul 30 2002 - 13:25:29 CDT

Dear Amber users:

We are analyzing some trajectories with mm_pbsa and are getting some
strange numbers. The trajectories are sufficiently large that we have
to use mm_pbsa to sample the trajectory every X number of frames. For
example, we started with frame 500 or 501 or 502 ... and sample every 14
frames. Regardless of the starting point we get the same set of
energies. Not surprisingly, when we use rdparm to generate some
coordinate files (based on the starting frame) they are all they same,
as though regardless of where we tell it to start (500, 501...) it is
starting at the same point. If we start near the beginning or end of
the trajectory, this doesn't happen.

Below is the control file we are using and then samples of the energy
output. Any suggestions for what we are screwing up would be most
appreciated.

Pete Gannett

Here is the control file:
#
# Input parameters for mm_pbsa.pl
#
# Holger Gohlke
# 08.01.2002
#
################################################################################
@GENERAL
#
# General parameters
# 0: means NO; >0: means YES
#
# mm_pbsa allows to calculate (absolute) free energies for one
molecular
# species or a free energy difference according to:
#
# Receptor + Ligand = Complex,
# DeltaG = G(Complex) - G(Receptor) - G(Ligand).
#
# PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added
during
# generation of snapshots as well as during mm_pbsa
calculations.
# PATH - Specifies the location where to store or get snapshots.
#
# COMPLEX - Set to 1 if free energy difference is calculated.
# RECEPTOR - Set to 1 if either (absolute) free energy or free
energy
# difference are calculated.
# LIGAND - Set to 1 if free energy difference is calculated.
#
# COMPT - parmtop file for the complex (not necessary for option
GC).
# RECPT - parmtop file for the receptor (not necessary for option
GC).
# LIGPT - parmtop file for the ligand (not necessary for option GC).
#
# GC - Snapshots are generated from trajectories (see below).
# AS - Residues are mutated during generation of snapshots from
trajectories.
# DC - Decompose the free energies into individual contributions
# (only works with MM and GB).
#
# MM - Calculation of gas phase energies using sander.
# GB - Calculation of desolvation free energies using the GB models
in sander
# (see below).
# PB - Calculation of desolvation free energies using delphi (see
below).
# MS - Calculation of nonpolar contributions to desolvation using
molsurf
# (see below).
# If MS == 0, nonpolar contributions are calculated with the
LCPO method
# in sander.
# NM - Calculation of entropies with nmode.
#
PREFIX z_phg_501
PATH ./
#
COMPLEX 0
RECEPTOR 1
LIGAND 0
#
#COMPT ./parm_com.top
RECPT ./parm_z_phg_nosi.top
#LIGPT ./parm_lig.top
#
GC 1
AS 0
DC 0
#
MM 1
GB 1
PB 0
MS 0
#
NM 1
#
################################################################################
@DECOMP
#
# Energy decomposition parameters (this section is only relevant if DC
= 1 above)
#
# Energy decomposition is performed for gasphase energies,
desolvation free
# energies calculated with GB, and nonpolar contributions to
desolvation
# using the LCPO method.
# For amino acids, decomposition is also performed with respect to
backbone
# and sidechain atoms.
#
# DCTYPE - Values of 1 or 2 yield a decomposition on a per-residue
basis,
# values of 3 or 4 yield a decomposition on a pairwise
per-residue
# basis. For the latter, so far the number of pairs must
not
# exceed the number of residues in the molecule
considered.
# Values 1 or 3 add 1-4 interactions to bond contributions.
# Values 2 or 4 add 1-4 interactions to either electrostatic
or vdW
# contributions.
#
# COMREC - Residues belonging to the receptor molecule IN THE
COMPLEX.
# COMLIG - Residues belonging to the ligand molecule IN THE COMPLEX.
# RECRES - Residues in the receptor molecule.
# LIGRES - Residues in the ligand molecule.
# {COM,REC,LIG}PRI - Residues considered for output.
# {REC,LIG}MAP - Residues in the complex which are equivalent to the
residues
# in the receptor molecule or the ligand molecule.
#
DCTYPE 2
#
COMREC 1-166 254-255
COMLIG 167-253
COMPRI 1-255
RECRES 1-168
RECPRI 1-168
RECMAP 1-166 254-255
LIGRES 1-87
LIGPRI 1-87
LIGMAP 167-253
################################################################################
@DELPHI
#
# Delphi parameters (this section is only relevant if PB = 1 above)
#
# The first group of the following parameters are passed to delphi.
# Additional parameters of delphi (e.g. SALT 0.10) may be added
here.
# For further details see the delphi documentation.
#
# FOCUS - If FOCUS > 0, subsequent (multiple) PERFIL and SCALE
parameters are
# used for multiple delphi calculations using the focussing
technique.
# The # of _focussing_ delphi calculations thus equals the value of
FOCUS.
# INDI - Dielectric constant for the molecule.
# EXDI - Dielectric constant for the surrounding solvent.
# PERFIL - Percentage of the lattice that the largest linear
dimension of the
# molecule will fill.
# SCALE - Lattice spacing in no. of grids per Angstrom.
# LINIT - No. of iterations with linear PB equation.
# BNDCON - Type of boundary condition.
# CHARGE - Name of the charge file.
# SIZE - Name of the size (radii) file.
#
# SURFTEN / SURFOFF - Values used to compute the nonpolar
contribution Gnp to
# the desolvation according to Gnp = SURFTEN * SASA +
SURFOFF.
#
FOCUS 0
INDI 1.0
EXDI 80.0
PERFIL 80.0
SCALE 2
LINIT 1000
BNDCON 4
CHARGE ./my_amber94_delphi.crg
SIZE ./my_parse_delphi.siz
#
SURFTEN 0.00542
SURFOFF 0.092
#
################################################################################
@GB
#
# GB parameters (this section is only relevant if GB = 1 above)
#
# The first group of the following parameters are passed to sander.
# For further details see the sander documentation.
#
# IGB - Switches between Tsui's GB (1), Onufriev's GB (2),
# Jayaram's et al. GB (3) or Jayaram's et al. MGB (4).
# SALTCON - Concentration (in M) of 1-1 mobile counterions in
solution.
# EXTDIEL - Dielectricity constant for the surrounding solvent.
#
# SURFTEN / SURFOFF - Values used to compute the nonpolar
contribution Gnp to
# the desolvation according to Gnp = SURFTEN * SASA +
SURFOFF.
#
IGB 4
SALTCON 0.00
EXTDIEL 80.0
#
SURFTEN 0.0072
SURFOFF 0.00
#
################################################################################
@MS
#
# Molsurf parameters (this section is only relevant if MS = 1 above)
#
# PROBE - Radius of the probe sphere used to calculate the SAS.
# RADII - Name of the radii file.
#
PROBE 1.4
RADII ./atmtypenumbers
#
#################################################################################
@NM
#
# Parameters for sander/nmode calculation (this section is only
relevant if NM = 1 above)
#
# The following parameters are passed to sander (for minimization)
and nmode
# (for entropy calculation using gasphase statistical mechanics).
# For further details see documentation.#
# DIELC - (Distance-dependent) dielectric constant
# MAXCYC - Maximum number of cycles of minimization.
# DRMS - Convergence criterion for the energy gradient.
#
DIELC 4
MAXCYC 5000
DRMS 0.1
#
#################################################################################
@MAKECRD
#
# The following parameters are passed to make_crd_hg, which extracts
snapshots
# from trajectory files. (This section is only relevant if GC = 1 OR
AS = 1 above.)
#
# BOX - "YES" means that periodic boundary conditions were used
during MD
# simulation and that box information has been printed in the
# trajecotry files; "NO" means opposite.
# NTOTAL - Total number of atoms per snapshot printed in the
trajectory file
# (including water, ions, ...).
# NSTART - Start structure extraction from NSTART snapshot.
# NSTOP - Stop structure extraction at NSTOP snapshot.
# NFREQ - Every NFREQ structure will be extracted from the
trajectory.
#
# NUMBER_LIG_GROUPS - Number of subsequent LSTART/LSTOP combinations
to
# extract atoms belonging to the ligand.
# LSTART - Number of first ligand atom in the trajectory entry.
# LSTOP - Number of last ligand atom in the trajectory entry.
# NUMBER_REC_GROUPS - Number of subsequent RSTART/RSTOP combinations
to
# extract atoms belonging to the receptor.
# RSTART - Number of first receptor atom in the trajectory entry.
# RSTOP - Number of last receptor atom in the trajectory entry.

# Note: If only one molecular species is extracted, use only the
receptor
# parameters (NUMBER_REC_GROUPS, RSTART, RSTOP).
#
BOX YES
NTOTAL 648
NSTART 501
NSTOP 1901
NFREQ 14
#
NUMBER_LIG_GROUPS 0
LSTART 0
LSTOP 0
NUMBER_REC_GROUPS 1
RSTART 1
RSTOP 648
#
#################################################################################
@ALASCAN
#
# The following parameters are additionally passed to make_crd_hg in
conjunction
# with the ones from the @MAKECRD section if "alanine scanning" is
requested.
# (This section is only relevant if AS = 1 above.)
#
# The description of the parameters is taken from Irina Massova.
#
# NUMBER_MUTANT_GROUPS - Total number of mutated residues. For each
mutated
# residue, the following four parameters must
be given
# subsequently.
# MUTANT_ATOM1 - If residue is mutated to Ala then this is a pointer
on CG
# atom of the mutated residue for all residues except
Thr,
# Ile and Val.
# A pointer to CG2 if Thr, Ile or Val residue is
mutated to Ala
# If residue is mutated to Gly then this is a pointer
on CB.
# MUTANT_ATOM2 - If residue is mutated to Ala then this should be
zero for
# all mutated residues except Thr and VAL.
# A pointer on OG1 if Thr residue is mutated to Ala.
# A pointer on CG1 if VAL or ILE residue is mutated to
Ala.
# If residue is mutated to Gly then this should be
always zero.
# MUTANT_KEEP - A pointer on C atom (carbonyl atom) for the mutated
residue.
# MUTANT_REFERENCE - If residue is mutated to Ala then this is a
pointer on
# CB atom for the mutated residue.
# If residue is mutated to Gly then this is a
pointer on
# CA atom for the mutated residue.
# Note: The method will not work for a smaller residue mutation to a
bigger
# for example Gly -> Ala mutation.
# Note: Maximum number of the simultaneously mutated residues is 40.
#
#NUMBER_MUTANT_GROUPS 3
#MUTANT_ATOM1 1480
#MUTANT_ATOM2 0
#MUTANT_KEEP 1486
#MUTANT_REFERENCE 1477
#MUTANT_ATOM2 1498
#MUTANT_ATOM1 1494
#MUTANT_KEEP 1500
#MUTANT_REFERENCE 1492
#MUTANT_ATOM1 1552
#MUTANT_ATOM2 0
#MUTANT_KEEP 1562
#MUTANT_REFERENCE 1549
#
#################################################################################
@TRAJECTORY
#
# Trajectory names
#
# The following trajectories are used to extract snapshots with
"make_crd_hg":
# Each trajectory name must be preceeded by the TRAJECTORY card.
# Subsequent trajectories are considered together; trajectories may
be
# in ascii as well as in .gz format.
# To be able to identify the title line, it must be identical in all
files.
#
TRAJECTORY ./z_phg_nosi_prod.traj
#TRAJECTORY ../prod_II/md_nvt_prod_pme_02.mdcrd.gz
#TRAJECTORY ../prod_II/md_nvt_prod_pme_03.mdcrd.gz
#TRAJECTORY ../prod_II/md_nvt_prod_pme_04.mdcrd.gz
#TRAJECTORY ../prod_II/md_nvt_prod_pme_05.mdcrd.gz
#
################################################################################
@PROGRAMS
#
# Program executables
#
DELPHI /home/gohlke/src/delphi.98/exe.R10000/delphi
#
################################################################################

:
# MEAN STD
# =======================
ELE -75.92 62.16
VDW -194.13 7.80
INT 977.92 20.92
GAS 707.88 62.69
GBSUR 26.32 0.22
GB -4836.75 58.36
GBSOL -4810.43 58.43
GBELE -4912.66 15.86
GBTOT -4102.56 20.09
TSTRA 15.59 0.00
TSROT 15.21 0.01
TSVIB 496.44 3.49
TSTOT 527.24 3.50
502:
# MEAN STD
# =======================
ELE -75.92 62.16
VDW -194.13 7.80
INT 977.92 20.92
GAS 707.88 62.69
GBSUR 26.32 0.22
GB -4836.75 58.36
GBSOL -4810.43 58.43
GBELE -4912.66 15.86
GBTOT -4102.56 20.09
TSTRA 15.59 0.00
TSROT 15.21 0.01
TSVIB 496.44 3.49
TSTOT 527.24 3.50
503:
# MEAN STD
# =======================
ELE -75.92 62.16
VDW -194.13 7.80
INT 977.92 20.92
GAS 707.88 62.69
GBSUR 26.32 0.22
GB -4836.75 58.36
GBSOL -4810.43 58.43
GBELE -4912.66 15.86
GBTOT -4102.56 20.09
TSTRA 15.59 0.00
TSROT 15.21 0.01
TSVIB 496.44 3.49
TSTOT 527.24 3.50

Also this is the output file when I started at frame 1 and went to
frame 10 skipping every other one:
# MEAN STD
# =======================
ELE 41.07 32.31
VDW -187.69 8.61
INT 973.15 23.13
GAS 826.53 32.39
GBSUR 26.27 0.15
GB -4953.70 22.27
GBSOL -4927.43 22.16
GBELE -4912.64 15.05
GBTOT -4100.90 22.21
TSTRA 15.59 0.00
TSROT 15.19 0.00
TSVIB 499.44 2.76
TSTOT 530.23 2.77