Molecular Dynamics

A List of Solutions to Some WARNINGs in CHARMM Simulation Package

2010-05-17T18:40:00Z

There are many types of warnings that appear in CHARMM. Here are the possible solutions for some of them.

1. "Element negative. No operation" from the SCALAR module

e.g.

CHARMM>    scalar sca4 sum sca5 ! + <y^2>

CHARMM>    scalar sca4 sum sca2 ! - <y>^2

CHARMM>    scalar sca4 sum sca6 ! + <z^2>

CHARMM>    scalar sca4 sum sca3 ! - <z>^2

CHARMM>    scalar sca4 sqrt

      ***** LEVEL 0 WARNING FROM <SCALAR> *****
      ***** Element negative. No operation
      ******************************************
      BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

[[ Solution ]]: This is most likely due to the fact that some atoms' coordinates are missing somehow. Try to find out these atoms and fill their internal coordinates.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2. "zero selection specified" when using HBUILD

e.g.

CHARMM>    !Rebuild hydrogen coordinates
CHARMM>    HBUIld SELE ALL END
SELRPN>    667 atoms have been selected out of    667

      ***** LEVEL 4 WARNING FROM <GTHBUI> *****
      ***** Zero selection specified
      ******************************************
      BOMLEV ( 0) IS NOT REACHED. WRNLEV IS 5

[[ Solution ]]: change the selection to hydrogen atoms only, i.e.

CHARMM>    !Rebuild hydrogen coordinates
CHARMM>    HBUIld SELE hydrogen END
SELRPN>    344 atoms have been selected out of    667
PRNHBD: CUToff Hydrogen Bond distance =    0.5000   Angle =   90.0000
         CuT switching ON HB dist. =     3.5000 OFf HB dist. =    4.0000
         CuT switching ON Hb Angle =    50.0000 OFf Hb Angle =   70.0000
         ACCEptor antecedents included
         All hydrogen bonds for each hydrogen will be found
         Hydrogen bonds between excluded atoms will be kept

[[ Reason ]]: (from CHARMM hbuild.doc)

"<atom-selection> specify the hydrogens to be (re-)constructed."

"By default (if no selection is specified) these are all unknown hydrogens and lone pairs (this is equivalent to a selection "SELEction (LONE .OR. HYDRogen) .AND..NOT INITial")."

[[ Remarks ]]: What is said in hbuild.doc may be true. However, if one does not make any selections after the HBUILD keyword, the same warning will still appear.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3. "Failed to allocate memory for natomingp, gptoalst" from gbsw.src

e.g.

***** LEVEL -1 WARNING FROM <GB_lookup>gbsw.src *****
      ***** Failed to allocate memory for natomingp,gptoalst
      ******************************************
      BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

[[ Solution ]]; Most likely this is due to the missing coordinates for some atoms. Try to find these atoms and fill their coordinates.

4. "Zero length string being converted to 0" from PHMD

e.g.

CHARMM>    label PERLDONE

CHARMM>    phmd par 23 wri 25 ph 7.0 npri 100 mass 10 barr 1.75 bartau 2 temp 300 lam sele .not. ( resn TIP3 .or. RESN CLA .or. RESN SOD ) end

PHMD> Continous constant pH molecular dynamics

SELRPN>     36 atoms have been selected out of   2163
TITLE> * STATES FILE FOR HYBRID-PHMD
TITLE> * SYNTAX:
TITLE> * RESNAME MODEL_PKA PARA PARB
TITLE> * ATOM_TYPE CHARGE(1) CHARGE(2) [RAD(1) RAD(2)]
TITLE> *
WARNING from DECODF -- Zero length string being converted to 0.
PHMD> Grp # 1:ARG (ARG    1) pKa= 12.50 Prot
PHMD> Grp # 1:ARG (ARG    1) HH12

[[ Solution ]]: double check the phmd titration parameter file. Most likely its format is not correct, e.g. leaving out the barrier potential.

results:

CHARMM>    phmd par 23 wri 25 ph 7.0 npri 100 mass 10 barr 1.75 bartau 2 temp 300 lam sele .not. ( resn TIP3 .or. RESN CLA .or. RESN SOD ) end

PHMD> Continous constant pH molecular dynamics

SELRPN>     36 atoms have been selected out of   2163
TITLE> * STATES FILE FOR HYBRID-PHMD
TITLE> * SYNTAX:
TITLE> * RESNAME MODEL_PKA PARA PARB BARR
TITLE> * ATOM_TYPE CHARGE(1) CHARGE(2) [RAD(1) RAD(2)]
TITLE> *
PHMD> Grp # 1:ARG (ARG    1) pKa= 12.50 Prot
PHMD> Grp # 1:ARG (ARG    1) HH12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

e.g.

CHARMM>    label PERLDONE

CHARMM>    nbond elec atom cdie shift vdw vswitch ctonnb 10.0 ctofnb 12.0 cutnb 14.0 cutim 14.0 inbfrq -1 imgfrq -
1 ewald pmew fftx 64 ffty 64 fftz 64 kappa .34 spline order 6
<PME> Total heap storage needed =     1034887
Fill Ewald table: Number of points=     10000 EWXmax=    4.250000
fill erfc table: linear inter has rms error = 0.979220D-08 maximum error = 0.218740D-07
fill erfc table: cubic spline has rms error = 0.360331D-11 maximum error = 0.108438D-10

NONBOND OPTION FLAGS:
     ELEC     VDW      ATOMs    CDIElec SHIFt    VATOm    VSWItch
     BYGRoup NOEXtnd EWALd
CUTNB = 14.000 CTEXNB =999.000 CTONNB = 10.000 CTOFNB = 12.000
WMIN   = 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS    = 1.000
NBXMOD =      5
PME EWALD OPTIONS: KAPPA = 0.340 QCOR = 0.000 Bspline order = 6
FFTX= 64 FFTY= 64 FFTZ= 64
                Using Pub FFT
                Real-to-Complex FFT
There are        0 atom pairs and        0 atom exclusions.
There are        0 group pairs and        0 group exclusions.
<MAKINB> with mode   5 found 27483 exclusions and   2472 interactions(1-4)
<MAKGRP> found    733 group exclusions.
Generating nonbond list with Exclusion mode = 5
== PRIMARY == SPACE FOR 7400454 ATOM PAIRS AND        0 GROUP PAIRS
VEHEAP> Expanding heap size by        11436032 words.
== PRIMARY == SPACE FOR 11100701 ATOM PAIRS AND        0 GROUP PAIRS
VEHEAP> Expanding heap size by        16990208 words.
== PRIMARY == SPACE FOR 16651072 ATOM PAIRS AND        0 GROUP PAIRS

General atom nonbond list generation found:
11314422 ATOM PAIRS WERE FOUND FOR ATOM LIST
   641446 GROUP PAIRS REQUIRED ATOM SEARCHES

**** Warning **** The following extraneous characters
were found while command processing in CHARMM
CUTIM 14.0 IMGFRQ -1

There are many types of warnings that appear in CHARMM. Here are the possible solutions for some of them.

1. "Element negative. No operation" from the SCALAR module

e.g.

CHARMM>    scalar sca4 sum sca5 ! + <y^2>

CHARMM>    scalar sca4 sum sca2 ! - <y>^2

CHARMM>    scalar sca4 sum sca6 ! + <z^2>

CHARMM>    scalar sca4 sum sca3 ! - <z>^2

CHARMM>    scalar sca4 sqrt

      ***** LEVEL 0 WARNING FROM <SCALAR> *****
      ***** Element negative. No operation
      ******************************************
      BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5

[[ Solution ]]: This is most likely due to the fact that some atoms' coordinates are missing somehow. Try to find out these atoms and fill their internal coordinates.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2. "zero selection specified" when using HBUILD

e.g.

[[ Solution ]]: change the selection to hydrogen atoms only, i.e.

[[ Reason ]]: (from CHARMM hbuild.doc)

"<atom-selection> specify the hydrogens to be (re-)constructed."

"By default (if no selection is specified) these are all unknown hydrogens and lone pairs (this is equivalent to a selection "SELEction (LONE .OR. HYDRogen) .AND..NOT INITial")."

[[ Remarks ]]: What is said in hbuild.doc may be true. However, if one does not make any selections after the HBUILD keyword, the same warning will still appear.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3. "Failed to allocate memory for natomingp, gptoalst" from gbsw.src

e.g.

[[ Solution ]]; Most likely this is due to the missing coordinates for some atoms. Try to find these atoms and fill their coordinates.

4. "Zero length string being converted to 0" from PHMD

e.g.

[[ Solution ]]: double check the phmd titration parameter file. Most likely its format is not correct, e.g. leaving out the barrier potential.

results:

5. "extraneous characters" from NBOND settings

e.g.

[[ Comments ]]: CUTIM should not appear in NBOND settings . The actual value of CUTIm should be calculated by "crystal build" based on the information on the water box dimensions (see below).

[[ Solution ]]: Thus CUTIM should be removed from the NBOND settings and IMGFRQ can be kept. After this modification, no more warning is issued.

CHARMM>    label PERLDONE

CHARMM>    crystal defined octa 70.488 70.488 70.488 109.4712206344907 109.4712206344907 109.4712206344907
Crystal Parameters : Crystal Type = OCTA
           A     =   70.48800 B    =   70.48800 C     =   70.48800
           Alpha = 109.47122 Beta = 109.47122 Gamma = 109.47122

CHARMM>    label PERLDONE

CHARMM>    crystal build cutoff 35.244 noper 0
XBUILD> Building all transformations with a minimum atom-atom
         contact distance of less than   35.24 Angstroms.

Range of Grid Search for Transformation     1 :
Lattice Vector A    -3 TO     3
Lattice Vector B    -3 TO     3
Lattice Vector C    -2 TO     2

Note: the actual value of CUTIm is 35.244 in this case.

   NSTEP =      500    NSAVC =      500    NSAVV =        0
ISCALE =        0   ISCVEL =        0   IASORS =        1
IASVEL =        1   ICHECW =        1   NTRFRQ =      500
IHTFRQ =        0   IEQFRQ =       50   NPRINT =      500
INBFRQ =       -1   IHBFRQ =        0   IPRFRQ =      500
ILBFRQ =       50   IMGFRQ =       -1    ISEED =          2023018546
Note: the value for IMGFRQ is "-1".

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Secondary Structure Bias in Generalized Born Solvent Models

2010-03-13T14:04:00Z

Normal 0 false false false EN-US ZH-CN X-NONE

Secondary Structure Bias in Generalized Born Solvent Models: Comparison of Conformational Ensembles and Free Energy of Solvent Polarization from Explicit and Implicit Solvation

Daniel R. Roe, Asim Okur, Lauren Wichstrom, Viktor Hornak, and Carlos Simmerling. J. Phys. Chem. B, 2007, 111, 1846-1857

To simulate the protein behavior in an aqueous environment, both implicit and explicit solvent models can be used. While explicit solvent models are more realistic and physically rigorous, implicit solvent models have their own attractiveness: 1) Without including solvent molecules, the implicit simulations can considerably reduce the system size and thus significantly lower computational cost. 2) Conformational sampling is increased due to accelerated molecular motions.

In the implicit solvent model, the free energy of solvation contains nonpolar and polar components. While the nonpolar part only need to deal with the surface area, the polar part has to consider the solute charges. The most accurate method of calculating polar solvation energy is by solving the Poisson equation (PE) but at a higher cost. The generalized Born (GB) implicit solvent models are based on the PE model but with several approximations to increase the speed of calculation.

The question is the accuracy of the GB models. Studies have shown that GB calculations can give significant errors due to Coulomb-field approximation, overstabilize ion pair interactions and α-helical conformations. In this paper, detailed comparisons of conformational ensembles and free energy of solvent polarization were provided. The model protein is Ala10 with two different starting conformations (extended and collapsed). Replica exchange molecular dynamics (REMD) were performed using TIP3P explicit solvent model, PE and three GB implicit solvent models (GBHCT, GBOBC, and GBNeck) using Amber 9 force field.

Results showed that explicit solvent model had Ala10 predominantly adopt PPII conformation, which is in agreement with experimental observations. But the GB models predicted significant different secondary structure populations. In particular, GBHCT and GBOBC models contained overabundance of helical structure. These results have significant implications for the use of GB models for structure prediction or characterization of folding landscapes.

Using the TIP3P model as the standard, the free energies of solvent polarization were compared according to four different conformations of Ala10: α-helix (alpha), left-hand α-helix (left), β-hairpin (hairpin), and PPII. Overall, the PE model has the best performance while the GB models are found to be conformational dependent. For the GB models, agreement with TIP3P was best for the well-solvated PP2 conformation, while growing worse for more compact conformations (hairpin, left, and alpha). The calculated effective Born radii for GB and PE shows small deviations for PPII while significant deviations were found for the more compact conformations. In general, an atom whose effective radius has been underestimated will have an overestimated solvation free energy and vice versa.

Both the secondary structure and the free energy of solvation are very important for the folding studies. The secondary structure bias was also observed in our simulations on (AAQAA)₃ with implicit solvent model which indicated high α-helical percentage. So does the HP36 fragments. Alternatively, by optimization of the GB input radii through examining pairwise interaction between amino acid polar groups, the conformational equilibria of (AAQAA)₃ peptide and the the GB1p series β-hairpins can be successfully reproduced (ref).

Ref: Jianhan Chen, Wonpil Im, and Charles L. Brooks III. Balancing Solvation and Intramolecular Interactions: Toward a Consistent Generalized Born Force Field, JACS, 2006, 128, 3728-3736

Monte-carlo sampling of alchemical degrees of freedom

2010-02-21T17:21:00Z

Molecular dynamics is routinely used to calculate free-energy changes associated with alchemical transformations. The free-energy differences can be used for many practical purposes such as drug screening, pKa calculations, and determining solvation energies, to mention but a few. Although the typical free-energy techniques (free-energy pertubartion and thermodynamic integration) are straight forward and deliver accurate free-energy changes, in many cases it may not be necessary to obtain accurate free-energy differences, but instead only the order of binding strength may be desired. For instance, one may want to determine which molecule out of several binds most strongly to a specific site in a protein for drug screening. For this type of comparison, it is necessary to know the free-energy changes of both removing the molecules from solvent and then binding to the host site. For two target molecules this would require four independent free-energy calculations.

To reduce the computational burden imposed by such calculations, it is desirable to couple binding and desolvation simulations of both target molecules. By coupling the manifolds of all end states it may then be possible to determine the most efficient binder directly from our simulation more efficiently. One method to couple the required states introduced by Kong and Brooks (J. Chem. Phys. 1996), is to first determine desolvation free-energies for each compound of interest and build these free-energies into a Hamiltonian that is then used to evaluate the preferential binding of the ligand. In the method of Kong and Brooks (Lambda-dynamics) the alchemical degree of freedom was treated as a dynamic variable, and forces applied to the lambda degree of freedom were used to determine the order of binding strength.

The general idea of building in desolvation energies into a competitive binding simulation has also been applied by Pitera and Kollman, but cast in a slightly different form. In the the chemical MC/MD method of Pitera and Kollman (J. Am. Chem. Soc. 1998), desolvation energies are first calculated. MD simulations are then run for a short period of time and periodically interrupted, and the question is asked whether a new ligand should replace the old ligand. The criteria for making the move is determined by the typical monte-carlo criteria, but the energy difference is augmented by the known desolvation energy differences. By this approach the simulation is analogs to an experimental binding competition experiment. Their approach can naturally be used for multiple possible ligands. In testing their methodthe binding strength for several ligands was tested for Rebek's "Tennis ball". The method quickly determined the correct order of binding, and with longer simulations delivered quantitative free-energy differences. In contrast to the continuous lambda-dynamics approach the MC/MD approach saves time by avoiding intermediate states, but may require longer simulations to get converged statistics due to possible low ligand exchange acceptance. Another method that uses a hybrid MC/MD protocol and coupled multiple free energy simulations is the work of Jarque and Tidor (J. Phys. Chem. B 1997).

In constrast to the MC/MD work of Pitera and Kollman this procedure required the simultaneous simulation of both environments, in this case the gas phase and solution, to evaluate the relative solvation energies of small molecules (methane/methanol). These simulations were periodically interrupted and a MC move was attempted. Instead of exchanging end states only however, intermediate states were attempted to improve the acceptance ratio. The acceptance criteria now was based on the combined energy difference from simulations in both environments. Another trick used by Jarque and Tidor to determine the most stable chemical state from both simulations was to use simulated annealing in the MC-move. This allowed very rapid determination of the most stable species (the species with the greatest solubility in aqueous solution).

The method of coupling MD to monte-carlo sampling of the alchemical degrees of freedom has been shown to be successful for small uncomplicated systems, and can provide a fast ordering of species stability. These methods generate an ensemble with MD and then use MC so sample the different energy surfaces generated by different ligands in the ensembles. It remains to be shown how this method would perform for more complicated systems, and systems with large energy differences between the end states. It is likely that for more complicated systems more traditional free-energy techniques may be superior due to low acceptance of exchange.