1、Ch121a Atomic Level Simulations of Materials and Molecules,William A. Goddard III, wagwag.caltech.edu Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology,BI 115 Hours: 2:30-3:30 Monday and Wednesday Lecture or Lab: Friday 2-3pm (
2、+3-4pm),Teaching Assistants Wei-Guang Liu, Fan Lu, Jose Mendoza, Andrea Kirkpatrick,Lecture 9, April 18, 2011 Monte Carlo Polymer distributions,Monte Carlo Method and applications in chemistry, biological, polymer related fields,Youyong Li and William A. Goddard III,Continuous self-avoiding walk wit
3、h application to the description of polymer chains; Youyong Li and WAG; J. Phys. Chem. B, 110: 18134-18137 (2006),Efficient Monte Carlo Method for Free Energy Evaluation of Polymer Chains; Jiro Sadanobu and WAG; Fluid Phase Equilibria 144, 415 (1998),The Continuous Configurational Boltzmann Biased D
4、irect Monte Carlo Method for Free Energy Properties of Polymer Chains; Jiro Sadanobu and WAG; J. Chem. Phys. 106, 6722 (1997),references,The Monte Carlo method provides approximate solutions to a variety of mathematical problems by performing statistical sampling experiments (on a computer). The met
5、hod applies to problems with no probabilistic content as well as to those with inherent probabilistic structure.,What is Monte Carlo?,Partition the space into grids for numerical evaluation Tedious, thorough study: enumerate all N2 grid points and evaluate whether satisfies rR Monte Carlo study: sam
6、ple part of the grids (e.g. 10%) to get the conclusion. For interesting problems, tedious and thorough study is seldom practical.,Simple example: Evaluate p from the area of a circle from numerical evaluation (dartboard method),Named after the city in the Monaco principality of France, famous for it
7、s gambling casinos. Roulette is a simple random number generator. The name and the systematic development of Monte Carlo methods dates from LASL 1944.,In 1873: A. Hall performed Needle-board experiment to determine p ( “On an experimental determination of Pi”).,In 1899, Lord Rayleigh showed that a o
8、ne-dimensional random walk without absorbing barriers could provide an approximate solution to a parabolic differential equation.,In 1931 Kolmogorov showed the relationship between Markov stochastic processes and certain integro-differential equations.,In 1948, Harris and Herman Kahn systematically
9、developed the ideas used in the study of random neutron diffusion in fissile material of the atomic bomb., 1948 Fermi, Metropolis, and Ulam obtained Monte Carlo estimates for the eigenvalues of Schrodinger equation,History of Monte Carlo method,Major components of Monte Carlo Algorithm,Probability d
10、istribution functions (pdfs) - the physical (or mathematical) system must be described by a set of pdfs. Random number generator - a source of random numbers uniformly distributed on the unit interval must be available. Sampling rule - a prescription for sampling from the specified pdfs, assuming th
11、e availability of random numbers on the unit interval, must be given. Scoring (or tallying) - the outcomes must be accumulated into overall tallies or scores for the quantities of interest. Error estimation - an estimate of the statistical error (variance) as a function of the number of trials and o
12、ther quantities must be determined.,Refinements,Variance reduction techniques - methods for reducing the variance in the estimated solution to reduce the computational time for Monte Carlo simulation Parallelization and vectorization - algorithms to allow Monte Carlo methods to be implemented effici
13、ently on advanced computer architectures.,Two main reasons: Global minimization (avoiding being trapped) Ensemble description: Ergodicity (avoiding being trapped),MD: straightforward algorithm, but only samples the high frequency region efficiently, may not be ergotic,MC: Complicated algorithm, but
14、much more practical to realize ergodicity in difficult problems,Why we need MC or MD?,Side chain torsions- Number of side chain torsions vary for each aa,Main chain torsions (peptide bond is planar) 180o,Backbone and Side chain torsions in Proteins - Flexible and complicated system,Ramachandran plot
15、,If we sample 10 points for each torsion of 40 monomers or residue, there will be 10400 possible conformations!,Systems with Many Degrees of Freedom,The complexity of polymer, protein, and other macromolecules,Markov Processes,Let us set up a so-called Markov chain of configurations Ct by the introd
16、uction of a fictitious dynamics. The time t is computer time (marking the number of iterations of the procedure), NOT real time - our statistical system is considered to be in equilibrium, and thus time invariant. Let P(A,t) be the probability of being in configuration A at time t.Let W(A-B) be the
17、probability per unit time, or transition probability, of going from A to B. Then:,At large t, once the arbitrary initial configuration is forgotten, want,Detailed balance,Clearly a sufficient (but not necessary) condition for an equilibrium (time independent) probability distribution is the so-calle
18、d detailed balance condition,This method can be used for any probability distribution, but if we choose the Boltzmann distribution,Note: Z does not appear in this expression! It only involves quantities that we know or can easily calculate,The Metropolis Algorithm,This dynamic method of generating a
19、n arbitrary probability distribution was invented by Metropolis, Teller , and Rosenbluth in 1953 (supposedly at a Los Alamos dinner party).,There are many possible choices of the Ws which will satisfy detailed balance. They chose a very simple one:,So if,And if,Realize Metropolis Algorithm,So we hav
20、e a valid Monte Carlo algorithm if: We have a means of generating a new configuration B from a previous configuration A such that the transition probability W(A-B) satisfies detailed balance The generation procedure is ergodic, i.e. every configuration can be reached from every other configuration i
21、n a finite number of iterations The Metropolis algorithm satisfies the first criterion for all statistical systems. The second criterion is model dependent, and not always true (e.g. at T=0).,Folding Kinetics,New View: Energy Bias,Simple sampling,Independent Rotational Sampling,Monte Carlo Sampling
22、for protein structure,The probability of finding a protein (biomolecule) with a total energy E(X) is:,Estimates of any average quantity of the form:,using uniform sampling would therefore be extremely inefficient.,Metropolis et al. developed a method for directly sampling according to the actual dis
23、tribution.,Metropolis et al. Equation of state calculations by fast computing machines. J. Chem. Phys. 21:1087-1092 (1953),Monte Carlo for the canonical ensemble,The canonical ensemble corresponds to constant NVT. The total energy (Hamiltonian) is the sum of the kinetic energy and potential energy:
24、E=Ek(p)+Ep(X) If the quantity to be measured is velocity independent, it is enough to consider the potential energy:,The kinetic energy depends on the momentum p; it can be factored and canceled.,Monte Carlo for the canonical ensemble,Let:,let be transition probability from state X to state Y.,Suppo
25、se we carry out a large number of Monte Carlo simulations, such that the number of points observed in conformation X is proportional to N(X). The transition probability must satisfy one obvious condition: it should not destroy this equilibrium once it is reached. Metropolis proposed to realize this
26、using the detailed balance condition:,or,Monte Carlo for the canonical ensemble,There are many choices for the transition probability that satisfy the balance condition. The choice of Metropolis is:,The Metropolis Monte Carlo algorithm: Select a conformation X at random. Compute its energy E(X) Gene
27、rate a new trial conformation Y. Compute its energy E(Y) Accept the move from X to Y with probability:,Pick a random number RN, uniform in 0,1. If RN P, accept the move.,4. Repeat 2 and 3.,Monte Carlo for the canonical ensemble,Notes: There are many types of Metropolis Monte Carlo simulations, chara
28、cterized by the generation of the trial conformation. The random number generator is crucial Metropolis Monte Carlo simulations are used for finding thermodynamics quantities, for optimization, An extension of the Metropolis algorithm is often used for minimization: the simulated annealing technique
29、, where the temperature is lowered as the simulation evolves, in an attempt to locate the global minimum.,Grand Canonical Monte Carlo,A classical force-field method where molecules interact and are .,. such that the proper distribution of positions, energies, and numbers of molecules is achieved for
30、 a system at fixed chemical potential (or fugacity or pressure) and temperature. GCMC is often used to study the adsorption of small molecules in inorganic materials such as zeolites.,Molecular Size via Computational Sieving,Use a Force-Field and Grand Canonical Monte Carlo to adsorb a gas molecule
31、into a confined system at a standard T & P. If significant numbers of the molecule fit into the system, the characteristic size of the system is related to the size of the molecule. An unbiased method which uses the informationinherent in the FF.,D = Dslit - d,Static Monte Carlo: CCBB,SS-DMC (Simple
32、 sampling direct Monte Carlo),IRS (Independent rotational sampling),fixed bond lengths and angles, random sampling of all torsions,Drastic sampling attrition,Rotational biased sampling, using the normalized torsion weighting function (TWF), WIRS:,NonBond bias corrected TWF,Excludes high torsion ener
33、gies Still have spatial overlaps,W* updated every step Too expensive for large system,CCB (Continuous configurationally biased direct Monte Carlo),Off-lattice, Rc introduced into TWF,CCB (Continuous configurationally biased direct Monte Carlo),United atomistic FF Isolated chain,Sadanobu, J. and Godd
34、ard, W. A., III. J. Chem. Phys. 1997, 106, 6722,CCBTX (Continuous configurationally biased TX direct Monte Carlo),Use CCB works for isolated UA chain, also for full atomistic model, multiple chains, dense polymer/ dendrimer,Description of CCBTX method,We developed the efficient Continuous Configurat
35、ional Biased TX (CCBTX) method to generate high quality amorphous polymer and dendrimer atomistic structures directly. The method developed here can also be used for protein fold prediction. The code is implemented in python environment and CMDF package, which provides friendly interface.,CCBTX Mont
36、e Carlo method,Application of CCBB,Thermodynamic Functions and Critical Exponents for Polymer Chains from CCBB,1. Polymer chain dimensions,Flory theta temperature: 2n=1,The radius of gyration, ln, as a function of polymer length, ln(N-1) for a range of temperatures (CCBB calculations using SKSFF). T
37、he uncertainties are shown, but are less than the size of the symbols.,Where n = 1/2,This confirms that ln is proportional to ln(N-1), as in following equation. Here the slope, 2n, is a function of temperature,Free energy theta temperature: g=1,Application of CCBB: Thermodynamic functions of polymer
38、s,Old methods of MC assumed atoms on a lattice. This gives bad estimates of entropy CCBB allows a continuous description (no lattice),Standard deviation s among 10 independent samples of CCBB calculations as a function of the number of clusters sampled.,Free energy A/kT for polymer chains of length
39、N at different temperatures from CCBB calculations with SKS-FF. The uncertainties are shown, but are less than the size of the symbols.,Linearity is expected from mean field description.,We see that A/kT is proportional to N, where the slope gives -lnm,Internal energy E/kT for polymer chains of leng
40、th N at different temperatures from CCBB calculations with SKS-FF. The uncertainties are shown, but are less than the size of the symbols. We find that E/kT is proportional to N, where the slope represents e0,Entropy S/k for polymer chains of length N at different temperatures from CCBB calculations
41、 with SKSFF. The uncertainties are shown, but are less than the size of the symbols. We find that S/k is proportional to N, where the slope is s0.,The radius of gyration, ln, as a function of polymer length, ln(N-1) for a range of temperatures (CCBB calculations using SKSFF). The uncertainties are s
42、hown, but are less than the size of the symbols.,This confirms that ln is proportional to ln(N-1), as in following equation. Here the slope, 2n, is a function of temperature,Critical attrition as a function of 1/T from CCBB with various force fields. The uncertainties are shown but are less than the
43、 size of the symbols.,We find that increases as the temperature (except for NoTorSKSFF) reflecting the increasing number of available conformations and partition function.,e0 as a function of 1/T from CCBB with different force fields. The uncertainties are shown, but are less than the size of the sy
44、mbols. We find that e0 decreases as temperature increases, except for NoTorSKSFF.,As the temperature increases, average energy increases slower than kT, except for NoTorSKSFF.,Radius of gyration, as a function of temperature from CCBB MC with various force fields, all for a chain length of N=300. Wh
45、en only the repulsive van der Waals energy is present (NoTorNoAttSKSFF), the polymer chains shrink with increasing temperature, which reflects the increasing number of available conformations.,When the attractive van der Waals energy is present (the other five force fields), the polymer chains colla
46、pse at low temperature due to the favorable enthalpy.,CR, Cn, CA derived from e0 at different temperatures from CCBB with NoTorNoAttSKSFF and NoTorSKSFF. We find that CR, Cn, CA is linear with 1/T, indicating that ER, En, EA are independent of temperature.,The critical exponent 2 as a function of 1/
47、T from CCBB with various force fields. As temperature decreases, the polymer chains collapse and 2n decreases except for NoTorNoAttSKSFF.,The result of NoTorNoAttSKSFF shows that including only self-avoidance leads to 2n 1.2, as deduced by Flory.,The Free energy increment from CCBB calculations (SKS
48、FF) for adding m=100 monomers A(N-m)-AN as a function of length N (plotted as lnN/(N-m). A(N-m)-AN/kT is proportional to lnN/(N-m), where the slope is (g-1).,Free energy theta temperature,The internal energy increment from CCBB calculations (SKSFF) for adding m=100 monomers E(N-m)-EN as a function o
49、f length N (plotted as lnN/(N-m). We see that E(N-m)-EN/kT is independent to lnN/(N-m) at high temperature limit, while dependent at low temperature.,At low temperature, the chain collapses and the average energy in the long chain is less than in the short chain.,The free energy critical exponent ( - 1) as a function of 1/T from CCBB MC using various force fields. increases with the temperature, reflecting the increased number of available conformations and partition function.,
