第十章分子模拟
第十章-中枢神经系统脱髓鞘疾病
第十章中枢神经系统脱髓鞘疾病(Demyelinating Diseases of theCentral Nervous System) 第一节概述脱髓鞘疾病(demvelinative diseases) 是一组脑和脊髓以髓鞘破坏或脱髓鞘病变为主要特征的疾病，脱髓鞘是其病理过程中具有特征性的突出表现。
脱髓鞘疾病通常公认的病理标准是：①神经纤维髓鞘破坏，呈多发性小的播散性病灶，或由一个或多个病灶融合而成的较大病灶；②脱髓鞘病损分布于中枢神经系统(CNS)白质，沿小静脉周围的炎症细胞浸润；③神经细胞、轴突及支持组织保持相对完整，无华勒变性或继发传导束变性。
因本组疾病并非采用病因学分类，故不完全符合上述标准，如Schilder 病和坏死性出血性白质脑炎的轴突损害几乎与髓鞘同样严重，仍在本组讨论；反之，某些疾病脱髓鞘病损较突出，但因病因已经清楚而未到A本组疾病，如缺氧性脑病由于慢性缺氧，Bin ·swanger病因慢性缺血，伴恶性贫血的脊髓亚急性联合变性系因维生素Bt2 缺乏，热带痉挛性截瘫(TsP) 因逆转录病毒所致，进行性多灶性白质脑病(PML) 是发生在免疫缺陷患者的少突胶质细胞病毒感染。
第二节多发性硬化多发性硬化( ㈣Itipl sclerosis ，MS)是以中枢神经系统白质脱髓鞘病变为特点的自身免疫病，可能是遗传易感个体与环境因素作用而发生的自身免疫过程。
由于其发病率较高、呈慢性病程和倾向于年轻人罹患，而成为最重要的神经系统疾病之一o 【病因学及发病机制】1．病毒感染与自身免疫反应MS的确切病因及发病机制迄今不明。
许多流行病学资料均提示，MS与儿童期接触的某种环境因素有关，推测这种因素可能是病毒性感染，但尚未从MS 患者的脑组织中发现或分离出病毒，包括20 世纪60 年代曾高度怀疑的作为嗜神经病毒的麻疹病毒，以及20 世纪80 年代的逆转录病毒，ePA~qf T 淋巴细胞病毒I 型(human T-lymphotroplc virus typeI ，HTLV· I) 。
分子动力学模拟的基本步骤
分子动力学模拟的基本步骤
碳纳米管管壁上C-C 原子之间形成强的相互作用σ健, 而在垂直管壁的方向为弱的π健相互作用, 使得C-Si 之间相互作用很弱, 因而碳纳米管在低温下可以作为硅的惰性容器, 从而也可以作为合成硅纳米线的模板.
由于碳纳米管-硅纳米线复合材料具有更高的热稳定性, 预计将来在纳米电子学和微电子器件中将有潜在的应用价值.
分子动力学模拟的基本步骤
1．确定研究对象
2．分子的初始位置和速度
3．势能模型——分子间作用力
4．分子运动方程的建立
5．周期性边界条件
6．位能截断
对于分子数为N的模拟体系，原则上任何两个分子问都存在相互作用，那么在计算体系位能时须进行N(N-I)／2次运算，一般情况下要占总模拟时间的80％左右，非常消耗机时。
为提高计算效率，在实际模拟过程应进行势能截断，最为常用的方法是球形截断法，截断半径一般取2．6D（D为分子的直径)，这样对截断距离之外分子间的相互作用就可以忽略，模拟过程中减少了计算量。
7．实施模拟
在周期性边界条件、时间平均等效于系综平均等基本假设之上，
通过求解体系的运动方程组得到各粒子在不同时刻的位置和速度。
体系达到充分平衡后，再经过几千、几万甚至几十万步的运算，体系的一些热力学参量可以通过统计平均得出。
生物物理学中的分子动力学模拟
生物物理学中的分子动力学模拟生物物理学是生物学与物理学的交叉学科，旨在研究生物大分子的结构与功能。
分子动力学模拟是生物物理学中的重要工具，用于研究分子在不同环境下的动力学行为。
本文将介绍分子动力学模拟的基本概念、应用和未来发展方向。
一、分子动力学模拟基本概念分子动力学模拟是利用计算机模拟分子在经典牛顿力学下的运动轨迹的过程。
分子动力学模拟的基本思想是将分子看作一组球体，通过求解牛顿运动方程，模拟它们在空间中的运动轨迹。
在模拟过程中，通过设置合适的势函数来描述分子之间的相互作用。
势函数主要包括键能、库伦势、范德华力、电子偶极子相互作用等。
模拟过程中还需要考虑分子的初始构象、温度、压力和场强等因素的影响。
二、分子动力学模拟的应用分子动力学模拟在生物物理学中的应用非常广泛，以下是一些常见的应用：1. 蛋白质动力学模拟蛋白质是生命体系中最重要的大分子之一，其结构与功能密切相关。
通过蛋白质动力学模拟，可以研究蛋白质的构象变化、动态行为及其与其他分子之间的相互作用。
例如，在研究药物靶点时，可以通过模拟药物分子与靶点蛋白之间的相互作用，来预测药物的活性及其副作用。
2. 生物膜模拟生物膜是生物体内各种细胞和细胞器之间的界面结构，是细胞膜的重要组成部分。
生物膜由脂质分子和蛋白质构成，其特殊的物理化学特性使得模拟其行为是非常具有挑战性的。
通过模拟生物膜的形成和变化，可以研究生物分子在膜内的运动与相互作用，为疾病治疗等领域提供理论基础。
3. RNA模拟RNA与DNA一样都是核酸分子，但其在功能和结构上有着巨大的差异。
通过RNA分子的模拟，可以研究RNA的三维结构、相互作用和在转录和翻译过程中的生物学功能等方面的问题，为生物医药领域的研究提供重要支撑。
三、分子动力学模拟的未来发展方向分子动力学模拟的应用领域不断扩大，未来其在生物物理学领域的应用将更为广泛。
以下是一些未来的发展方向：1. 强化学习算法在模拟中的应用强化学习是一种机器学习方法，在分子动力学模拟中，可以将其应用于动力学过程的控制和优化。
分子模拟第九、十章
分 子 模 拟牛继南njn0516@2011.3第玖章限制性分子动力学虽然分子动力模拟的计算能力随着计算机速度日益剧增，但通常能计算只有数纳秒内上万个分子体系的变化，与实际的体系相距还很远，因此，改进计算方法为分子动力模拟中最重要的问题之一。
9.1 刚性线性分子的计算•执行分子动力学计算时，若系统的热力学性质或宏观性质为主要的研究对象，则可以忽略分子内部的相对运动。
•通常采用刚性分子计算法，即固定分子中原子的相对位置，移动时分子作为一个整体。
•通过这种方法可以增加计算的效率，并可采用较大的步长。
平动转动•分子质心的坐标为：其中n为原子个数，m a和分别为a原子的质量和其位置矢量，M为分子的总质量。
•对于线性分子，如HCN和CO2，其角速度和力矩必定垂直其分子轴，设为分子轴向的单位矢量，则：•其中取决于分子间作用力。
具有以下特性：可以依此来求分子力矩平方的平均值。
•分子中原子相对于质心的位置为：其中 为a原子相对于质心的距离。
•则：• 为原子a所受的力。
式子中的 仅有垂直于分子轴的分量 起作用，即：•对于 ，可由下式求出：•线性分子的转动可写成：其中为角速度，为分子的转动惯量。
为拉格朗日未定乘子（Lagrange undetermined multiplier）。
由于和，以LF法解此微分式可得:•分子的移动即为质心的移动，而质心所受的力为分子中原子所受力的总和，即：•由LF法解得的质心运动方程为：•式中M为分子的总质量。
9.2 刚性非线性分子的计算•对于非线性刚性分子，造成分子转动的力矩为：•分子的位向定义了分子固定轴和空间固定轴之间的关系。
一般分子固定轴的选取是使分子的转动惯量矩阵成对角化状态。
•设空间中一个单位向量在分子轴坐标系中的位置是 ，在空间坐标系中的位置是 ，两者之间的关系为：•式中，转置矩阵和欧拉角(Ruler angle) 相关。
两坐标系O-ξ η ζO-x y z交线ON称为节线ON-Oξφ进动角ON-Oxψ自转角Oξ -Oz θ 章动角•由欧拉角确定的矩阵 为：•得到分子坐标轴的运动方程为：•式中 为分子固定轴的3个转动惯量。
分子动力学模拟
分子动力学模拟分子动力学模拟是一种重要的计算方法，用来研究分子体系的运动和相互作用。
该方法基于牛顿力学和统计力学的原理，通过数值模拟来预测和描述分子在不同条件下的行为。
在分子动力学模拟中，通过计算每个分子的受力和相互作用，可以得到关于分子位置、速度和能量等物理量的时间演化。
这些信息可以被用来研究分子体系的动力学、热力学和结构性质等。
为了进行分子动力学模拟，需要确定分子的力场和初始状态。
力场是一组描述分子分子间相互作用的数学函数，包括键的强度、键角的刚度、电荷分布等。
初始状态则是给定分子的初始位置和速度。
在分子动力学模拟中，分子受到的力主要来自于势能函数的梯度。
通过运用牛顿运动方程，可以计算得到每个分子的加速度，并进一步更新位置和速度。
这个过程重复进行，直到达到所需的模拟时间。
分子动力学模拟可以用来研究各种不同类型的分子体系。
例如，可以模拟液体中分子的运动和结构，以研究其流变性质和相变行为。
还可以模拟气体中分子的运动和相互作用，以研究化学反应和传输过程。
此外，分子动力学模拟还可以用来研究固体材料的力学性质和热导率等。
通过模拟材料内部原子的动力学行为，可以计算材料的弹性模量、杨氏模量等力学性质。
同时，还可以计算材料的热导率，从而了解其热传导性能。
分子动力学模拟已经成为了许多领域的重要工具。
它在材料科学、生物科学、化学工程和环境科学等领域中都得到了广泛应用。
通过模拟和理解分子体系的行为，我们可以更好地设计新材料、药物和催化剂，以及解决各种科学和工程问题。
然而，分子动力学模拟也有一些局限性。
首先，模拟的时间尺度受到限制，通常只能模拟纳秒或微秒级别的时间。
其次，模拟的精度也受到一定的限制，特别是在处理量子效应和极化效应等方面。
为了克服这些限制，研究人员正在发展和改进分子动力学模拟的方法。
例如，开发更精确的势能函数和更高效的计算算法，可以提高模拟的时间尺度和精度。
同时，与实验相结合，通过验证和修正模型，也可以提高模拟的可靠性和预测能力。
分子动力学模拟方法
分子动力学模拟方法分子动力学模拟是一种用于研究分子系统在原子尺度上运动规律的计算方法。
通过模拟分子在一定时间范围内的运动轨迹，可以揭示分子在不同条件下的结构、动力学和热力学性质，为理解分子系统的行为提供重要信息。
本文将介绍分子动力学模拟的基本原理、常用方法和应用领域。
分子动力学模拟的基本原理是利用牛顿运动方程描述分子系统中原子的运动。
根据牛顿第二定律，分子系统中每个原子受到的力可以通过势能函数求得，从而得到原子的加速度，再通过数值积分方法求解原子的位置和速度随时间的演化。
通过大量的时间步长积分，可以得到分子系统在一段时间内的运动轨迹。
在实际应用中，分子动力学模拟可以采用不同的数值积分方法，如Verlet算法、Leap-Frog算法等。
这些算法在计算效率和数值稳定性上有所差异，根据模拟系统的特点和研究目的选择合适的数值积分方法至关重要。
此外，分子动力学模拟还需要考虑原子间相互作用的描述方法，如分子力场、量子力场等，以及边界条件和初值设定等参数的选择。
分子动力学模拟方法在材料科学、生物物理、化学反应动力学等领域有着广泛的应用。
在材料科学中，可以通过模拟材料的力学性能、热学性质等，为新材料的设计和开发提供参考。
在生物物理领域，可以研究蛋白质、核酸等生物大分子的结构和功能，揭示生物分子的运动规律和相互作用机制。
在化学反应动力学研究中，可以模拟分子在化学反应中的动力学过程，为理解反应机理和优化反应条件提供理论支持。
总之，分子动力学模拟方法是一种强大的研究工具，可以深入理解分子系统的运动规律和性质。
随着计算机硬件和软件的不断发展，分子动力学模拟在科学研究和工程应用中的地位将更加重要，为解决现实世界中的科学和工程问题提供重要的理论和技术支持。
通过本文的介绍，相信读者对分子动力学模拟方法有了更深入的了解。
希望本文可以为相关领域的研究工作提供一定的参考和帮助，促进分子动力学模拟方法在更多领域的应用和发展。
《分子模拟方法》课件
加速研发进程
分子模拟可以大大缩短药 物研发、材料合成等领域 的实验周期，降低研发成 本。
揭示微观机制
通过模拟，可以揭示分子 间的相互作用机制和反应 过程，有助于深入理解物 质的性质和行为。
分子模拟的发展历程
经典力学模拟
基于牛顿力学，适用于 较大分子体系，但精度
较低。
量子力学模拟
适用于小分子体系，精 度高，但计算量大，需
详细描述
利用分子模拟方法，模拟小分子药物与生物大分子（如蛋白质、核酸等）的相 互作用过程，探究药物的作用机制和药效，为新药研发提供理论支持。
高分子材料的模拟研究
总结词
研究高分子材料的结构和性能，优化 材料的设计和制备。
详细描述
通过模拟高分子材料的结构和性能， 探究高分子材料的物理和化学性质， 优化材料的设计和制备过程，为新材 料的研发提供理论指导。
分子动力学方法需要较高的计算资源和 精度，但可以获得较为准确的结果，因 此在计算化学、生物学、材料科学等领
域得到广泛应用。
介观模拟的原理
介观模拟是一种介于微观和宏观之间的模拟方 法，通过模拟一定数量的粒子的相互作用和演 化来研究介观尺度的结构和性质。
介观模拟方法通常采用格子波尔兹曼方法、粒 子流体动力学等方法，适用于模拟流体、表面 、界面等介观尺度的问题。
分子模拟基于量子力学、经典力 学、蒙特卡洛等理论，通过建立 数学模型来描述分子间的相互作
用和运动。
分子模拟可以用于药物研发、材 料科学、环境科学等领域，为实 验研究和工业应用提供重要支持
。
分子模拟的重要性
01
02
03
预测分子性质
通过模拟，可以预测分子 的性质，如稳定性、溶解 度、光谱等，为实验设计 和优化提供指导。
分子模拟课件
ME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout1.An Overview of Molecular SimulationWei Caic All rights reservedSeptember26,2005Contents1Examples of molecular simulations1 2What consists of a molecular dynamics?2 3Newton’s equation of motion2 4A simple numerical integrator:Verlet algorithm5 1Examples of molecular simulationsPictures and movies of modern large scale molecular simulations can be downloaded from the course web site /∼me346(click“Course Materials”).Exam-ples include plastic deformation in metals containing a crack,phase transformation induced by shock waves,nano-indentation of metalfilm,laser ablation on metal surface,etc.In each of these examples,molecular simulations not only produce beautiful pictures,but they also offer new understanding on the behavior of materials at the nano-scale.Many times such understanding provides new insights on how the materials would behave at the micro-scale and nano-scale.With unique capabilities to offer,molecular simulation has become a widely used tool for scientific investigation that complements our traditional way of doing science: experiment and theory.The purpose of this course is to provide thefirst introduction to molecular level simulations—how do they work,what questions can they answer,and how much can we trust their results—so that the simulation tools are no longer black boxes to us.Students will develop a set of simulation codes in Matlab by themselves.(A simple Matlab tutorial written by Sigmon[5]may help you refresh your memory.)We will also give tutorial and exercise problems based on MD++,a molecular simulation package used in my research group.MD++will be used for illustration in class as well as tools to solve homework problems.2What consists of a molecular dynamics?In a molecular simulation,we view materials as a collection of discrete atoms.The atoms interact by exerting forces on each other and they follow the Newton’s equation of motion. Based on the interaction model,a simulation compute the atoms’trajectories numerically. Thus a molecular simulation necessarily contains the following ingredients.•The model that describes the interaction between atoms.This is usually called the interatomic potential:V({r i}),where{r i}represent the position of all atoms.•Numerical integrator that follows the atoms equation of motion.This is the heart of the ually,we also need auxiliary algorithms to set up the initial and bound-ary conditions and to monitor and control the systems state(such as temperature) during the simulation.•Extract useful data from the raw atomic trajectory pute materials properties of interest.Visualization.We will start our discussion with a simple but complete example that illustrates a lot of issues in simulation,although it is not on the molecular scale.We will simulate the orbit of the Earth around the Sun.After we have understood how this simulation works(in a couple of lectures),we will then discuss how molecular simulations differ from this type of “macro”-scale simulations.3Newton’s equation of motionThroughout this course,we limit our discussions to classical dynamics(no quantum mechan-ics).The dynamics of classical objects follow the three laws of Newton.Let’s review the Newton’s laws.Figure1:Sir Isaac Newton(1643-1727England).First law:Every object in a state of uniform motion tends to remain in that state of motion unless an external force is applied to it.This is also called the Law of inertia.Second law:An object’s mass m,its acceleration a,and the applied force F are related byF=m a.The bold font here emphasizes that the acceleration and force are vectors.Here they point to the same direction.According to Newton,a force causes only a change in velocity(v)but is not needed to maintain the velocity as Aristotle claimed(erroneously).Third law:For every action there is an equal and opposite reaction.The second law gives us the equation of motion for classical particles.Consider a particle (be it an atom or the Earth!depending on which scale we look at the world),its position is described by a vector r=(x,y,z).The velocity is how fast r changes with time and is also a vector:v=d r/dt=(v x,v y,v z).In component form,v x=dx/dt,v y=dy/dt,v z=dz/dt. The acceleration vector is then time derivative of velocity,i.e.,a=d v/dt.Therefore the particle’s equation of motion can be written asd2r dt2=Fm(1)To complete the equation of motion,we need to know how does the force F in turn depends on position r.As an example,consider the Earth orbiting around the Sun.Assume the Sun isfixed at origin(0,0,0).The gravitational force on the Earth is,F=−GmM|r|2·ˆe r(2)whereˆe r=r/|r|is the unit vector along r direction.|r|=√x2+y2+z2is the magnitudeof r.m is the mass of the Earth(5.9736×1024kg),M is mass of the Sun(1.9891×1030kg), and G=6.67259×10−11N·m2/kg2is gravitational constant.If we introduce a potential field(the gravitationalfield of the Sun),V(r)=−GmM|r|(3)Then the gravitational force can be written as minus the spatial derivative of the potential,F=−dV(r)d r(4)In component form,this equation should be interpreted as,F x=−dV(x,y,z)dxF y=−dV(x,y,z)dyF z=−dV(x,y,z)dzNow we have the complete equation of motion for the Earth:d2r dt2=−1mdV(r)d r(5)in component form,this reads,d2x dt2=−GMx|r|3d2y dt2=−GMy|r|3d2z dt2=−GMz|r|3In fact,most of the system we will study has a interaction potential V(r),from which the equation of motion can be specified through Eq.(5).V(r)is the potential energy of the system.On the other hand,the kinetic energy is given by the velocity,E kin=12m|v|2(6)The total energy is the sum of potential and kinetic energy contributions,E tot=E kin+V(7) When we express the total energy as a function of particle position r and momentum p=m v, it is called the Hamiltonian of the system,H(r,p)=|p|22m+V(r)(8)The Hamiltonian(i.e.the total energy)is a conserved quantity as the particle moves.To see this,let us compute its time derivative,dE tot dt =m v·d vdt+dV(r)d r·d rdt=md rdt·−1mdV(r)d r+dV(r)d r·d rdt=0(9)Therefore the total energy is conserved if the particle follows the Newton’s equation of motion in a conservative forcefield(when force can be written in terms of spatial derivative of a potentialfield),while the kinetic energy and potential energy can interchange with each other.The Newton’s equation of motion can also be written in the Hamiltonian form,d r dt =∂H(r,p)∂p(10)d p dt =−∂H(r,p)∂r(11)Hence,dH dt =∂H(r,p)∂r·d rdt+∂H(r,p)∂p·d pdt=0(12)Figure2:Sir William Rowan Hamilton(1805-1865Ireland).4A simple numerical integrator:Verlet algorithmA dynamical simulation computes the particle position as a function of time(i.e.its trajec-tory)given its initial position and velocity.This is called the initial value problem(IVP). Because the Newton’s equation of motion is second order in r,the initial condition needs to specify both particle position and velocity.To solve the equation of motion on a computer,thefirst thing we need is to discretize the time.In other words,we will solve for r(t)on a series of time instances t ually the time axis is discretized uniformly:t i=i·∆t,where∆t is called the time step.Again,the task of the simulation algorithm is tofind r(t i)for i=1,2,3,···given r(0)and v(0).The Verlet algorithm begins by approximating d2r(t)/dt2,d2r(t) dt2≈r(t+∆t)−2r(t)+r(t−∆t)∆t2(13)Thusr(t+∆t)−2r(t)+r(t−∆t)∆t2=−1mdV(r(t))dr(14)r(t+∆t)=2r(t)−r(t−∆t)−∆t·1m·dV(r(t))dr(15)Therefore,we can solve r(t+∆t)as along as we know r(t)and r(t−∆t).In other words,if we know r(0)and r(−∆t),we can solve for all r(t i)for i=1,3,4,···.Question:Notice that to start the Verlet algorithm we need r(0)and r(−∆t).However, the initial condition of a simulation is usually given in r(0)and v(0).How do we start the simulation when this initial condition is specified?References1.Alan and Tildesley,Computer Simulation of Liquids,(Oxford University Press,1987)pp.71-80.2.Frenkel and Smit,Understanding Molecular Simulation:From Algorithms to Applica-tions,2nd ed.(Academic Press,2002).3.Bulatov and Cai,Computer Simulation of Dislocations,(Oxford University Press,2006)pp.49-51.ndau and Lifshitz,Mechanics,(Butterworth-Heinemann,1976).5.Kermit Sigmon,Matlab Primer,Third Edition,/Matlab/matlab-primer.pdfME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout2.Numerical Integrators IWei Caic All rights reservedJanuary12,2007Contents1Which integrator to use?1 2Several versions of Verlet algorithm2 3Order of accuracy3 4Other integrators4 1Which integrator to use?In the previous lecture we looked at a very simple simulation.The Newton’s equation of motion is integrated numerically by the original Verlet algorithm.There are many algorithms that can be used to integrate an ordinary differential equation(ODE).Each has a its own advantages and disadvantages.In this lecture,we will look into more details at the numerical integrators.Our goal is to be able to make sensible choices when we need to select a numerical integrator for a specific simulation.What questions should we ask when choosing an integrator?Obviously,we would like something that is both fast and accurate.But usually there are more issues we should care about.Here are a list of the most important ones(Allen and Tildesley1987):•Order of accuracy•Range of stability•Reversibility•Speed(computation time per step)•Memory requirement(as low as possible)•Satisfying conservation law•Simple form and easy to implementItem1and2(especially2)are important to allow the use of long time steps,which makes the simulation more efficient.2Several versions of Verlet algorithmThe Verlet algorithm used in the previous lecture can be written as,r(t+∆t)=2r(t)−r(t−∆t)+a(t)·∆t2(1)wherea(t)=−1mdV(r(t))d r(t)(2)While this is enough to compute the entire trajectory of the particle,it does not specify the velocity of the particle.Many times,we would like to know the particle velocity as a function of time as well.For example,the velocity is needed to compute kinetic energy and hence the total energy of the system.One way to compute velocity is to use the followingapproximation,v(t)≡d r(t)dt≈r(t+∆t)−r(t−∆t)2∆t(3)Question:if the initial condition is given in terms of r(0)and v(0),how do we use the Verlet algorithm?Answer:The Verlet algorithm can be used to compute r(n∆t)for n=1,2,···provided r(−∆t)and r(0)are known.The task is to solve for r(−∆t)given r(0)and v(0).This is achieved from the following two equations,r(∆t)−r(−∆t)2∆t=v(0)(4)r(∆t)−2r(0)+r(−∆t)∆t2=−1mdV(r(t))d r(t)t=0(5)Exercise:express r(−∆t)and r(∆t)in terms of r(0)and v(0).Therefore,r(−∆t)and r(∆t)are solved together.During the simulation,we only know v(t)after we have computed(t+∆t).These are the(slight)inconveniences of the(original) Verlet algorithm.As we will see later,in this algorithm the accuracy of velocity calculation is not very high(not as high as particle positions).Several modifications of the Verlet algorithms were introduced to give a better treatment of velocities.The leap-frog algorithm works as follows,v(t+∆t/2)=v(t−∆t/2)+a(t)·∆t(6)r(t+∆t)=r(t)+v(t+∆t/2)·∆t(7) Notice that the velocity and position are not stored at the same time slices.Their values are updated in alternation,hence the name“leap-frog”.For a nice illustration see(Allen and Tildesley1987).This algorithm also has the draw back that v and r are not known simultaneously.To get velocity at time t,we can use the following approximation,v(t)≈v(t−∆t/2)+v(t+∆t/2)2(8)The velocity Verlet algorithm solves this problem in a better way.The algorithm reads,r(t+∆t)=r(t)+v(t)∆t+12a(t)∆t2(9)v(t+∆t/2)=v(t)+12a(t)∆t(10)a(t+∆t)=−1mdV(r(t+∆t))d r(t+∆t)(11)v(t+∆t)=v(t+∆t/2)+12a(t+∆t)∆t(12)Exercise:show that both leap-frog and velocity Verlet algorithms gives identical results as the original Verlet algorithm(except forfloating point truncation error).3Order of accuracyHow accurate are the Verlet-family algorithms?Let’s take a look at the original Verletfirst. Consider the trajectory of the particle as a continuous function of time,i.e.r(t).Let us Taylor expand this function around a given time t.r(t+∆t)=r(t)+d r(t)dt∆t+12d2r(t)dt2∆t2+13!d3r(t)dt3∆t3+O(∆t4),(13)where O(∆t4)means that error of this expression is on the order of∆t4.In other words, (when∆t is sufficiently small)if we reduce∆t by10times,the error term would reduce by 104times.Replace∆t by−∆t in the above equation,we have,r(t−∆t)=r(t)−d r(t)dt∆t+12d2r(t)dt2∆t2−13!d3r(t)dt3∆t3+O(∆t4),(14)Summing the two equations together,we have,r(t+∆t)+r(t−∆t)=2r(t)+d2r(t)dt2∆t2+O(∆t4)≡2r(t)+a(t)∆t2+O(∆t4)(15)r(t+∆t)=2r(t)−r(t−∆t)+a(t)∆t2+O(∆t4)(16) Therefore,the local accuracy for the Verlet-family algorithms is to the order of O(∆t4). This is quite accurate for such a simple algorithm!The trick here is that,by adding t+∆t and t−∆t terms,the third order term cancels out.Notice that the accuracy analysis here deals with the local error.That is,assuming r(t)is exact,how much error we make when computing r(t+∆t).In practice,the total error of computed r at time t is the accumulated (and possibly amplified)effect of all time steps before time t.Because for a smaller timestep ∆t more steps is required to reach a pre-specified time t,the global order of accuracy is always local order of accuracy minus one.For example,the global order of accuracy for Verlet algorithms is3.Since all Verlet-family algorithms generate the same trajectory,they are all O(∆t4)in terms of particle positions.But they differ in velocity calculations.How accurate is the velocity calculation in the original Verlet algorithm?Let us subtract Eq.(13)and(14),r(t+∆t)−r(t−∆t)=2v(t)∆t+13d3r(t)dt3∆t3+O(∆t4),(17)v(t)=12∆tr(t+∆t)−r(t−∆t)−16d3r(t)dt2∆t3+O(∆t3)(18)v(t)=r(t+∆t)−r(t−∆t)2∆t+O(∆t2)(19)Thus the velocity is only accurate to the order of O(∆t2)!This is much worse than O(∆t4) of the positions.A reason for the lose of accuracy is that the velocities are computed by taking the difference between two large numbers,i.e.r(t+∆t)and r(t−∆t),and then divide by a small number.1Intuitively,higher order algorithms allows the use of bigger time steps to achieve the same level of accuracy.However,the more important determining factor for the choice of time steps is usually not the order of accuracy,but instead is the algorithm’s numerical stability (see next section).4Other integratorsBefore closing this lecture,let us briefly look at two other well known integrators.Gear’s predictor-corrector method is a high-order integrator.It keeps track of several time deriva-tives of the particle position.Definer0(t)≡r(t)(20) 1Strictly speaking,Eq.(17)is valid when we have exact values for r(t+∆t),r(t)and r(t−∆t).But in doing the order-of-accuracy analysis,we can only assume that we have exact values for r(t)and r(t−∆t) and we know that we are making an error in r(t+∆t).Fortunately,the error we make in r(t+∆t)is also O(∆t4).Thus Eq.(17)still holds.r1(t)≡∆t d r(t)dt(21)r2(t)≡∆t22!d2r(t)dt2(22)r3(t)≡∆t33!d3r(t)dt3(23)···Sincer(t+∆t)=r(t)+∆t d r(t)dt+∆t22!d2r(t)dt2+∆t33!d3r(t)dt3+O(∆t4)=r(t)+r1(t)+r2(t)+r3(t)+O(∆t4)r1(t+∆t)=r1(t)+2r2(t)+3r3(t)+O(∆t4)r2(t+∆t)=r2(t)+3r3(t)+O(∆t4)r3(t+∆t)=r3(t)+O(∆t4)we can predict the values of r and r i at time t+∆t simply based on their values at time t,r P(t+∆t)=r(t)+r1(t)+r2(t)+r3(t)r P1(t+∆t)=r1(t)+2r2(t)+3r3(t)r P2(t+∆t)=r2(t)+3r3(t)r P3(t+∆t)=r3(t)Notice that we have not even calculated the force at all!If we compute the force at time t+∆t based on the predicted position r P(t+∆t),we will obtain the second derivative of r at time t+∆t—let it be r C(t+∆t).Mostly likely,r C(t+∆t)will be different from the predicted value r P(t+∆t).Thefinal results are obtained by adding a correction term to each predictor,r n(t+∆t)=r Pn (t+∆t)+C n[r C2(t+∆t)−r P2(t+∆t)](24)The Gear’s6th order scheme(uses r0through r5)takes the following C n values,C0=3 16C1=251 360C2=1C3=11 18C4=1 6C5=1 60The global accuracy is of the order O(∆t)6,hence the method is called the Gear’s6th order algorithm.Exercise:show that the local truncation error of the above scheme is of the order of O (∆t 7).In a way,the Runge-Kutta method is similar to predictor-corrector methods in that the final answer is produced by improving upon previous predictions.However,the Runge-Kutta method requires multiple evaluation of forces.The Runge-Kutta method is most conveniently represented when the equation of motion is written as a first order differential equation,˙η=g (η)≡ω∂H (η)∂η(25)whereη≡(r ,p )(26)andω= 0I−I 0 (27)and I is the 3×3identity matrix and H is the Hamiltonian.The fourth-order Runge-Kutta method goes as (from ),k 1=∆t ·g (η(t ))k 2=∆t ·g η(t )+12k 1 k 3=∆t ·g η(t )+12k 2 k 4=∆t ·g (η(t )+k 3)η(t +∆t )=η(t )+16k 1+13k 2+13k 3+16k 4(28)Exercise:show that the local truncation error of the above scheme is of the order of O (∆t 5)and the global accuracy is of the order of O (∆t 4).Suggested reading1.(Alan and Tildesley 1987)pp.71-84.2.J.C.Butcher,The numerical analysis of ordinary differential equations:Runge-Kutta and general linear methods ,(Wiley,1987),pp.105-151.ME346–Introduction to Molecular Simulations–Stanford University–Winter2007 Handout3.Symplectic ConditionWei Caic All rights reservedJanuary17,2007Contents1Conservation laws of Hamiltonian systems1 2Symplectic methods4 3Reversibility and chaos5 4Energy conservation7 1Conservation laws of Hamiltonian systemsThe numerical integrators for simulating the Newton’s equation of motion are solvers of ordinary differential equations(ODE).However,the Newton’s equation of motion that can be derived from a Hamiltonian is a special type of ODE(not every ODE can be derived from a Hamiltonian).For example,a Hamiltonian system conserves it total energy,while an arbitrary ODE may not have any conserved quantity.Therefore,special attention is needed to select an ODE solver for Molecular Dynamics simulations.We need to pick those ODE solvers that obeys the conservation laws of a Hamiltonian system.For example,if an ODE solver predicts that the Earth will drop into the Sun after100cycles around the Sun,it is certainly the ODE solver to blame and we should not immediately start worrying about the fate of our planet based on this prediction.A conservation law is a(time invariant)symmetry.A Hamiltonian system possess more symmetries than just energy conservation.Ideally the numerical integration scheme we pick should satisfy all the symmetries of the true dynamics of the system.We will look at these symmetries in this and following sections.When discussing a Hamiltonian system,it is customary to use letter q to represent the particle position(instead of r).Recall that the equation of motion for a system with Hamiltonian H(q,p)is,˙q=∂H(q,p)∂p(1)˙p=−∂H(q,p)∂q,(2)where˙q≡d q/dt,˙p≡d p/dt.If the system has N particles,then q is a3N-dimensional vector specifying its position and p is a3N-dimensional vector specifying its momentum. (In the above Earth orbiting example,N=1.)Notice that in terms of q and p,the equation of motion is a set of(coupled)first-order differential equations(whereas F=m a is second order).The equation can be reduced to an even simpler form if we define a6N-dimensional vectorη≡(q,p)(3) Then the equation of motion becomes,˙η=ω∂H(η)∂η(4)whereω=0I−I0(5)and I is the3×3identity matrix.The6N-dimensional space whichη≡(q,p)belongs to is called the phase space.The evolution of system in time can be regarded as the motion of a pointηin the6N-dimensional phase space following thefirst-order differential equation(4).(a)(b)Figure1:(a)Motion of a point(q,p)in phase space is confined to a shell with constant energy.(b)Area of a phase space element(shaded)is conserved.Because the Newton’s equation of motion conserves total energy,the motion of a point in phase space must be confined to a subspace(or hyperplane)with constant energy,i.e.a constant energy shell,as illustrated in Fig.1(a).Think about all the points in the constant energy shell.As time evolves,the entire volume that these points occupy remains the same —the points simply rearrange within the constant energy shell.Let us now consider the evolution of a small element in phase space over time,as illus-trated in Fig.1(b).Think of this as the ensemble of trajectories of many Molecular Dynamics simulations with very similar initial conditions.We will show that the area enclosed by this element remains constant,even though the element inevitably experiences translation anddistortion.1Let the element at time t be a hypercube around pointη,whose area is,|dη|=|d q|·|d p|=dq1···dq3N dp1···dp3N(6) For a pointηat time t,letξbe its new position at time t+δt.In the limit ofδt→0,wehave,ξ=η+˙ηδt=η+ω∂H(η)∂ηδt(7)Because the above equation is true only in the limit ofδt→0,we will ignore any higher order terms(e.g.δt2)in the following discussions.Let M be the Jacobian matrix of the transformation fromηtoξ,i.e.,M≡∂ξ∂η=1+ω∂2H(η)∂η∂ηδt(8)For clarity,we write M explicitly in terms of q and p,M=1+0I−I0·∂2H∂q∂q∂2H∂q∂p∂2H∂p∂q∂2H∂p∂p·δt(9)=1+∂2Hδt∂2Hδt−∂2H∂q∂qδt1−∂2H∂q∂pδt(10)The area of the new element is related to the the determinant of the Jacobian matrix,|dξ|=|dη|·|det M|=|dη|·(1+O(δt2))(11) Because thefirst order term ofδt in det M vanishes,we can show that the area of the element remains constant for an arbitrary period of time(i.e.forever).2This is an important property of a Hamiltonian system.Because the area of any element in phase space always remains constant,the evolution of phase space points is analogous to that in an incompressiblefluid.Exercise:Write down the explicit expression for the Jacobian matrix M for problem of the Earth orbiting around the Sun.The Hamiltonian dynamics has even more symmetries.Notice that the transpose of M is,M T=1−∂2H∂η∂ηωδt(12)1It then follows that the area enclosed by an arbitrarily shaped element in the phase space also remains constant.2To show that area remains constant after afinite period∆t,we only need to divide the time interval into N subintervals,each withδt=∆t/N.The area change per subinterval is of the order of1/N2.The accumulated area change over time period∆t is then of the order of1/N,which vanishes as N→∞.we haveM ωM T = 1+ω∂2H ∂η∂ηδt ω 1−∂2H ∂η∂ηωδt =ω+ω∂2H ∂η∂ηωδt −ω∂2H ∂η∂ηωδt +O (δt 2)(13)=ω+O (δt 2)(14)(15)Therefore the Jacobian matrix M thus satisfies the so called symplectic condition (up to O (δt 2)),M ωM T =ω(16)Again,we can show that the symplectic condition is satisfied for an arbitrary period of time.To be specific,let ξbe the point in phase space at time t ;at time 0this point was at η.ξcan then be regarded as a function of t and η(initial condition),i.e.ξ=ξ(η,t ).Define M as the Jacobian matrix between ξand η,we can show that M satisfies the symplectic condition Eq.(16).3Obviously,the Hamiltonian dynamics is also time reversible.This means that if ηevolves to ξafter time t .Another phase space point η that has the same q as ηbut with reversed p (momentum)will exactly come to point ξafter time t .To summarize,the dynamics of a system that has a Hamiltonian has the following symmetry properties:•Conserves total energy•Conserves phase space area (incompressible flow in phase space)•Satisfies symplectic condition Eq.(16)•Is time reversibleIdeally,the numerical integrator we choose to simulate Hamiltonian system should satisfy all of these symmetries.In the following,let us take a look at how specific integrators are doing in this regard.In fact,if an algorithm satisfies the symplectic condition,it automatically conserves phase space area 4,is reversible,and should have good long-term energy conserva-tion properties.This is why symplectic methods are good choices for simulating Hamiltonian systems.2Symplectic methodsWhat do we mean if we ask whether a numerical integrator satisfies the symplectic condition or not?For a numerical integrator,the system moves forward in time by discrete jumps.Let3Again,we devide the time period t into N subintervals,each with δt =∆t/N .Let M i be the Jacobian matrix from time (i −1)δt to iδt .We have M = N i =1M i .We can show that M T ωM =ωin the limit of N →∞.4Since M T ωM =ω,det ω=det(M T ωM )=(det M )2det ω.Thus det M =±1,which means area conservation.η(i)=(q(i),p(i))be the place of the system in the phase space at time step i,i.e.t i=i∆t. At the next time step,the algorithm brings the system to pointη(i+1)=(q(i+1),p(i+1)). Let M be the Jacobian matrix that connectsη(i)withη(i+1).The method is symplectic if MωM T=ω.The good news is that the Verlet-family algorithms are symplectic.To show this is the case,let us consider the leap-frog algorithm.Let q(i)correspond to q(i·∆t)and p(i) correspond to p((i−12)·∆t).Then the leap-frog algorithm is,p(i+1)=p(i)−∂V∂q(i)∆t(17)q(i+1)=q(i)+p(i+1)m∆t=q(i)+1mp(i)−∂V∂q(i)∆t∆t(18)Thus,M=∂q(i+1)∂q(i)∂q(i+1)∂p(i)∂p(i+1)(i)∂q(i+1)(i)=1−1m∂2V∂q(i)∂q(i)∆t2∆tm−∂2V∂q(i)∂q(i)∆t1(19)Thus M is symplectic.Exercise:show that for M defined in Eq.(19),MωM T=ω(symplectic)and det M=1 (area preserving).Exercise:derive the Jacobian matrix M for velocity Verlet algorithm and show that it is symplectic and area preserving.3Reversibility and chaosBecause the Verlet-family algorithms(original Verlet,leap-frog and velocity Verlet)are sym-plectic,they are also time reversible.Nonetheless,it is instructive to show the time reversibil-ity explicitly.If we run the leap-frog algorithm in reverse time,then the iteration can be written as,˜q(i+1)=˜q(i)+˜p(i)m∆t(20)˜p(i+1)=˜p(i)−∂V∂q(i+1)∆t(21)(22)where i is the iteration number,˜q(i)=q(−i∆t),˜p(i)=p((−i−12)∆t).It is not difficult toshow that if˜q(i)=q(i+1)and˜p(i)=−p(i+1),then˜q(i+1)=q(i)and˜p(i+1)=−p(i).Exercise:show that the leap-frog algorithm is reversible.Exercise:show that the velocity Verlet algorithm is reversible.Exercise:show that the Forward Euler method:q(i+1)=q(i)+∆t·p(i)/m,p(i+1)=p(i)−∆t·∂V/∂q(i),is not reversible and not symplectic.Since every integration step using the Verlet algorithm is time reversible,then in principle an entire trajectory of Molecular Dynamics simulation using the Verlet algorithm should be time reversible.To be specific,let the initial condition of the simulation be q(0)and p(0). After N steps of integration,the system evolves to point q(N)and p(N).Time reversibility means that,if we run the reverse Verlet iteration starting from q(N)and−p(N),after N steps we should get exactly to point q(0)and−p(0).In practice,however,we never have perfect reversibility.This is due to the combined effect of numerical round offerror and the chaotic nature of trajectories of many Hamiltonian systems.The chaotic nature of many Hamiltonian system can be illustrated by considering the two trajectories with very close initial conditions.Let the distance between the two points in phase space at time0be d(0).At sufficiently large t,the two phase space points will diverge exponentially,i.e.d(t)∼d(0)eλt,whereλis called the Lyapunov exponent and this behavior is called Lyapunov instability.Given that we can only represent a real number on a computer withfinite precision,we are in effect introducing a small(random) perturbation at every step of the integration.Therefore,sooner or later,the numerical trajectory will deviate from the analytical trajectory(if the computer had infinite precision) significantly5.While the analytical trajectory is reversible,the one generated by a computer is not.This is why there exist a maximum number of iterations N,beyond which the original state cannot be recovered by running the simulation backwards.This by itself does not present too big a problem,since we seldom have the need to recover the initial condition of the simulation in this way.However,this behavior does mean that we will not be able to follow the“true”trajectory of the original system forever on a computer.Sooner or later, the“true”trajectory and the simulated one will diverge significantly.This is certainly a problem if we would like to rely entirely on simulations to tell us where is our satellite or space explorer.However,this usually does not present a problem in Molecular Simulations, in which we are usually not interested in the exact trajectories of every atom.We will return to this point in future lectures when we discuss Molecular Dynamics simulations involving many atoms.5Notice that the analytical trajectory here is not the same as the“true”trajectory of the original system. The analytical trajectory is obtained by advancing the system in discrete jumps along the time axis.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第十章 计算化学软件
计算机商业计算软件之介绍
几种常用的计算软件介绍算软件之介绍
涉及领域广泛
所有与化学有关 的科研领域都有 许多利用计算机 成功模拟的研究。
计算机商业计算软件之介绍
优点
计算机商业计算软件之介绍
目前流行最广的计算软件为美国的Accelrys公司所 推出的各项产品。
几种常用的计算软件介绍——Materials Studio
几种常用的计算软件介绍——Materials Studio
几种常用的计算软件介绍——VASP
几种常用的计算软件介绍——VASP
VASP：(Vienna Ab-initio Simulation Package)
VASP是使用赝势和平面波基组，进行从头量子力 学分子动力学计算的软件包，是目前材料模拟和计算 物质科学研究中最流行的商用软件之一。
几种常用的计算软件介绍——GAUSSIAN
Gaussian：
是一个功能强大的 量子化学综合软件包。
几种常用的计算软件介绍——GAUSSIAN
计算可以对体系的 基态或激发态执行。可 以预测周期体系的能量， 结构和分子轨道。因此， Gaussian可以作为功能 强大的工具，用于研究 许多化学领域的课题， 例如取代基的影响，化 学反应机理，势能曲面 和激发能等等。
总能量和部分能量； 原子力； 应力张量； 电偶极矩； 原子，轨道和键分析 ； 电子密度； 常温分子动力学； 可变晶胞动力学 ； 其它的线性标度方法； 增强的MD历史框架； BZ区的k-取样； 自旋极化计算(共线或者非共线)； 几何松弛，固定或者改变晶胞； 态的局域和轨道投影密度； 能带结构； HF和混和泛函； 通过过滤或移到原子格点的方法平滑“蛋箱效应”； 多格点方法计算溶剂中分子的Poisson-Boltzman方程；
几种常用的计算软件介绍——SIESTA
SIESTA：
(Spanish Initiative for Electronic Simulations with Thousands of Atoms) 用于分子和固体的电子结构计算和分子动力学 模拟。
几种常用的计算软件介绍——SIESTA
SIESTA主要有以下功能:
应用领域涵盖
计算机商业计算软件之介绍
另外几种常用的计算软件：
VASP
GAUSSIAN
SIESTA
计算机商业计算软件之介绍
几种常用的计算软件介绍
计算机计算的回顾与发展
几种常用的计算软件介绍——Materials Studio
Materials Studio：是专门为材料科学领域研究者
开发的一款可运行在PC上的模拟软件。该软件能使任 何研究者达到与世界一流研究部门相一致的材料模拟的 能力。模拟的内容包括了催化剂、聚合物、固体及表面、 晶体与衍射、化学反应等材料和化学研究领域的主要课 题。
几种常用的计算软件介绍——Materials Studio
几种常用的计算软件介绍——Materials Studio
此模块的主要功能为高品质的图形显示，可显 示各种分子，晶体与高分子。亦可以提供分子模拟 的功能。可由所建立的分子图形，结合该软件内含 的化学数据库，找出最低能量构形。
下图为利用此模块所得的一些分子与晶体的结 构的图形示例：
该模块主要功能为执行相平衡的模拟计算，有助 于化工与石化工业的制成设计。
几种常用的计算软件介绍——Materials Studio
此模块主要应用于衍射光谱，晶体结构的模拟计 算。Reflex模块由所建立的晶体结构计算可能的X-光衍 射，中子衍射与电子衍射的粉状样品图谱。
该模块为Reflex模块的推广，主要功能为由高分 辨率实验的X-光衍射图谱推出晶体的结构。
几种常用的计算软件介绍——Materials Studio
该模块提供特殊的力场，此力场适用于凝态材质 的模拟计算。此力场由精确的量子计算及各种实验的 数据所得出，其目的在于准确的计算出分子及系统于 不同的温度及压力条件下的结构，振动光谱，热力学 特性。此力场目前适用于有机分子，无机氧化物分子， 高分子及金属离子等系统。
计算机商业计算软件之介绍 几种常用的计算软件介绍
计算机计算的回顾与发展
计算机计算的回顾与发展
利用计算机计算进行研究约从上世纪六十年代开始
计算机飞速发展
计算机已成为高质量研究所不可少的重要工具
计算机计算的回顾与发展
计算机计算
优 点 一般实验
因此，可预见计算机模拟将可取代许多传统的实 验，主导日后科研的方向。
几种常用的计算软件介绍——Materials Studio
此模块的主要功能为建立，模拟计算与分析非 晶形聚合物系统。该模块可用以预测聚合物的性质。
该模块为Material Studio软件的核心计算模块。 此模块以经典力学的方法结合适当的全原子力场以 模拟计算各种化学与材料体系。可用以研究催化， 分离，结晶及高分子体系的特性。利用此模块可预 测固体，液体及气体的重要特性。