利用PWmat与VASP做分子动力学.docx

利用PWmat与VASP做分子动力学.docx

《利用PWmat与VASP做分子动力学.docx》由会员分享，可在线阅读，更多相关《利用PWmat与VASP做分子动力学.docx（21页珍藏版）》请在冰豆网上搜索。

利用PWmat与VASP做分子动力学

利用PWmat与VASP做分子动力学

一、资料

1、VASP

vasp做分子动力学的好处，由于vasp是近些年开发的比较成熟的软件，在做电子scf速度方面有较好的优势。

缺点：

可选系综太少。

尽管如此，对于大多数有关分子动力学的任务还是可以胜任的。

主要使用的系综是NVT和NVE。

一般做分子动力学的时候都需要较多原子，一般都超过100个。

当原子数多的时候，k点实际就需要较少了。

有的时候用一个k点就行，不过这都需要严格的测试。

通常超过200个原子的时候，用一个k点，即Gamma点就可以了。

主要参数：

INCAR：

EDIFF一般来说，用1E-4或者1E-5都可以，这个参数只是对第一个离子步的自洽影响大一些，对于长时间的分子动力学的模拟，精度小一点也无所谓，但不能太小。

IBRION=0分子动力学模拟

IALGO=48一般用48，对于原子数较多，这个优化方式较好。

NSW=1000多少个时间步长。

POTIM=3时间步长,单位fs,通常1到3.

ISIF=2计算外界的压力.

NBLOCK=1多少个时间步长，写一次CONTCAR，CHG和CHGCAR，PCDAT.

KBLOCK=50NBLOCK*KBLOCK个步长写一次XDATCAR.

ISMEAR=-1费米迪拉克分布.

SIGMA=0.05单位:

电子伏

NELMIN=8一般用6到8,最小的电子scf数.太少的话,收敛的不好.

LREAL=A

APACO=10径向分布函数距离,单位是埃.

NPACO=200径向分布函数插的点数.

LCHARG=F尽量不写电荷密度,否则CHG文件太大.

TEBEG=300初始温度.

TEEND=300终态温度。

不设的话，等于TEBEG.

SMASS=-3NVEensemble；-1用来做模拟退火；大于0NVT系综。

///////////////////////////////////////////////////////////////////////////////////////////////

1）收敛判据的选择

结构弛豫的判据一般有两种选择：

能量和力。

这两者是相关的，理想情况下，能量收敛到基态，力也应该是收敛到平衡态的。

但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大，力收敛速度绝大部分情况下都慢于能量收敛速度。

这是因为力的计算是在能量的基础上进行的，能量对坐标的一阶导数得到力。

计算量的增大和误差的传递导致力收敛慢。

到底是以能量为收敛判据，还是以力为收敛判据呢？

关心能量的人，觉得以能量为判据就够了；关心力相关量的人，没有选择，只能用力作为收敛标准。

对于超胞体系的结构优化，文献大部分采用Gamma点做单点优化。

这个时候即使采用力为判据（EDIFFG=-0.02），在做静态自洽计算能量的时候，会发现，原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。

这个时候，是不是该怀疑对超胞仅做Gamma点结构优化的合理性呢？

是不是要提高K点密度再做结构优化在我看来，这取决于所研究的问题的复杂程度。

我们的计算从原胞开始，到超胞，到掺杂结构，到吸附结构，到反应和解离。

每一步都在增加复杂程度。

结构优化终点与初始结构是有关的，如果遇到对初始结构敏感的优化，那就头疼了。

而且，还要注意到，催化反应不仅与原子本身及其化学环境有关，还会与几何构型有关。

气固催化反应过程是电子的传递过程，也是分子拆分与重新组合的过程。

如果优化终点的构型不同，可能会导致化学反应的途径上的差异。

仅从这一点来看，第一性原理计算的复杂性，结果上的合理性判断都不是手册上写的那么简单。

对于涉及构型敏感性的结构优化过程，我觉得，以力作为收敛判据更合适。

而且需要在Gamma点优化的基础上再提高K点密度继续优化，直到静态自洽计算时力也是达到收敛标准的。

（2）结构优化参数设置

结构优化，或者叫弛豫，是后续计算的基础。

其收敛性受两个主要因素影响：

初始结构的合理性和弛豫参数的设置。

初始结构

初始结构包括原子堆积方式，和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。

比如掺杂，表面吸附，空位等结构，初始原子的距离，角度等的设置需要有一定的经验积累。

DFT计算短程强相互作用（相对于德华力），如果初始距离设置过远（如超过4埃），则明显导致收敛很慢甚至得到不合理的结果。

比较好的设置方法可以参照键长。

比如CO在O顶位的吸附，可以参照CO2中C-O键长来设置（如增长20%）。

也可以参照文献。

记住一些常见键长，典型晶体中原子间距离等参数，有助于提高初始结构设置的合理性。

实在不行，可以先在小体系上测试，然后再放到大体系中算。

弛豫参数

弛豫参数对收敛速度影响很大，这一点在计算工作没有全部铺开时可能不会觉察到有什么不妥，反正就给NSW设置个“无穷大”的数，最后总会有结果的。

但是，时间是宝贵的，恰当的设置3小时就收敛的结果，不恰当的设置可能要一个白天加一个黑夜。

如果你赶文章或者赶着毕业，你就知道这意味这什么。

结构优化分电子迭代和离子弛豫两个嵌套的过程。

电子迭代自洽的速度，有四个影响很大的因素：

初始结构的合理性，k点密度，是否考虑自旋和高斯展宽（SIGMA）；离子弛豫的收敛速度，有三个很大的影响因素：

弛豫方法（IBRION）,步长（POTIM）和收敛判据（EDIFFG）.

一般来说，针对理论催化的计算，初始结构都是不太合理的。

因此一开始采用很粗糙的优化（EDIFF=0.001，EDIFFG=-0.2），很低的K点密度（Gamma），不考虑自旋就可以了，这样NSW<60的设置就比较好。

其它参数可以默认。

经过第一轮优化，就可以进入下一步细致的优化了。

就我的经验，EDIFF=1E-4,EDIFFG=-0.05，不考虑自旋，IBRION=2，其它默认，NSW=100;跑完后可以设置IBRION=1，减小OPTIM（默认为0.5，可以设置0.2）继续优化。

优化的时候让它自己闷头跑是不对的，经常看看中间过程，根据情况调节优化参数是可以很好的提高优化速度。

这个时候，提交两个以上的任务排队是好的方式，一个在调整的时候，下一个可以接着运行，不会因为停下当前任务导致机器空闲。

无论结构优化还是静态自洽，电子步的收敛也常常让新手头痛。

如果电子步不能在40步收敛，要么是参数设置的问题，要么是初始模型太糟糕。

静态自洽过程电子步不收敛一般是参数设置有问题。

这个时候，改变迭代算法（ALGO），提高高斯展宽（SIGMA增加）,设置自洽延迟（NELMDL）都是不错的方法。

对于大体系比较难收敛的话，可以先调节AMIN,BMIX跑十多步，得到电荷密度和波函数，再重新计算。

实在没办法了，可以先放任它跑40步，没有收敛的迹象的话，停下来，得到电荷密度和波函数后重新计算。

一般都能在40步收敛。

对于离子弛豫过程，不调节关系也不大。

开始两个离子步可能要跑满60步（默认的），后面就会越来越快了。

（3）优化结果对初始结构和“优化路径”的依赖

原子吸附问题不大，但是小分子吸附，存在初始构型上的差异。

slab上水平放置，还是垂直放置，可能导致收敛结果上的差异。

根据H-K理论，理想情况下，优化得到的应该是全局最小，但在数值计算的时候可能经常碰到不是全局最小的情况。

实际操作中发现，多个不同初始结构优化收敛后在能量和结构上存在一定差异。

为了加快收敛速度，特别是对于表面-分子吸附结构，初始放松约束，比如EDIFF=1E-3，EDIFFG=-0.3,NSW=30可能是很好的设置。

但是下面的情况应当慎重：

EDIFF=1E-3;

EDIFFG=-0.1;！

或者更小

NSW=500；！

或者更大

电子步收敛约束较小，而离子步约束偏大，离子步数又很多，这种情况下，可能导致的结果是结构弛豫到严重未知的区间。

再在这个基础上提高约束来优化，可能就是徒劳的了——结果不可逆转的偏向不正常的区间。

好的做法，是对初始结构做比较松弛的约束，弛豫离子步NSW应该限制在一个较小的数值。

EDIFF=1E-3的话，EDIFFG也最好是偏大一些，如-0.3而不是-0.1.这样可以在较少的步数达到初步收敛。

对于远离基态的初始结构，一开始在非常松弛的约束下跑若干离子步，时间上带来的好处是很大的。

对于100个原子的体系用vasp做Gamma点优化，如果一开始就是正常优化（EDIFF=1E-4,EDIFFG=-0.02）设置，开始十个离子步可能都要花上几个小时。

如果这个时候才发现输入文件有错误，那下午的时间就白费了，顺便带上晚上机器空转。

所以，我习惯的做法，是在初始几步优化后，会用xcrysden检查一下XDATCAR中的数据，用xdat2xyz.pl生成movie.xyz，然后看看弛豫过程是不是按照设想的那样。

后续过程跑完一个收敛过程，就再检查一下movie.xyz。

如此这般，才放心的展开后续计算。

//////////////////////////////////////////////////////////////////////////////////////////////////////

结构完全相同的、处于各种运动状态的、各自独立的系统的集合。

全称为统计系综。

系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式；系综并不是实际的物体，构成系综的系统才是实际物体。

研究气体热运动性质和规律的早期统计理论是气体动理论。

统计物理学的研究对象和研究方法与气体动理论有许多共同之处，为了避免气体动理论研究中的困难，它不是以分子而是以由大量分子组成的整个热力学系统为统计的个体。

系综理论使统计物理成为普遍的微观统计理论。

系统的一种可能的运动状态，可用相宇中的一个相点表示，随着时间的推移，系统的运动状态改变了，相应的相点在相宇中运动，描绘出一条轨迹，由大量系统构成的系综则可表为相宇量相点的集合，随着时间的推移，各个相点分别沿各自的轨迹运动，类似于流体的流动。

系综并不是实际的物体，构成系综的系统才是实际物体。

约束条件是由一组外加宏观参量来表示。

在平衡统计力学畴下，可以用来处理稳定系综。

一、常用系综分类

根据宏观约束条件，系综被分为以下几种：

1.正则系综（canonicalensemble），全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数（N）、体积（V）、温度（T）。

正则系综是蒙特卡罗方法模拟处理的典型代表。

假定N个粒子处在体积为V的盒子，将其埋入温度恒为T的热浴中。

此时，总能量（E）和系统压强（P）可能在某一平均值附近起伏变化。

平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。

正则系综的特征函数是亥姆霍兹自由能F（N,V,T）。

2.微正则系综（micro-canonicalensemble），简写为NVE，即表示具有确定的粒子数（N）、体积（V）、总能量（E）。

微正则系综广泛被应用在分子动力学模拟中。

假定N个粒子处在体积为V的盒子，并固定总能量（E）。

此时，系综的温度（T）和系统压强（P）可能在某一平均值附近起伏变化。

平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。

微正则系综的特征函数是熵S（N,V,E）。

3.等温等压（constant-pressure,constant-temperature），简写为NPT，即表示具有确定的粒子数（N）、压强（P）、温度（T）。

一般是在蒙特卡罗模拟中实现。

其总能量（E）和系统体积（V）可能存在起伏。

体系是可移动系统壁情况下的恒温热浴。

特征函数是吉布斯自由能G（N,P,T）。

4.等压等焓（contant-pressure,constant-enthalpy），简写为NPH，即表示具有确定的粒子数（N）、压强（P）、焓（H）。

由于由于H=E+PV，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少遇到了。

5.巨正则系综（grandcanonicalensemble），简写为VTμ，即表示具有确定的粒体积（V）、温度（T）和化学势（μ）。

巨正则系综通常是蒙特卡罗模拟的对象和手段。

此时、系统能量（E）、压强（P）和粒子数（N）会在某一平均值附近有一个起伏。

体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的T，。

特征函数是马休（Massieu）函数J（μ，V，T）。

二、系综调节

系综调节主要是指在进行分子动力学计算过程中，对温度和压力参数的调节，分为调温技术和调压技术。

1.调温技术

在NVT系综或NPT系综中，即使在NVE系综模拟的平衡态中，也经常调整温度到期望值。

如果希望知道系统的平衡态性质怎样依赖于温度，那么就必须在不同的温度下进行模拟。

目前实现对温度的调节有4种方式：

速度标度、Berendsen热浴、Gaussian热浴、Nose—Hoover热浴。

2.调压技术

在等压模拟中，可以通过改变模拟原胞的三个方向或一个方向的尺寸来实现体积的变化．类似于温度控制的方法，也有许多方法用于压力控制，总的来说有以下3种技术：

Berendsen方法、Anderson方法、Parrinello-Rahman方法。

三、系综选择

原则上巨正则系综应用最广，但却不一定是最方便的。

因为可以看到三种系综的演化过程既是约束解除的过程，却是以增加变量为代价的，这也就增加了数学上的复杂性。

因此一般情况下如果不需求解。

，则不必使用巨正则系综。

系综选择的基本原则为：

1.微正则系综能够简单的求得近独立，全同，定域粒子系统，并且每个粒子只能有两个不同的可能状态，例如简单的铁磁，顺磁模型

2.微正则系综难求的系统，可用正则系综求解

3.当微正则和正则系综均难求时，可用巨正则系综求解

系综（ensemble）：

在一定的宏观条件下，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统的集合。

全称为统计系综。

系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式。

以下是介绍系综的容。

系综的性质

研究对象

常用的三个系综

原理

系综的性质：

系综是假想的概念，并不是真实的客观实体。

真正的实体是组成系综的一个个系统，这些系统具有完全相同的力学性质。

每个系统的微观状态可能相同，也可能不同，但是处于平衡状态时，系综的平均值应该是确定的。

研究对象：

研究气体热运动性质和规律的早期统计理论是气体动理论。

统计物理学的研究对象和研究方法与气体动理论有许多共同之处，为了避免气体动理论研究中的困难，它不是以分子而是以由大量分子组成的整个热力学系统为统计的个体。

系综理论使统计物理成为普遍的微观统计理论。

常用的三个系综：

J.W.吉布斯把整个系统作为统计的个体，提出研究大量系统构成的系综在相宇中的分布，克服了气体动理论的困难，建立了统计物理。

在平衡态统计理论中，对于能量和粒子数固定的孤立系统，采用微正则系综（NVE）；对于可以和大热源交换能量但粒子数固定的系统，采用正则系综（NVT）；对于可以和大热源交换能量和粒子的系统，采用巨正则系综（mVT）。

这是三种常用的系综，各系综在相宇中的分布密度函数均已得出。

量子统计与经典统计的研究对象和研究方法相同，在量子统计中系综概念仍然适用。

区别在于量子统计认为微观粒子的运动遵循量子力学规律而不是经典力学规律，微观运动状态具有不连续性，需用量子态而不是相宇来描述。

原理

系统的一种可能的运动状态，可用相宇中的一个相点表示，随着时间的推移，系统的运动状态改变了，相应的相点在相宇中运动，描绘出一条轨迹，由大量系统构成的系综则可表为相宇量相点的集合，随着时间的推移，各个相点分别沿各自的轨迹运动，类似于流体的流动。

　　若系统具有s个自由度，则相宇是以s个广义坐标p（详写为p、p2……ps）和s个广义动量q（详写为q1、q2……qs）为直角坐标构成的2s维空间。

在相宇任一点（p，q）附近单位相体积元的相点数目D（p，q，t）称为密度函数。

D（p，q，t）在整个相宇的积分等于全部相点数，即等于系综所包含的全部系统数N，与时间t无关。

定义ρ（p,q,t）=D（p,q,t）/N，称为系综的概率密度函数。

ρ（p，q，t）dpdq表示在t时刻出现在（p，q）点附近相体积元dpdq的相点数在全部相点数中所占的比值，即表示任一系统在t时刻其运动状态处于（p，q）附近的相体积元dpdq的概率。

显然，概率密度函数ρ（p，q，t）满足归一化条件∫ρ（p,q,t）dpdq=1。

　　统计物理学的认为系统的任意宏观量I（t）是相应微观量L（p，q）在一定宏观条件下对系统一切可能的微观运动状态的统计平均值，即I（t）＝∫L（p，q）ρ（p，q，t）dpdq。

由此可见，经典统计物理的基本课题是确定各种条件下系综的概率密度函数ρ（p，q，t），ρ确定后，即可对相应的热力学系统的宏观性质作出统计描述。

这就是统计系综的方法。

　　ρ（p，q，t）的具体形式与系统所处的宏观状态有关。

如果系统处于平衡态，则ρ＝ρ（p，q）不显含时间t，在平衡态的系综理论中，由能量、体积和粒子数都固定的系统构成的统计系综称为微正则系综；由与温度恒定的大热源接触，具有确定粒子数和体积的系统构成的统计系综称为正则系综；由与温度恒定的大热源和化学势恒定的大粒子源接触，具有确定体积的系统构成的统计系综称为巨正则系综；由与温度恒定的大热源接触并通过无摩擦的活塞与恒压强源接触，具有确定粒子数的系统构成的统计系综称为等温等压系综。

上述各种统计系综都有各自的概率密度函数。

在微正则系综中，系统处于所有可能的微观状态上的概率都相等，即概率密度是不随时间改变的常数，这就是等概率原理。

等概率原理是平衡态统计物理的基本假设，它的正确性由它的推论与实际相符而得到肯定。

由微正则系统可以推导出其它系综的概率分布函数的形式。

微正则系综是由许多具有相同能量，粒子数，体积的体系的集合。

它是统计力学系综的一种。

其配分函数Ω是在能量E_0上的能态密度。

微正则系综是个简并度下的正则系综，正则系综可以被分开进入子系综，每个子系综被对应到可能的能量值且自身为另一些微正则系综。

很长一段时间以来，一直有人在问BOMD与CPMD究竟有什么不同？

这里我就简单说一些二者的区别。

实际上，我想大部分人理解和接触的第一原理分子动力学方法以CPMD居多。

CPMD，就是Car和Parrinello两个人作出的基于密度泛函的分子动力学方法，其特点是在引入电子虚拟质量，将电子运动耦合到了运动方程中，每一步分子动力学计算后，对电子结构的计算就不再需要自洽场迭代这一过程，因此可以大大节省计算资源，在计算机还不是那么好的80，90年代是弥足珍贵的。

CPMD的建立开创了第一原理分子动力学方法的新时代，使其真正开始了实用化进程。

而BOMD，顾名思义，就是Born-Oppenheimer分子动力学。

Born-Oppenheimer近似，也就是绝热近似，指的是将电子和离子的求解分离开来，只处理离子的动力学部分，而认为电子可以快速跟上电子的运动。

其特点是每一步分子动力学计算之后都需要对电子结构进行自洽场迭代，使电子达到基态。

正是由于这个自洽场迭代过程需要的计算量巨大，致使其一直没有达到广泛运用，直到90年代中后期，计算机技术的发展，才使BOMD开始逐步被人们重视。

CPMD和BOMD各有优劣。

CPMD虽然计算速度更快，不需要进行自洽场迭代计算，计算量小，但是由于计算中并没有使体系真正达到基态，而只是尽量靠近基态，因此其准确性对电子的虚拟质量这个参数的选取依赖程度很大，一旦计算参数不对，得到的体系很可能远远偏离真实的势能面，得不到正确的动力学轨迹。

同时，为了保证其尽量靠近基态，分子动力学时间步长一般选得较小。

而BOMD虽然计算量大，但是由于每一步都保证系统达到基态，因此其分子动力学步长可以取得较大，一般1fs到5fs都有可能。

综合来说，如果能结合二者的优势，第一原理分子动力学计算效率将大大提高，这也成为近年来其发展的重要方向。

一般做分子动力学的时候都需要较多原子，一般都超过100个。

当原子数多的时候，k点实际就需要较少了。

有的时候用一个k点就行，不过这都需要严格的测试。

通常超过200个原子的时候，用一个k点，即Gamma点就可以了。

INCAR：

EDIFF一般来说，用1E-4或者1E-5都可以，这个参数只是对第一个离子步的自洽影响大一些，对于长时间的分子动力学的模拟，精度小一点也无所谓，但不能太小。

IBRION=0分子动力学模拟

IALGO=48一般用48，对于原子数较多，这个优化方式较好。

NSW=1000多少个时间步长。

POTIM=3时间步长,单位fs,通常1到3.

ISIF=2计算外界的压力.

NBLOCK=1多少个时间步长（离子步），写一次CONTCAR，CHG和CHGCAR，PCDAT.

KBLOCK=50NBLOCK*KBLOCK个步长写一次PCDAT

ISMEAR=-1费米迪拉克分布.

SIGMA=0.05单位:

电子伏

NELMIN=8一般用6到8,最小的电子scf数.太少的话,收敛的不好.

LREAL=A

APACO=10径向分布函数距离,单位是埃.

NPACO=200径向分布函数插的点数.

LCHARG=F尽量不写电荷密度,否则CHG文件太大.

TEBEG=300初始温度.

TEEND=300终态温度。

不设的话，等于TEBEG.

SMASS-3NVEensemble;-1用来做模拟退火。

大于0NVT系综。

SMASS=1,2,3是没有区别的。

都是NVTensemble。

SMASS只要是大于0就是NVT系综。

CONTCAR是每个离子步之后都会写出来的，但是会用新的把老的覆盖

CHG是在每10个离子步写一次，不会覆盖

CHGCAR是在任务正常结束之后才写的。

1、分子动力学模拟是不是分为classicMolecularDynamics和abinitiomoleculardynamics？

这两种MD有什么区别呢？

classical的力和能量是直接给的原子距离,键角的表达式abinitio的是从第一性原理计算得到的,解Kohn-Sham方程

2、LAMMPS中MD是classicMD，计算的结果可靠吗？

计算量大不大？

如果计算100个原子左右的体系，需要的存和核数大概为多少？

看你的力场,力场好就准,100个原子体系对于classical很小,笔记本随便跑,对abinitio很大,即使是CH系也要上服务器,一般估测16-32G存可以,核心越多越好,8线程至少吧

3、vasp软件中的MD应该是abinitiomoleculardynamics吧？

但看到VASP手册说计算结果不一定准确，所以很纠结要不要用vasp算。

而且好像vasp做MD计算所需计算量蛮大的。

VASPMD的NVT模块没什么问题,NPT不行

VASP的MD计算量很大，

二、计算

1、VASP做分子动力学

按照上述资料，我选了GeS，即2Ge+2S的一个小系统试处理。

按照结构优化、静态自洽、分子动力学这三步进行。

结构优化参数设置如下：

其中SYSTEM表示该体系及其计算命名；ISTART用于设置从头开始计算抑或是从上次结果开始计算，这里设置从头开始计算；ICHARG=2表示如何建立初始电荷密度；PREC表示精度；LREAL=F表示投影于倒易空间；IBRION=2决定原子如何移动，2是CG方法（共轭梯度算法）；ISIF=2表示计算离子所受力、计算原胞的应力量、粒子位置驰豫。

NSW表示最多进行100步。

POTIM表示离子移动步长？

；EDIFF、EDIFFG两个参数表示离子驰豫中断条件；ISMEAR=0表示是半导体材料；S