几种常见的优化方法ppt课件.ppt

1几种常见的优化方法电子结构几何机构函数稳定点最小点Taylor 展开：V(x) = V(xk) + (x-xk)V(xk) +1/2 (x-xk)2 V(xk)+. 当x是3N个变量的时候， V(xk)成为3Nx1的向量，而V(xk)成为3Nx3N的矩阵，矩阵元如：jixxVHessian21.一阶梯度法 a. Steepest descendentSk = -gk/|gk|directiongradient知道了方向，如何确定步长呢？最常用的是先选择任意步长l，然后在计算中调节kkkkSXXl1用体系的能量作为外界衡量标准，能量升高了则逐步减小步长。robust, but slow最速下降法 3最陡下降法（SD） 4b. Conjugate Gradient （CG) 共轭梯度1kkkkvgvl第k步的方向11kkkkkggggl标量Usually more efficient than SD, also robust不需要外界能量等作为衡量量利用了上一步的信息52。二阶梯度方法这类方法很多，最简单的称为Newton-Raphson方法，而最常用的是Quasi-Newton方法。Newtons method for finding an extreme point isxk+1 = xk - H-1(xk) y(xk)Quasi-Newton方法： use an approximation of the inverse Hessian.Form of approximation differs among methods牛顿-拉夫逊法 BFGS methodBroyden-Fletcher-Golfarb-ShannoDFP methodDavidon-Fletcher-Powell6Molecular dynamics 分子动力学分子动力学HistoryIt was not until 1964 that MD was used to study a realistic molecular system, in which the atoms interacted via a Lennard-Jones potential.After this point, MD techniques developed rapidly toencompass diatomic species, water (which is still thesubject of current research today!), small rigid molecules,flexible hydrocarbons and now even macromoleculessuch as proteins and DNA.These are all examples of continuous dynamicalsimulations, and the way in which the atomic motion iscalculated is quite different from that in impulsivesimulations containing hard-core repulsions.7What can we do with MD Calculate equilibrium configurational properties in a similarfashion to MC. Study transport properties (e.g. mean-squared displacement anddiffusion coefficients). MD in the NVT, NpT and NpH ensembles The united atom approximation Constraint dynamics and SHAKE Rigid body dynamics Multiple time step algorithmsExtend the basic MD algorithm8Impulsive molecular dynamics1.Dynamics of perfectly hard particles can be solved exactly, but process becomes involved for many part (N-body problem). 2.Can use a numerical scheme that advances the system forward in time until a collision occurs. 3.Velocities of colliding particles (usually a pair!) then recalculated and system put into motion again. 4. Simulation proceeds by fits and starts, with a mean time between collisions related to the average kinetic energy of the particles. 5.Potentially very efficient algorithm, but collisions between particles of complex shape are not easy to solve, and cannot be generalised to continuous potentials. 9Continuous time molecular dynamics1. By calculating the derivative of a macromolecular forcefield, we can find the forces on each atomas a function of its position.2. Require a method of evolving the positions of the particles in space and time to produce a true dynamical trajectory.3. Standard technique is to solve Newtons equations ofmotion numerically, using some finite difference scheme,which is known as integration.4. This means that we advance the system by some smalltime step t, recalculate the forces and velocities, and thenrepeat the process iteratively.5. Provided t is small enough, this produces an acceptable approximate solution to the continuous equations of motion.10Example of integrator for MD simulation One of the most popular and widely used integrators isthe Verlet leapfrog method: positions and velocities ofparticles are successively leap-frogged over each otherusing accelerations calculated from force field. The Verlet scheme has the advantage of high precision(of order t4), which means that a longer time step canbe used for a given level of fluctuations. The method also enjoys very low drift, provided anappropriate time step and force cut-off are used.r(t+D Dt)=r(t)+v(t+D Dt/2)D Dtv(t+D Dt/2)=v(t-D Dt/2)+a(t+D Dt/2)D Dt11Other integrators for MD simulations Although the Verlet leapfrog method is not particularlyfast, this is relatively unimportant because the timerequired for integration is usually trivial in comparison tothe time required for the force calculations. The most important concern for an integrator is that itexhibits low drift, i.e. that the total energy fluctuates aboutsome constant value. A necessary (but not sufficient)condition for this is that it is symplectic. Crudely speaking, this means that it should be timereversible (like Newtons equations), i.e. if we reverse themomenta of all particles at a given instant, the systemshould trace back along its previous trajectory.12Other integrators for MD simulations The Verlet method is symplectic, but methods such aspredictor-corrector schemes are not. Non-symplectic methods generally have problems withlong term energy conservation. Having achieved low drift, would also like the energyfluctuations for a given time step to be as low as possible. Always desirable to use the largest time step possible. In general, the trajectories produced by integration willdiverge exponentially from their true continuous paths dueto the Lyapunov instability. However, this does not concern us greatly, as the thermalsampling is unaffected expectation values unchanged.13Choosing the correct time step 1. The choice of time step is crucial: too short and phase space is sampled inefficiently, too long and the energy will fluctuate wildly and the simulation may become catastrophically unstable (“blow up”). 2. The instabilities are caused by the motion of atoms being extrapolated into regions where the potential energy is prohibitively high (e.g. atoms overlapping). 3. A good rule of thumb is that when simulating an atomic fluid, the time step should be comparable to the mean time between collisions (about 5 fs for Ar at 298K). 4. For flexible molecules, the time step should be an order of magnitude less than the period of the fastest motion (usually bond stretching: CH around 10 fs so use 1 fs).14For classic MD, there could be many tricks to speed up calculations, all centering around reducing the effort involved in the calculation of the interatomic forces, as this is generally much more time-consuming than integration.For example1. Truncate the long-range forces: charge-charge, charge-dipole2. Look-up tablesFor first principles MD, as forces are evaluated from quantum mechanics, we are only concerned with the time-step.15Because the interactions are completely elastic andpairwise acting, both energy and momentum areconserved. Therefore, MD naturally samples from themicrocanonical or NVE ensemble.As mentioned previously, the NVE ensemble is not very useful for studying real systems. We would like to be able to simulate systems at constant temperature or constant pressure.The simplest MD, like verlet method, is a deterministic simulation technique for evolving systems to equilibrium by solving Newtons lawsnumerically.16MD in different thermodynamic ensembles In this lecture, we will discuss ways of using MD tosample from different thermodynamic ensembles, whichare identified by their conserved quantities. Canonical (NVT) Fixed number of particles, total volume and temperature.Requires the particles to interact with a thermostat. Isobaric-isothermal (NpT) Fixed number of particles, pressure and temperature. Requiresparticles to interact with a thermostat and barostat. Isobaric-isenthalpic (NpH) Fixed number of particles, pressure and enthalpy. Unusual, butrequires particles to interact with a barostat only.17Advanced applications of MD We will then study some more advanced MD methodsthat are designed specifically to speed up, or makepossible, the simulation of large scale macromolecularsystems. All these methods share a common principle: they freezeout, or decouple, the high frequency degrees of freedom. This enables the use of a larger time step withoutnumerical instability. These methods include: United atom approximation Constraint dynamics and SHAKE Rigid body dynamics Multiple time step algorithms18Revision of NVE MD Lets start by revising how to do NVE MD. Recall that we calculated the forces on all atoms from thederivative of the force field, then integrated the e.o.m.using a finite difference scheme with some time step t. We then recalculated the forces on the atoms, andrepeated the process to generate a dynamical trajectoryin the NVE ensemble. Because the mean kinetic energy is constant, the averagekinetic temperature TK is also constant. However, in thermal equilibrium, we know thatinstantaneous TK will fluctuate. If we want to sample fromthe NVT ensemble, we should keep the statisticaltemperature constant.19Extended Lagrangians There are essentially two ways to keep the statisticaltemperature constant, and therefore sample from the trueNVT ensemble. Stochastically, using hybrid MC/MD methods Dynamically, via an extended LagrangianWe will describe the latter method in this lecture An extended Lagrangian is simply a way of including adegree of freedom which represents the reservoir, andthen carrying out a simulation on this extended system. Energy can flow dynamically back and forth from thereservoir, which has a certain thermal inertia associatedwith it. All we have to do is add some terms to Newtonsequations of motion for the system.20Extended Lagrangians The standard Lagrangian L is written as the difference ofthe kinetic and potential energies: Newtons laws then follow by substituting this into theEuler-Lagrange equation: Newtons equations and Lagrangian formalism areequivalent, but the latter uses generalised coordinates.VxmLNiii221.0iixLxLdtd.21Canonical MD So, our extended Lagrangian includes an extracoordinate , which is a frictional coefficient that evolves in time so as to minimise the difference between the instantaneous kinetic and statistical temperatures. The modified equations of motion are:The conserved quantity is the Helmholtz free energy.1/ )(1pFp/pr2SKTiiiiiiTtTm(modified form of Newton II)22Canonical MD By adjusting the thermostat relaxation time T (usuallyin the range 0.5 to 2 ps) the simulation will reach anequilibrium state with constant statistical temperature TS. TS is now a parameter of our system, as opposed to themeasured instantaneous value of TK which fluctuatesaccording to the amount of thermal energy in the systemat any particular time. Too high a value of T and energy will flow very slowlybetween the system and the reservoir (overdamped). Too low a value of T and temperature will oscillate aboutits equilibrium value (underdamped). This is the Nos-Hoover thermostat method.23Canonical MD There are many other methods for achieving constanttemperature, but not all of them sample from the trueNVT ensemble due to a lack of microscopic reversibility.We call these pseudo-NVT methods, and they include: Berendsen methodVelocities are rescaled deterministically after each step so that the system is forced towards the desired temperature Gaussian constraintsMakes the kinetic energy a constant of the motion by minimising the least squares difference between the Newtonian and constrained trajectories These methods are often faster, but only converge on thetrue canonical average properties as O(1/N).24Isothermal-isobaric MD We can apply the extended Lagrangian approach to simulations at constant pressure by simply adding yet another coordinate to our system. We use , which is a frictional coefficient that evolves intime to minimise the difference between the instantaneous pressure p(t), measured by a virial expression, and the pressure of an external reservoir pext. The equations of motion for the system can then be obtained by substituting the modified Lagrangian into the Euler-Lagrange equations. These now include two relaxation times: one for thethermostat T, and one for the barostat p.25Isothermal-isobaric MDThe is known as the Nos-Hoover method (Melchionnatype) and the equations of motion are:26Isothermal-isobaric MD27Constraint dynamicsfreeze the bond stretching motions of the hydrogens (or any other bond, in principle). We apply a set of holonomic constraints to the system, which are just algebraic equations connecting the coordinates of the individual atoms. These constraints contain information about the bond lengths. The most widely used algorithm for implementing these constraints is called SHAKE. This works on the principle of applying fictitious forces to all the particles such that they follow appropriately constrained trajectories. There is also a variant of SHAKE called RATTLE.28We can calculate what the constraint forces are (obviously they must be directed parallel to the bonds) by introducing Lagrange multipliers lij:These multipliers can be determined by substituting the modified expressions for Newton II into our Verlet integration scheme and imposing the constraints that:29 This leads to a set of simultaneous equations for lij whichcan be recast into a matrix form and solved by invertingthe matrix. For larger molecules, this is very time consuming and soan iterative procedure is used (hence the “shaking” of themolecule as each constraint is enforced in turn). Angle and higher order constraints can easily beincorporated into the algorithm by enforcing fictitiousconstraint bonds between atoms such that the anglebetween them is held constant. However, there are problems with convergence of thealgorithm if the molecule is overconstrained.