当前位置:首页 » 网页前端 » 分子动力学能量脚本
扩展阅读
webinf下怎么引入js 2023-08-31 21:54:13
堡垒机怎么打开web 2023-08-31 21:54:11

分子动力学能量脚本

发布时间: 2023-04-01 08:11:41

㈠ 纵轴为RMSD值,横轴为时间的图是利用什么软件做出来的如下图(与分子动力学模拟有关)

按照教程来,先跑完动力学,然后用脚本获得RMSD值,然后根据你模拟的时间计算出每一帧的对应的时间步长,然后将数据拿到origin里面作图

㈡ 分子反应动力学的举例

例如,为了使吸能反应 I+HCl─→HI+Cl
能够发嫌闷生,增加 HCl的振动能比增加其平动能更为有效。它的逆反应 Cl+HI─→HCl+I
是一个放能反应,分子反应动力学能够提供产物分子HCl振动态“布居反转”的信息,从而为寻找化学激光工作物质提供了依据。它还能提供反应体系“碰撞对”真实碰撞过程的信息──“碰撞对”是直接反应还是经过一个络合物的反应。
理论计算方法 20世纪30年代,以美国物理化学家H.艾林为代表的学派,用海特勒-伦敦计算H2的方法建立了H+H2反应体系的第一个势能面,借助统计力学方法计算了在该势能面上的热平衡反应速率常数,称为绝对芹亮弯反应速率理论或过渡态理论。
分子反应动力学的理论计算方法分为三部分:①化学反应体系势能面的量子化学计算;②反应截面(或几率)的计算;③由反应截面计算反应速率常数。因此,也可以说分子反应动力学是研究反应体系在热能面上运动过程的学科。在确定的势能面上求解核的运动方程,既可以用经典力学方法,也可以用量子力学方法。
理论 严格的理论是量子力学散射理论。分子反应过程的全部信息包含在波函数中,在给定能量下,求解满足一定渐近条件的薛定谔方程得到波函数,借助入射波和出射波的几率流密度守恒的关系,键核就可以得到反应截面(或几率)。
以A+BC─→AB+C双分子共线交换反应为例 (共线反应是指反应体系的三个原子沿直线相互接近的反应),该反应体系的坐标系见图1。
在非相对论近似下,反应体系的哈密顿算符H写作(图2):
式中μA,BC和μBC分别为A和BC,B和C之间相对运动的约化质量;mA、mB、mC分别为原子A、B、C的质量;h为普朗克常数;Vα和Vγ为有效势函数。
核运动的薛定谔方程为: Hψ=Eψ (3)
渐近条件为:(图3) 式中α为反应体系的初始排布,即A+BC;nα或n为BC的内量子数,nα为始态,n为反射态;γ表示终态排布,即C+AB;n为AB分子的内量子数,每一种排布和分子的一组内量子数(如α,nα)称为反应体系的一个通道;kα或kγ为原子与双原子分子相对运动的波数(图4);
为双原子分子的内态波函数(图5);
称为散射幅。能量守恒条件要求(图6): (5)
式中啚=h/2π;E为能量。由入射波和出射波几率流密度守恒的条件,就可以得到由通道(α,nα)到通道(λ,nλ)的反应几率为(图7): 式中v为(λ,nλ)通道中反应体系的相对运动速度。
H+H2(n)─→H2+H共线交换反应几率的数值计算结果见图2。
对于实际的三维化学反应,用上面的方法可以得到反应截面随碰撞能变化的关系。用量子散射理论求反应截面(或几率)的关键是求散射幅,一般是在自然反应坐标中用数值求解耦合微分方程。这是一项十分复杂的计算工作。
当反应体系的质量较大,德布罗意波长很短时,用经典轨迹法或者用准经典轨迹法,即对反应物初态分布和产物终态分布作量子校正的经典轨迹法研究反应体系沿势能面的运动,往往也能得到比较满意的定性或半定量的结果。

㈢ 分子动力学简介

分子动力学(MD)是一种用于分析原子和分子物理运动的计算机模拟方法。原子和分子被允许在一段固定的时间内相互作用,从而可以看到系统的动态“进化”。
如下图所示,当一个原子沉积在一个晶体上时,由于原子之间的吸引力,它并没有反弹,而是保持附着状态。随着时间的迁移,我们可以看到从顶部接近的原子的动能在其他原子之间进行了重新分配。

一般来说,上述方程会转化为谈厅牛顿方程进行计算,原子和分子的轨迹是通过数值求解粒子相互作用系统的牛顿运动方程来确定的。其中粒子之间的力和它们的势能通常是通过原子间势或分子力学力场来计算的。

  键的伸缩能:成键原子间相互作用力
  键角弯曲能:成键角度的影响
  键的扭曲能:化学键旋转御侍裤时的能量变化
  非键作用:分子间相互作用力,
  非键作用:静电相互作用

分子力学力场方法的基本原理是最小作用量原理,在分子内部化学键都有“自然”的键长值和键角值。分子要调整它的几何形状(构象),以使其键长值和键角值尽可能接近自然值,同时也使非键作用处于最小的状态。
举个例子,原子之间靠近时相互排斥,为了减少这种作用,原子就趋于相互离开;但是这将使键长伸长或键角发生弯曲,又引起了相应的能量升高。最后的构型将是这两种力折衷的结果,并且是能量最低的构型”。

MD是一种基于牛顿力学确定论的热力学计算方法,与蒙特卡洛法相比在宏观性质计算上具有更高的准确度和有效性,可以广泛应用于物理,化学,生物,材料,医学等各个领域。
但是MD没有考虑电子转移效应,不能准确模拟化学成键;需要大量的参数;且获得的能量估计不是非常准确。因此,一种新的方法QM/镇简MM出现,它结合了QM(量子力学)的精度和MM(分子力学)的速度的优点。一般来说,它将系统的一小部分(如酶的活性位点)进行量子力学处理,其余部分进行经典处理。

㈣ 分子动力学的基本理论有哪些

分子动力学是一门结合物理,数学和化学的综合技术。
确定起始构型
进行分子动力学模拟的第一步是确定起始构型, 一个能量较低的起始构型是进行分子模拟 的基础 ,一般分子的起始构型主要来自实验数据或量子化学计算。 在确定起始构型之后要赋予构成分子的各个原子速度,这一速度是根据波尔兹曼分布随机生成的,由于速度的分布符合波尔兹曼统计,因此在这个阶段,体系的温度是恒定的。另外,在随机生成各个原子的运动速衫罩塌度之后须进行闷码调整,使得体系总体在各个方向上的动量之和为零,即保证体系没有平动位移。
进入平衡相
由上 一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。
进入生产相
进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个粒子的运动轨迹进行 计算 ,在这个过程中,体系总能量不变,但分子内部势能和动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正是在这个过程中抽取的。
+ 计算结果 用抽样所得体或圆系的各个状态计算当时体系的势能,进而计算构型积分。

㈤ 分子动力学的基本步骤

用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。作用势与动力学计算
作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动 和 分子内部运动的轨迹也会不同,进而影响到抽样的结果和抽样结果的势能计算,在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是相对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大的困难,因此在一般应用中,通过第一性原理计算结果拟合势函数的L-J,morse等势模型的应用仍然非常广泛。 以下是做模拟的一般性步骤,具体的步骤和过程依赖于确定的系统或者是软件,但这不影响我们把它当成一个入门指南:
1)首先我们需要对我们所要模拟的系统做一个简单的评估, 三个问题是我们必须要明确的:
做什么(what to do)为什么做(why to do)怎么做(how to do)
2)选择合适的模拟工具,大前提是它能够实现你所感兴趣的目标,这需要你非常谨慎的查阅文献,看看别人用这个工具都做了些什么,有没有和你相关的,千万不要做到一半才发现原来这个工具根本就不能实现你所感兴趣的idea,切记!
考虑1:软件的选择,这通常和软件主流使用的力场有关,而软件本身就具体一定的偏向性,比如说,做蛋白体系,Gromacs,Amber,Namd均可;做DNA, RNA体系,首选肯定是Amber;做界面体系,Dl_POLY比较强大,另外做材料体系,Lammps会是一个不错的选择
考虑2:力场的选择。力场是来描述体系中最小单元间的相互作用的,是用量化等方法计算拟合后生成的经验式,有人会嫌它粗糙,但是它确确实实给我们模拟大系统提供了可能,只能说关注的切入点不同罢了。常见的有三类力场:全原子力场,联合力场,粗粒化力场;当然还有所谓第一代,第二代,第三代力场的说法,这里就不一一列举了。
再次提醒注意:必须选择适合于我们所关注体系和我们所感兴趣的性质及现象的力场。
3)通过实验数据或者是某些工具得到体系内的每一个分子的初始结构坐标文件,之后,我们需要按我们的想法把这些分子按照一定的规则或是随机的排列在一起,从而得到整个系统的初始结构,这也是我们模拟的输入文件。
4)结构输入文件得到了,我们还需要力场参数输入文件,也就是针对我们系统的力场文件,这通常由所选用的力场决定,比如键参数和非键参数等势能函数的输入参数。
5)体系的大小通常由你所选用的box大小决定,我们必须对可行性与合理性做出评估,从而确定体系的大小,这依赖于具体的体系,这里不细说了。6)由于初始构象可能会存在两个原子挨的太近的情况(称之为bad contact),所以需要在正式模拟开始的第一步进行体系能量最小化,比较常用的能量最小化有两种,最速下降法和共轭梯度法,最速下降法是快速移除体系内应力的好方法,但是接近能量极小点时收敛比较慢,而共轭梯度法在能量极小点附近收敛相对效率高一些,所有我们一般做能量最小化都是在最速下降法优化完之后再用共轭梯度法优化,这样做能有效的保证后续模拟的进行。
7)以平衡态模拟为例,你需要设置适当的模拟参数,并且保证这些参数设置和力场的产生相一致,举个简单的例子,gromos力场是用的范德华势双截断来定范德华参数的,若你也用gromos力场的话也应该用双截断来处理范德华相互作用。常见的模拟思路是,先在NVT下约束住你的溶质(剂)做限制性模拟,这是一个升温的过程,当温度达到你的设定后, 接着做NPT模拟,此过程将调整体系的压强进而使体系密度收敛。
经过一段时间的平衡模拟,在确定系统弛豫已经完全消除之后,就可以开始取数据了。如何判断体系达到平衡,这个问题是比较技术性的问题,简单的讲可以通过以下几种方式,一,看能量(势能,动能和总能)是否收敛;二,看系统的压强,密度等等是否收敛;三看系统的RMSD是否达到你能接受的范围,等等。
8)运行足够长时间的模拟以确定我们所感兴趣的现象或是性质能够被观测到,并且务必确保此现象出现的可重复性。
9)数据拿到手后,很容易通过一些可视化软件得到轨迹动画,但这并不能拿来发文章。真正的工作才刚刚开始——分析数据,你所感兴趣的现象或性质只是表面,隐含在它们之中的机理才是文章中的主题。

㈥ 分子动力学能量和结构优化能量的区别

体系的中原子类型和位置一旦确定,体系的总势能也确定,因此通过移动体系的原子的位置,体系就在体系的势能面上运动。那无论是分子动力学还是结构优化都是如此。在分子动力学模拟时,体系达到平衡就是体系在势能面最低点附近振动。结构优化也是找这个最低点,

㈦ 分子动力学名词解释

分子动力学是一门结合物理,数学和化学的综合技术。分子动力学是一套分子模拟山此方法,该方法主要是依靠牛顿力学来模拟分子体系的运动渣唯斗,以在由分子体系的如磨不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。
进行分子动力学模拟的第一步是确定起始构型,一个能量较低的起始构型是进行分子模拟的基础,一般分子的起始构型主要来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度,这一速度是根据波尔兹曼分布随机生成的,由于速度的分布符合波尔兹曼统计,因此在这个阶段,体系的温度是恒定的。

㈧ 什么是从头计算分子动力学模拟

入门阶段,首先你要知道你想做什么,最好是找个看起来不太难的文章照着把里面的模拟自己重复一遍。因为全原子模拟大都是用一些来进行的,因此你首先需要的是学会一些的使用,常用的生物分子模拟包括:Gromacs、Amber 和 NAMD 等等,材料有关的模拟还有 Lammps 等。学这些东西的时候首先主要是要知道模拟的基本流程以及实现的方法,包括怎样搭建模拟的体系、各种文件格式的转换、系综与盒子的选择、水及离子、能量极小化等等,等到模拟的轨迹出来怎样对数据进行处理,等到之后还可以学习里面的一些插件,例如一些加速采样的方法等等。

自己学一种语言的话,在初期,做 MD 比较重要的是脚本语言,包括 Shell 脚本或者其它你自己喜欢的脚本。因为最终你还是不太可能完全在自己的电脑上跑程序的,所以要有一个你自己用得比较熟的、能对大规模的数据进行处理的语言,我觉得 Python 是很适合的,而且里面的 Prody,Matplotlib 等等各种包都非常好用。

入门之后,如果希望自己通过一些量子化学的计算结果去调整和修改现有的力场,那么需要能看懂其他人的代码,这种时候很可能会需要能读懂 Fortran 的代码。如果自己喜欢做一些简化模型自己弄着玩,用 Python 之类的写起来是简单,但是效率太低,还是需要会一点点 C 或者 C++,当然语言只是一方面,更重要的是自己要结合实际的体系做一些最简单的优化。

相比起书籍来,还可以关注一些做模拟的学术们聚集的论坛和社区,例如:小木虫、分子模拟论坛、ResearchGate 等等。

参考书的话,其实有很多,不过还是要看你自己需要哪方面的内容:
分子模拟方面的经典书籍:Understanding molecular simulation: From algorithms to applications 和 Molecular Modelling - Principles and Applications ,两本书的侧重点有些不同。

中文书籍:《分子模拟的理论与实践》《计算化学——从理论化学到分子模拟》中的部分章节;

偏统计和计算物理方面:Statistical Mechanics: Algorithms and Computations。

㈨ 谈谈分子模拟中的能量最小化,弛豫和平衡态

目前分子模拟在很多领域得到了广泛地应用。这主要得益于分子模拟理论的逐渐成熟和模拟工具的开发及完善。Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 是目前材料领域应用广泛的分子模拟工具包;而Gromacs则是生物体系模拟的一把利器。根据统计热力学的要求,只有对平衡后结构做分析才能得到合理的结果。[1] 平衡态分子模拟的一般运行流程是能量最小化,弛豫、平衡模拟及采样。下面我们将详细介绍这三个流程。

1        能量最小化

分子模拟的初始构型往往偏离平衡态较远。如果初始结构中有些原子靠的比较近,这样系统的能量就会过高而使得模拟崩溃。另外如果原子靠的很近,相互作用力就是非常大,可能会导致原子在一个时间步内移出盒子边界从而出现原子丢失的现象。一般分子动力学软件包都提供了能量最小化(minimize)工具,能够在一定程度上优化结构以减少原子间的坏的接触。除了能量最小化工具,LAMMPS还提供了一些独特的手段处理比较差的初始构形:

① pair_style soft :使用一个余弦对势将重叠的原子排开,这个势能确保原子重叠时,体系的能量不会爆炸。进行一段时间的模拟后,可能会得到比较合理的构型,从而可以用正常的力场参数继续进行模拟。

②  fix nve/limit 0.1 : 限制一个原子在一个时间步内的最大位移,从而避免原子丢失。

当然能量最小化工具只能在一定程度上缓解(alleviate)坏的构型冲搏带来的影响。为了顺利地模拟,在搭建构型时就应注意避开原子靠近的情况。周期性带来的原子重叠问题往往容易被忽略,因此需格外小心。

2        弛豫

弛豫(relaxation)指的是将非平衡的系统慢慢演变为平衡态。[2] 模拟刚开始时,系统的温度,压力等性质可能偏离设定值很远,因此需要经过一段时间的模拟趋于平衡。通常针对体系的不同,可以分几步弛豫不同的物理量。但是最终都需要一个预平衡过程中将系统弛豫到平衡。弛豫通常采用什么系综和控温、控压方法呢?

1       系综的选择

①   如果需要在NPT 系综的平衡态下抽样,可以考虑先将体系在NVT下模拟,把温度T弛豫到稳态,再转NPT系综进行预平衡,从而把压强P弛豫到稳态。因为温度相较于压强容易控制,在NVT系综下模拟能够在短时间内能进入平衡态。

②   如果需要在NVT 系综的平衡态下抽样,可以先在NPT系综下将系统散燃祥密度弛豫到稳态,然后再转NVT进行预平衡。

③   特殊体系处理: (A) 两相体系[1,3,4]:在构建初始构型前,可以对两相分别进行平衡,再组合成新的结构。在NVT系综下,限制固相的位置,让液体弛豫以适应固相的位置。然后放开限制,转到NPT系综下进行预平衡。 (B) 不同分子的尺寸相差较大,如蛋白质在溶液中[1,5]:在NVT系综下,限制蛋白的位置,让溶液弛豫以适应蛋白的位置。然后放开限制,转到NPT系综下进行预平衡。(C)对于平衡时间长的体系,如离子液体[1]:此类体系所需的平衡时间较长,可以考虑在NVT下高温跑一段时间,使体系空间分布变得均匀,然后退火到目标温度并弛豫温度,再转到NPT系综下进行预平衡。

通常的速度控制方法有

速度标度方法简单,但是调节的速度不能严格符合玻尔兹曼(Boltzmann)分布,也就是不属于严格的正则系综。Berendsen是一种弱耦合热浴,在系统远离平衡态时,对温度的调节较好,相较于Nosé-Hoover温度的震荡较小。因此可以先用Berendsen将温度弛豫到设定值,再转用Nosé-Hoover进行预平衡。

速度中等,精确度中等。Andersen和郎之万热浴都引入了非物理的随机力,干扰了系统正常的时间演化,粒子之间的速度相关性减弱, 因此算出的动力学性质是不可靠的,但是可以用来计算热力学性质 。

比较复杂,计算的速度相对慢了些。严格遵守正则系综,体系可以时间反演,通常用于平衡采样。但是在段洞系统远离平衡态时,可能温度涨落巨大,难以收敛,因此可以考虑先使用Berendsen、Langevin热浴先将系统温度弛豫到合理的值,再转Nosé-Hoover热浴进行预平衡、平衡采样。[8]

LAMMPS 中实现Berendsen和Langevin的命令分别为:fix temp/berendsen 和 fix langevin。但是这两种方式只会调节速度或力,并不会积分位置,因此需要配合 fix NVE 更新位置,因此其实系统是NVT正则系综。而实现Nosé-Hoover控温只需要一个fix nvt/npt就可以更新速度、力以及位置。

控压方法的选择原理与控温方法基本一致。Berendsen 控压适合用于最初的压力弛豫,不适用于平衡采样。而基于Nosé-Hoover算法的Nosé-Hoover以及Parrinello-Rahman 恒压器(barostat)适用于平衡态控压[9]。前者常用于LAMMPS,而后者是Gromacs中常用的控温器。

对于各向异性的体系,比如双脂膜,均匀压力缩放(isotropic)并不适合。双脂层在xy平面无限大,而在z方向尺寸有限。因此系统在x,y方向上的压力调节与z方向的调节应该脱钩。在Gromacs中通过设置pcoupltype=semiisotropic实现[9],而在LAMMPS中可使用fix npt 命令中分别指定x,y,z方向的参数并指定couple xy 实现。[10]

Nosé-Hoover控温和控压算法并不是万能的。Nosé-Hoover算法中温度和压力的标度并不是瞬时的,而是在一个关联时间tau的周期内弛豫温度或者压力。LAMMPS推荐的控温关联时间Tdamp为100个时间步(100.0*dt),而控压关联时间Pdamp推荐值为1000个时间步。[10]由于关联时间的存在,在计算当前压力时会有几个步骤用到前几个时间步的压力。因此对于需要获得准确压力波动,或者系统的温度、压力关联时间比较短的情况,Nosé-Hoover算法有可能是不合适的。[9]。

对于平衡态模拟,采样必须在确认系统进入平衡态之后才能进行。预平衡和平衡态模拟的参数(系综,控温控压方式)基本一致,但也可以通过减少预平衡过程数据输出(轨迹,日志等)以加快平衡的速度。平衡并不是意味着数值稳定。温度、压力波动的幅度都非常大,尤其是压力变化可能会出现符号的变化。压力为正说明系统有膨胀的趋势,而压力为负系统有缩小的趋势。因此判断系统是不是处于平衡态需要考察所关注的物理量的平均值是不是收敛在设定值附近。通常可以通过考察系统的势能(potential  energy)、温度(T)、压力、均方根波动(RMSD) [11] 等物理量是不是都达到了收敛从而判定系统是否进入了平衡态。再次强调,数据需要从到达平衡态的时间点开始收集。

采样的一个原则是样本数目足够多[12]。样本量越大,结果的可信度就越高。另一个原则是样本抽样要尽量无关联。合理地选取收集样本的间隔能够降低样本的关联性。[12] 一个经验做法是通过计算速度自相关函数确定大致的弛豫时间(tau)[13],采样的时间尺度至少跨越100*tau,或者采集至少5000个样本。两相体系以及非平衡态模拟的采样请参考教程[3,12]。

平衡态模拟所采用的系综往往是按照自己的模拟需求选用的。一般可以考虑选用NPT系综进行平衡模拟,这样系综的密度能够弛豫到与实验接近的状态。而扩散系数的模拟一般是在NVT系综下,因为NPT系综可能会产生原子坐标标度的问题。如果算一些物性,比如热容,压缩系数之类的,需要在公式对应系综下抽样。[1,14,15]

定容比热容 Cv 定义为恒容条件下能量对温度的偏导数,这就决定了必须在定容的系综下抽样。根据热力学公式可以分别推导得到 Cv 在不同定容系综下的公式[14]:

其中< >代表的是系综平均。

可以看到NVT是计算定容比热容 Cv 最简单的系综,因为温度得到了控制,因此只需要得到势能的涨落即可计算得到 Cv 。

计算定压热容 Cp 可以在NPT系综下进行,需要得到 H+PV 的涨落。