Computational Chemistry & Biophysics Connections

A List of Solutions to Some WARNINGs in CHARMM Simulation Package

StevenWang — Mon, 17 May 2010 18:40:00 GMT

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

cshi — Sat, 13 Mar 2010 14:04:00 GMT

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

Jason — Sun, 21 Feb 2010 17:21:00 GMT

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.

Book: Genetics and Molecular Biology 2nd Ed. by Robert Schleif (JHU)

StevenWang — Wed, 13 Jan 2010 19:10:00 GMT

With the generosity of Prof. Robert Schleif (Biology department, John Hopkins University), there is another book of him that can be accessed for free. The book "Genetics and Molecular Biology 2nd Ed" can be download from his website: http://gene.bio.jhu.edu/pub.html

Book: Analysis of Protein Structure and Function: A Beginner's Guide to CHARMM by Robert Schleif (JHU)

StevenWang — Wed, 13 Jan 2010 19:04:00 GMT

Prof. Robert Schleif (Biology department, John Hopkins University) has made his book "Analysis of Protein Structure and Function: A Beginner's Guide to CHARMM" available for free download: http://gene.bio.jhu.edu/pub.html

He also provided CHARMM input scripts free for download: http://gene.bio.jhu.edu/scriptindex.html

Steven

Molecular Driving Forces: statistical thermodynamics for chemistry and biology

Deepthi — Wed, 23 Dec 2009 20:25:00 GMT

I used this book for my course on statistical mechanics. It is a very good introduction to statistical mechanics. The use of lattice model in explaining various topics made it easy for us to understand the underlying concepts. It starts with a simple introduction to probability, equilibrium, math tools such as series and approximations, and then goes on to Boltzmann distribution law, entropy, free energies, statistical mechanics, phase transitions, the details of how water is different from other substances, polymer solutions etc. I liked the some of the chapters discussed in the book like the entropy, free energy of mixing, phase transitions and the topics on water because of the way they were interpreted in detail. I was fascinated by the anamolous properties of water due to its hydrogen bonding network explained using statistical thermodynamics. The chapters were well divided into microscopic and macroscopic properties of various substances.This book is a very good book for beginners tol earn about statistical mechanics. I had previous knowledge about many of the concepts discussed but still I liked the way the author explained the different concepts using various models. But I thought that the exercise problems towards the end of each chapter were a bit confusing due to the lack of sufficient information in some of the problems. Other than that the book was really good.

Comments on "Molecular driving forces: statistical thermodynamics in chemistry and biology" by Ken A. Dill and Sarina Bromberg.

Linh — Fri, 18 Dec 2009 17:21:00 GMT

"Molecular Driving Forces: Statistical Thermodynamics in Chemistry and Biology" is a moderately good book. It really depends on your purposes and expectations whether you should choose this book or not.

This book covers an impressive amount of material on statistical thermodynamics in both chemistry and biochemistry areas. The material that Ken and Sarina try to present in this book is a little overwhelming for a one semester course. They have gone over most of the classical concepts that need to develop the statistical mechanics understanding vary from easy to difficult level. On the top of that, they apply the statistical thermodynamics onto biological systems such as proteins and nucleic acids, which is tremendously helpful for people who are interested in studying the physical properties of biological systems. However, the book is too spread out since it tries to combine two large areas into one. I would think if they split one good book into two better books, one is focusing on the statistical thermodynamics concepts and theories itself and one is focusing on the applications particularly in biological systems, they can discuss both subjects at a more thoroughly detailed level. Another thing that bother me just as much is the book completely ignores the quantum mechanics and quantum statistics, not even mention them in introductions or reviews sections. The idea of quantum statistics might be frightening to me as well as many students, but in order for students to have a better appreciation, a better understanding of a bigger picture, I suppose they should dedicate a chapter to discuss about quantum mechanics and quantum statistic. The last thing that I dislike about this book is the problems set at the end of each chapter. The problems are ranging from easy "plug in numbers and get answers" problems to extreme difficult derivations problems. I didn't find many intermediate levels problems. A lots of the problems are just testing whether one can interpret a plot or figures. Some of the questions are almost helpless in trying to inspire you develop your critical thinking skills.

On the bright side, I have found the book of Ken and Sarina to be a very good book for self-study. It has some very well done reviews on mathematical tools that help you develop mathematical skill in studying statistical mechanics. Chapter 1 is called "principles of probability", which goes over the fundamental of probability concepts and thus lead to the foundation of entropy. Chapter 4 is another math review on series and approximations. Follow by chapter 5, which is a review on multivariate calculus. Chapter 17 is a vector calculus reviews and Gauss's theorem. These chapters makes the books qualify to be a very suitable self-study book. I am very grateful for these so called reviews; a lots of them were new to the students. The chapter about Boltzmann Distribution Law (chapter 10) was done very nicely. The detailed, well put and very organized explanations and derivations about Boltzmann Distribution Law and Partition Function was impressively presented. This chapter is one of many of the chapters in this book where statistical thermodynamics concepts are explained and derived slowly and clearly and easily to understand. Overall, I like the book. It is a moderately good book like I said in the beginning. It has the potential to be a great book. I was just slightly upset because there was so much material was mentioned in the book and that make it hard for learners to keep up; it was a little crowded and overwhelming.

Book Review: K. Dill and S. Bromberg, "Molecular driving forces: statistical thermodynamics in chemistry and biology", Garland Science Publisher, 2003

StevenWang — Fri, 18 Dec 2009 01:13:00 GMT

Molecular driving forces: statistical thermodynamics in chemistry and biology

by Ken A. Dill and Sarina Bromberg. Garland Science (Taylor & Francis Group), New York, 2003. 686 pp. $109

The book "Molecular driving forces: statistical thermodynamics in chemistry and biology" by Ken Dill and Sarina Bromberg, published by the Garland Science publisher in 2003, is the most comprehensive and useful book I have ever read on the the topic of thermodynamics.This book starts from an elementary introduction of probability calculation, and gradually extends to the free energy and entropy concepts, based on which a lot of phenomena are clearly explained (both microscopically and macroscopically), like surface tension, diffusion, catalysis, phase transition, electrostatics, cooperativity, etc. A special topic on polymer science is also discussed (Dr. Ken. A. Dill used to be a postdoc of Dr. Paul J. Flory (1, 2)), which is not common for a textbook on this topic.There are several merits that make it a very good textbook. First, the prerequisite requirement is low. All mathematical tools are included in the book chapters or appendixes. Secondly, the authors have a very clear thought on explaining the topic. Each chapter is based on the text in the previous chapters and all of them are well connected. The reader will feel no gap between them. Thirdly, the authors use a very friendly style in discussing the topics. Also, the authors are frank about the knowledge in discussion. No obscure expression (linguistically or mathematically) is used. Everything is expressed in a direct and easy to understand manner. But there are also some flaws of this book. Besides some typos (not intended by the authors), some exercise questions are not well stated. Sometimes the reader has to guess what the authors really want to say. Also, the book does not cover the non-equilibrium thermodynamics, therefore some recent advances in statistical thermodynamics are not discussed. Although with these little flaws, this is really a good book to study if one is interested in learning thermodynamics.

References:

1.The homepage of Ken. A. Dill's group (http://www.dillgroup.ucsf.edu/)

2. Dill and Flory. Molecular organization in micelles and vesicles. Proc. Natl. Acad. Sci. U S A. 1981, 78:676–680.

Steven Wang

e-mail: stevenaura [at] ou.edu

December 17, 2009

Free energy reconstruction from nonequilibrium single-molecule pulling experiments

kxia — Mon, 20 Apr 2009 20:07:00 GMT

Recently, the development of single molecule manipulation has led new sight into mechanical properties of individual molecules. However, the measurement in experiment drive systém away from equilibrium. How can one extract equilibrium properties? According to the second thermodynamics law, the mechanical work is always larger than the free energy, except that the experiment is performed reversibly or infinitely slowly. However, recently Jarzynski discovered the identity between thermodynamic free energy differences and the irreversible work. Thus one should be able to extract free energy surfaces from the atomic force measurements in repeated pulling single molecule experiments.

In the single molecular force measuring experiment, the sample is moved with a constant speed relative to the cantilever with spring constant k. The position of the cantilever with respect to the sample is recorded. After repeating this measurement a certain times, the free energy profile of the equilibrium systém can be obtained exactly according to Jarzynski's equality. This reconstruction result fits well with the reconstruction curve. Through measuring the pulling force, one can also estimate the free energy profile through integrating. However, this introduces some artifacts and should not be considered accurate.

In another experiment, the free energy profile from Jarzynski's equality is compared with that from fluctuation-dissipation relation. The results show that in reversible process both of them fits well with theoretical predictions for near-equilibrium conditions. In irreversible process, the free energy profile from Jarzynski's equality performs consistently over the entire range, while the free energy profile fails at large extension range. The effect of repeating times is also tested. The results shows that the more times the process repeats, the better the free energy profile recovers.

In conclusion, Jarzynski raised an equality relating the irreversible work to the equilibrium free energy difference. And this is approved by the recently developed single molecule manipulation experiments.

Reference: Gerhard Hummer and Attila Szabo, PNAS, 98, 7, 2001

Jan Liphardt etc, Science, 296, 1832, 2002

Modeling Materials Properties Without Experimental Input

kxia — Mon, 20 Apr 2009 20:01:00 GMT

Empirical physical models rely on parameters from experiments, and so induce various inaccuracies. However, a quantum mechanical model can offer independent data. Also, QM is able to predict everything theoretically. But QM has properties not easily addressed, which depends on complexities at larger length scales over heterogeneities. And, there is no universal QM method appropriate for all materials and phenomena.

Ab initio post-Hartree-Fock quantum chemistry and quantum Monte Carlo simulation are the most accurate. They make no other physical approximations and can get accurate ground and excited state properties. However, the main drawback of such methods is computational expense. Therefore, more common methods for QM modeling uses or builds upon density functional theory(DFT). It is much simpler and less costly. DFT is a formally exact ground state theory in which the material's energy is expressed as a fun of the electron density alone. The primary disadvantage of current DFT methods is the approximate XC functional. However, for most ground state properties, the generalized gradient approximation(GGA) XC functionals provide sufficient accuracy. Hybrid XC DFT-GGA techniques are developed, such as DFT+U, include some exact HF exchange, and are suitable for description of mid-to-late first row transition metal oxides and sulfides, but not appropriate for metals. TD-DFT can be used to calculate electronic or optical spectra of materials and GW method can be used to obtain ionization energies and electron affinities. BSE takes DFT and GW data as input and accounts for electron -hole interactions. With these methods in mind, we can use the appropriate method to predict a given property for a given materials class as a function of accuracy and expense.

Amorphous structures are difficult to model with QM because the usual 3D periodic boundary conditions introduce correlation length artifacts and it is certain that a random amorphous structure generated. Heterogeneous mixtures offer the most severe challenge for future materials modeling. Multiscale modeling aims to bridge length and time scales to make overarching predictions of materials behavior. Major unsolved issues in this area include how to transfer heat and mass across all scales and etc.

Now the typical simulation still starts with guidance from experiment regarding approximate initial structure and composition, but given such guidance, QM can provide sight to how properties change when composition and structure change, thereby furthering atomic-scale manipulation of material design.

Reference: Emily A. Carter, Science, 321, 800, 2008

Recent Advances in the Generalized Born Method

Jason — Mon, 20 Apr 2009 14:54:00 GMT

Due to the high computational costs of explicit solvent simulations, the development of accurate implicit solvent models is an attractive approach to reach longer timescales in biomolecular simulations. Possibly the most crucial aspect an implicit solvent model is the calculation of the electrostatic solvation energy. To completely describe the electrostatic solvation energy the response of the solvent to the solute and vice versa should be considered, but in empirical force-field calculations the solute charge response to the solvent is not considered. The benchmark of continuum electrostatic calculations is the Poisson-Boltzmann (PB) equation, where a finite difference approach is used to solve for the potential. This method, although accurate, is computationally demanding and does not lend itself to energy minimization or dynamics because the forces cannot be calculated analytically. Recent development of the generalized born approach, a pairwise-approximation to the PB method, has been shown to reproduce remarkably well the PB electrostatic solvation energy if one uses the same definition for the molecular volume. In 2002 Lee et. al. developed new grid-based and analytical born models that empirically correct for the non-spherical nature of the molecular volume. These models were shown to have excellent agreement with PB calculations. In 2003 an extension to the previous work of Lee et. al. refined their empirical correction for the molecular volume, and included a vector based scaling approach to correct error in the standard overlapping atomic functions approach. They also developed an accurate solvent-accessible surface area approximation to account for the nonpolar contribution to the solvation energy that is based on the same computational machinery as their GB model. This recent method by Lee et. al. has been shown to be one of the most accurate GB models to date. Lee, Salsbury and Brooks, J. Chem. Phys. 116: 10606-10614 (2002) Lee, Feig, Salsbury and Brooks, J. Comput. Chem 24:1348-1356 (2003)

Normal Mode Analysis (Tama and Sanejouand, Protein Eng. 2001; Dobbins et al, PNAS 2008)

Linh — Fri, 17 Apr 2009 18:33:00 GMT

Proteins with many transitional conformations when changing from one conformation to another usually involve relative motions of domains. A good understanding of these movements will help the study of protein functions such as catalysis and regulations. In these papers, a theoretical method for studying collective motions in proteins, normal mode analysis (NMA) is used. With this tool, conformational changes can be expressed in term of superposition of collective variables, ie normal mode coordinates.

Normal Mode Theory:

The displacements of atomic coordinate i, r_i(t) with near by stationary point of the potential energy surface can be calculated using Goldstein equation (1950). Within this approach, a simplified potential energy function proposed by Tirion (1996) is used. This simplified potential energy function is designed such that any configuration of any system is a minimum of the function. With this method, NMA no longer require energy minimization step and thus, cut off CPU time. Another two techniques of the NMA in this study are all the normal mode calculation used a cutoff of 8 Å which has been used by Bahar et al (1997) and only C_α atoms were taken into account. This model allowed the study of backbone motions, which characterizing low frequency normal modes of large proteins, to be achieved with small amount of CPU time and lower computational expense. The simple potentials and models have been used in previous studies, and the results strongly suggested that simple potentials yielded low frequency normal mode as accurate as those with standard NMA.

In order to compare with the experimental method, the authors looked at two qualities, the overlap I_j and the correlation coefficient C_j where I_j is the overlap between delta_vector_r = {Δr₁, ..., Δr_i, ..., Δr_3N}, the conformation change in crystal structures and the jth normal mode of the protein; C_j is the similarity of the patterns of atomic displacements in the conformational change and in mode j. Equations for quantitatively calculations of C_j and I_j are given in the paper (p2-3). Note that delta_vector_r is calculated for pairs of crystallographic conformations of “open” and “closed” form after both conformations were superposed. The overlap and correlation are considered as functions of degree of collectivity (к) of the conformational change where degree of collectivity is a measurement of how collective a motion of the protein is (Bruschweiler, 1995). Overlap of the most involved modes in conformational changes calculated using simplified model are found to be almost equivalent as those obtained using standard NMA and to have smaller RMSD 1.2 – 1.9 Å comparing to 2-3 Å obtained in standard approach. The results also suggested that NMA performed better in the open form than closed (refer to table III), better with highly collective motions (refer to table IV) and better with more localized motions (refer to the case of tryglyceride lipase in table IV and V). In conclusion, a single normal mode can carry a lot of information about conformational transitions of proteins. The single mode with best correlation and best overlap is more likely not the lowest frequency mode but often one of the three lowest frequency modes.

Carrying on the NMA to study conformational change of protein, the work in the second paper performed NMA calculations to investigate the extend to which thermal motions of proteins could provide the general nature of conformational change upon protein-protein docking.

A set of 20 proteins with flexibility which have been observed to have large conformational change under forming complexation by Tama et al, is used in this study. Their result showed high agreement with the experiment in term of average collectivity, 0.40 for observed and 0.38 for predicted. They also stated that the RMSD of binding sites are found to be 1.2 times larger than those over the entire protein, which suggested that upon complexation, the protein conformational change has highly flexible regions around binding site. Further investigation of mobilities over the binding site for both observed and predicted using NMA were performed. Splitting the binding site into the core and the periphery residues, they found that the observed motion of the peripheral residues is more than twice large than that of the core residues. This confirm the hypothesis that binding site has regions of highly flexibility along with rigid regions upon docking (refer to figure3). There were no obvious trend for the overlap. They were observed to be different for different frequency modes. However, their results are in agreement with previous study that large conformational change of proteins usually have low frequency modes but not the lowest frequency mode. The protein size are also examined in this paper. They proposed that the larger the protein size the more collective motions it has and thus might also be used to predict conformational change. However the size is less reliable than the normal mode frequency.

Overall, this particular study showed that almost all proteins undergo large conformational change has thermal motions in isolation and required some intrinsic flexibility. In their prediction using NMA, only one third of the proteins have the direction and location of motion among the lowest modes is similar to the observed conformational transition. They suggested this because of the bound condition which might altered structural changes. In my opinion, this is one of the limitations of this work, along with unclear trends for overlap. In conclusion, the significant contribution of this study is that it suggested a model to study the conformational change during molecular recognition will be desirable in future work.

References:

Tama and Sanejouand, Protein Eng., 2001

Dobbins, Lesk, and Sternberg, PNAS, 2008

λ-dynamics method and applications (Kong and Brooks, JCP 1996; Khandogin and Brooks, Biophys J 2005)

StevenWang — Thu, 16 Apr 2009 14:52:00 GMT

λ-dynamics is a new method developed on the basis of the widely used free energy perturbation
method (FEP) and the umbrella sampling method. It differs from the FEP method by introducing
multiple parameters λi, instead of only one λ, which are assigned to different types of energy
components. In order to surmount the possible energy barriers and thus get more samples,
special biasing potentials are used based on the idea of umbrella sampling. The forms of these
biasing potentials are quite innovative, which look like kinetic energy and potential energy for
factious particles where λi, instead of r, serve as the coordinates and a set of fake masses of
particles are introduced. The benefit from introducing more λ is that more reaction pathways are
accessible during the simulation. On the other hand, the umbrella potential term is useful in
flattening the energy curve of the reaction pathway, preventing the boundary crossing of λi (e.g.
λi < 0 or λi > 1) and limiting the range of {λi}. Additionally, the fake masses make the control of
λi more flexible. In a word, these new variables and terms make the simulation more adjustable
and easier to control for different purposes. The results of this new method are at least as
accurate as those by the single λ FEP method, which are verified by simulations of relative free
energies of hydration and relative binding affinity simulations.

Constant pH molecular dynamics (CPHMD) is one of the applications of the λ-dynamics
approach. In CPHMD method, n λi are assigned to n titrating residues, which varies from 0 to 1
when the titrating residue change from protonated form to deprotonated form and it takes the
mathematical form as an sine function. A two-dimensional λ-dynamics method is used here for
the purpose of take the proton tautomerism effect into consideration. The two dimensions are
due to the two reaction coordinates: one is for the deprotonation process (λ); the other is for the
tautomeric effect (x). These two coordinates are used throughout the potential energy
calculation via λ-dynamics method. When calculating van der Waals energy, pairwise
interactions are divided into two classes: interactions between titrating residues and other
resides and those between two titrating residues. For coulomb and GB electrostatic energies,
atomic partial charges on titrating residues are changed along the reaction coordinates. Biasing
potentials are also used to facilitate the simulation. One is like harmonic umbrella potential and
the other one originates from the calculation of free energy based on experimental pKa value.
In CPHMD, the unprotonated form has λ > 0.9, while that of the protonated form is smaller than
0.1. In both cases, the x value is either larger than 0 or smaller than 1 by a maximum magnitude
of 0.1, in other words, the tautomer state must be sufficiently pure. The enhance the sampling
near the two ends of the reaction coordinates, a barrier potential is introduced. This potential
function is symmetric considering the middle point of the reaction coordinate. The result is that
the energy is lower when the titrating residues are in pure states than that of the titrating
residues in mixed states and therefore more pure sate residues are sampled.

In my opinion, the CPHMD method makes a best use of the advantage of the λ-dynamics
method, making the titration simulation more flexible and closer to the real experiments. Thus , I
expect more repeatable and accurate data will come from this method. But there are still some
minor problems. The authors assume the van der Waals energy changes linearly or
quadratically during the change of reaction coordinates, which may not be true compared with
that in the real world. Besides, the rate of the changes of van der Waals energy is coherent with
that of electrostatic energy in this paper. It is possible that assigning different reaction
coordinates to them may produce better results, just like what is done in the first paper.

References:
Kong and Brooks III.1996. λ−dynamics: A new approach to free energy
calculations. J.Chem. Phys. 105:2414-2423

Khandogin and Brooks III. 2005. Constant pH Molecular Dynamics with Proton
Tautermerism. Biophysical Journal. 89:141-157

Challenges in Theoretical Chemistry - Intermolecular Potentials (Stone, Science, 2008)

StevenWang — Thu, 16 Apr 2009 02:49:00 GMT

Main points: This paper concerns about the basics, as well as the capabilities and limitations of the intermolecular potential calculation. As a novice, I like the basic part better. Because the small molecules cases are easy to cope with, the author cares more about the intermolecular potential calculation between two large molecules. In this situation, the author claims, intermolecular potential calculation is transformed into pairwise atom-atom potential calculation, which is the sum of all the potentials between any of the atoms in molecule A and another one in molecule B. This can cause a double count problem but the author does not mention too much about it. The calculation of the potential is through a general-form function: U {sub a & b}, whose variables are interatomic distance R and relative orientation Omega between the two atoms involved. For the Omega (relative orientation), I think this the result of concerning the electron distribution in the electron cloud of each atoms or the partial charge distribution, which is discussed later in the paper.

There are a total of five terms contributing to the potential function: 1) the electrostatic term; 2) the exchange-repulsion term or van de Waals repulsion; 3) the dispersion term or van de Waals attraction; 4) induction term; 5) charge transfer term. (1) For the electrostatic term, the author reviewed several methods that involve charges, dipoles, multipoles, etc, instead of using point charges. And, truly, these methods can interpret the electrostatic term better, especially in the hydrogen-bonded systems. (2) For the exchange-repulsion term, according to the author, the Born-Mayer form works better than the R {superscript -12} term in the Lennard-Jones potential. (3) As for the dispersion term or van de Waals attraction, it is calculated through a series with a introduction of a correction term (a damping function) for R near zero. (3) The induction term and charge transfer term are the most troublesome ones, from the perspective of the author. The induction term comes from the interaction between an electric-field-giving atom and the consequent induced dipole of another atom. The charge transfer term is caused by the donating of electron density of one atom to the acceptor molecular orbital (MO) of another atom (maybe I should say the MO of another molecule), which will give rise to charges in electron configuration. There are two ways of dealing with this induction term: inexplicit induction, by modifying the properties of each molecule and explicit induction, by assigning isotropic dipole-dipole polarizabilities to atoms and repeating solving induction term at each step during the optimization or simulation.

Up until now, it is already practical to accurately calculate the interatomic potentials for small molecules (using, say, ab nitio method) and molecules made up of 20-30 atoms (using, say, symmetry-adapted perturbation theory based on DFT). For macromolecules, the simplest way is employing atomic point charges, Lennard-Jones or Born-Mayer repulsions, and isotropic atom-atom dispersion, while other treatments are also available. However, there is currently no definitive all-purpose force field.

Two limitations are associated with intermolecular potential calculation. One is the difficulty in describing flexible molecules. The other is the limited available experimental data for the use of parametrizing the force field.

Generally, this is a good review of the development of intermolecular potential calculation and it is appropriate for all level reader interested in this field.

References:

Stone, Science, 321, 787-789, 2008.

Replica-Exchange MD (Sugita and Okamoto, Chem Phys Lett, 1999; Nymeyer et al, Methods Enzymol, 2004)

Kevin — Tue, 14 Apr 2009 22:22:00 GMT

The purpose of replica-exchange molecular dynamics is to provide adequate sampling of a wide conformational space MD simulation. In MD simulations, especially for biomolecules, the energy landscape often contains large energy barriers that will trap the system in saddle points, or local minima. In order to find the absolute minimum of the system we require a method that will allow us to cross over the large energy barriers in search of deeper minima. This is where replica-exchange helps us. In replica-exchange we create several copies of a system that are identical except for their temperature. The higher temperature replicas will be better able to sample over the entire conformational space and not be trapped as easily in local minima. Each replica is searching for energy minima. After a set number of steps in the simulation the energy of each replica is compared with its nearest neighbor(s). Exchanges of temperature are then determined according to the following probability: min{1,exp(ΔH*Δβ). Thus, if a higher temperature replica is at a lower energy than the lower temperature replica, the temperatures will be exchanged and the simulation continues. If the energy is greater there is still a chance, although small, that the temperature will be exchanged. The basic rundown is this: if a high temperature replica finds a deep energy well, the temperature will be lowered for the newly found well repeatedly until the absolute minimum is reached. If an even deeper energy well in the landscape is discovered by a higher temperature replica, the process repeats itself focusing on the new energy well. This continues until, ideally, the absolute minimum of the energy landscape is found.

References:

Sugita and Okamoto, Chem Phys Lett, 1999
Nymeyer et al, Methods Enzymol, 2004

Challeges in Theoretical Chemistry - Surface Scattering Simulations (Kroes, Science, 2008)

Kevin — Tue, 14 Apr 2009 19:10:00 GMT

The goal of surface scattering simulations is to accurately model how molecules interact with surfaces. This is a daunting challenge because molecule-surface interactions can result in many different processes including: dissociation, scattering, and various excitations. Currently most molecule-surface interactions are modeled using the Born-Oppenheimer approximation which limits the reaction to one potential energy surface (PES). However, the BO approximation does not generally hold for non-adiabatic processes. Recent research has studied whether or not there are certain reactions that can be accurately described using the BO approximation. In general the BO approximation has proven fairly accurate for light molecules (i.e. H₂) but heavier molecules tend to experience other surface-interactions that affect its scattering such as phonon excitation and electron-hole excitation. These scenarios ideally should be modeled using a diabatic PES or multiple PESs which is not the case with the Born-Oppenheimer. Many challenges remain in this field including: developing an accurate method to model interactions in the electronic ground state, excited electronic states, as well as accurately describing the effects of phonon excitation. Overall the article does not discuss surface scattering as much as I expected. Rather it focused heavily on modeling phonon and electron-hole excitations. These do present a challenge to calculating accurate scattering simulations, but I expected more discussion on how scattering interactions (vibrational, rotational, and diffractive) are calculated.

Reference: Kroes, Science, 321, 794-797, 2008

Challenges in Theoretical Chemistry - Quantum Dynamic of Chemical Reactions (Clary, Science, 2008)

Linh — Fri, 10 Apr 2009 21:25:00 GMT

Quantum Dynamic (QD) is used to study chemical reactions in the gas phase.Current computational and theoretical methods are used to calculate the rate of reactions and to predict the detail of the reactions.

Methods: For reactions with high energy barrier such as those involving bond breaking or formation, QD methods need to be used.

QD procedure: Solve Schrodinger equation to get the electronic structure. Compute Potential Energy Surface (PES) using ab initio electronic structure methods. This yield reliable reaction barriers. Use these PES in the solutions of Schrodinger equation for the nuclei to predict fundamental informations of chemical reactions.

For three atoms reactions, a computer code called ABC, which involves hyperspherical coordinates, is used. This theory can be extended to full dimensions by including functions that describe molecular rotations. Agreement between this theory and experiment is remarkable. For reactions with several maxima and minima in PES, "direct dynamic" methods are sometimes used. In these methods, electronic energies are calculated directly from an ab initio code whenever they are required. However direct dynamic are difficult to aplly in QD calculations due to large number of quantum wave functions. An alternative method is used which involves solving the time dependent Schrodinger equation. This is called the "wave packet" calculations.

For four or more atoms reactions, a method is developed that allows the construction of effective Potential Surface (PS) in reduced dimensionality through a small number of accurate electronic structure calculations. These PS are used in time independent quantum scattering calculations with hyperspherical coordinates, treating only bonds broken and formed in reactions explicitly. This method is used in reaction between H atom and propane which has result agree swell with experiment. Another approach is multiconfigurational time dependent Hartree (MCTDH) which treat quantum wave functions as a product of single particles and time dependent functions. MCTDH can be used with a method for calculating the quantum fluxes in a chemical reaction to calculate the reactive rate constant of multidimensional system.

Challenges in Theoretical Chemistry - Modeling Materilals Properties (Carter, Science, 2008)

Jason — Tue, 17 Feb 2009 15:52:00 GMT

This perspective article in Science Magazine outlines the goals and challenges in attempting to model materials properties from first principles. The article is a concise review of the specialized techniques that are most successful in modeling different types of properties. Reading this article would be a good first step to find the method of choice to model a particular property of a material. Of course, as a review, the article does not discuss how the different techniques are actually implemented or why one technique is better suited for a particular problem. You must go to the references to get a deeper level of understanding for each technique mentioned.

Reference: Carter, Science, 321, 800-803, 2008.

Hydrophobic Effect as Investigated by CPMD

Jason — Mon, 15 Dec 2008 23:33:00 GMT

The hydrophobic effect has traditionally been thought to be caused by energetic and entropic factors. A recent ab initio molecular mechanics simulation gave some surprising results that seems to suggest that the entropic component is not as powerful as one may think. The entropic component of the HE is based on the idea that around the hydrophobic entity a 'frozen' clathrate cage forms; when the hydrophobic regions aggregrate some water is freed from this order cage and causes a free energy benefit. A recent car-parrinello simulation showed that there is a floppy solvation shell around the smallest possible hydrophobic molecule, a hydrogen atom, that promotes diffusion. If this behavior is assumed to be the same for larger hydrophobic molecules, this means that the hydrophobic effect is solely an energetic phenomenon. However, there may well be dramatically different behavior of water surrounding larger hydrophobic molecules, and this should be investigated further. ( Barbara Kirchner, John Stubbs, and Dominik Marx, Phys. Rev. Lett., 89: 215901 2002)

Nelson, P., Biological Physics (2004)

Jason — Tue, 09 Dec 2008 03:12:00 GMT

Nelson's book covers the basics of biological physics. The book reviews fundamental topics such as thermodynamics and makes the connection to biological processes. The book is a good survey of biological physics and discusses recent research areas such as single molecule experiments. The book gives many examples and thorough explanation. It is appropriate for anyone who is beginning to learn about biological physics.

Pathria R.K., Statistical Mechanics 2nd Ed. (1996)

Jason — Tue, 09 Dec 2008 02:49:00 GMT

Pathria's book covers all the essential material related to equilibrium statistical mechanics. The material is covered in depth but, much of the learning must be done by working the numerous problems found at the end of each chapter. This book is great for anyone who wants to learn the tools of statistical mechanics. The book is concise but gives adequate explanation.