CP2K Tutorial - MD
这个教程是我上学期为统计力学准备的习题
原文链接exercises:2017_uzh_acpc2:l-j_flu [CP2K Open Source Molecular Dynamics ]
推荐书目
1. Understanding Molecular Simulation
2.Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods
书目1的作者是Daan Frenkel, Berend Smit，每年冬天在阿姆斯特丹都会有CECAM-Winter School, 感兴趣的同学可以去听一下，同时会送这本书。
书目2的作者是我老板Hutter，该书比较难懂，建议先读一些第一性原理的和第一本书。
本文的例子使用Ar原子的L-J potential来模拟Ar的流体性质。
Lennard-Jones_potential的相关介绍请在Wiki中自行阅读，简单说就是使用12次方相和6次方相来表示引力和斥力，这些当然也可以通过DFT计算得到。本文主要侧重于MD，具体的AIMD的从头计算法部分，可以通过改变计算中如何得到力的方式来进行不同精度的计算。
## It's highly recommended to go
## CP2K Input Reference Manual
## and learn how to set up CP2K
## calculation correctly using manual.
&GLOBAL
PROJECT ar108 #Project Name
RUN_TYPE md #Calculation Type : MD (molecular dynamics), GEO_OPT (Geometry Optimization), Energy (Energy Calculation)
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE #The ensemble for MD propagation, NVE (microcanonical), NVT (canonical), NPT_I (NPT with isotropic cell)
STEPS 100 #The number of MD steps to perform
TIMESTEP 5. #The length of an integration step (fs)
TEMPERATURE 85.0 #The temperature in K used to initialize the velocities with init and pos restart, and in the NPT/NVT simulations
&END MD
&END MOTION
&FORCE_EVAL
METHOD FIST #Method to calculate force: Fist (Molecular Mechanics), QS or QUICKSTEP (Electronic structure methods, like DFT)
&MM
&FORCEFIELD
&CHARGE #charge of the MM atoms
ATOM Ar #Defines the atomic kind of the charge
CHARGE 0.0 #Defines the charge of the MM atom in electron charge unit
&END
&NONBONDED
&LENNARD-JONES #LENNARD-JONES potential type.Functional form: V(r) = 4.0 * EPSILON * [(SIGMA/r)^12-(SIGMA/r)^6]
atoms Ar Ar #Defines the atomic kind involved in the nonbonded potential
EPSILON 119.8 #Defines the EPSILON parameter of the LJ potential (K_e)
SIGMA 3.405 #Defines the SIGMA parameter of the LJ potential (Angstrom)
RCUT 8.4 #Defines the cutoff parameter of the LJ potential
&END LENNARD-JONES
&END NONBONDED
&END FORCEFIELD
&POISSON # Poisson solver
&EWALD
EWALD_TYPE none
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 17.1580 17.158 17.158 #Simulation box size
&END CELL
&COORD
Ar -8.53869012951987116 -15.5816257770688615 2.85663672298278293
Ar 1.53007304829383051 9.28528179040142554 11.1777824543317941
Ar 11.9910225119590699 -7.48825329565798015 -9.96545306345559823
Ar -12.6782400030290496 -3.34105872014234606 4.07471097818485806
Ar -1.77046254278594462 -0.232459464264201887 13.2012946017273016
Ar 8.01761371186688443 -2.57249587730733298 -4.12720554747711432
Ar 8.57849517232300052 4.01396664624232002 5.57368821983998419
Ar -3.89200679277030925 -10.2930917801117356 -6.98640232289045482
Ar -3.35457160564444568 -16.1119619276890056 16.1358515626317427
Ar 9.78957155103081966 -16.2628264194939263 -5.69790857071688350
Ar 0.505143495414835719 -4.22978415759568183 12.4854171634357307
Ar 15.5632243939617503 -7.98048905093276240 2.20994708545912832
Ar -5.40741643995084953 -2.64764457113743079 -0.681485212640798199
Ar -0.983719068448489081E-01 -1.73674004862212694 -7.11915545117132265
Ar 7.52655781331927187 -5.52969969672439632 -12.8886150439489313
Ar -5.45655410995716128 0.564445754429787061 2.03902510096247536
Ar -11.8590998267164665 3.40407446386207724 3.72687933934436399
Ar 16.7175362589401821 -7.47132377347522780 -1.02274476672697889
Ar -20.4572129717055340 -5.73700807719791683 4.81845086375497811
Ar 14.8485522289272627 -1.41608633045414667 -16.0839111490847451
Ar 8.04379470511429595 -8.14033814842439263 -4.75543123809189261
Ar 12.2738439612049568 -1.70589834674486429 12.9622486199573572
Ar -0.421851806372696092 -11.1177490353157999 20.4545363332536283
Ar 2.28194341698637571 5.92083917539752136 -11.1732449877738436
Ar -13.9648466918215064 8.77923885764231926 8.07373370482465091
Ar -10.3147439499058429 6.38529561240966004 -15.3411964215061527
Ar -2.71899964647918457 -21.4890074469143855 10.8678096818980006
Ar -17.7923879123397271 -10.7840901151121251 -4.83954996524571968
Ar 5.23494138507746420 -6.79222906792632841 -6.07187690814296133
Ar 3.52448750638480446 -10.1225951872349782 2.96829048662758721
Ar -16.1586602901979361 -5.18274316385346445 8.57072694078649455
Ar -5.80982824422251287 4.32640193501643733 2.55599101868223322
Ar 6.29160109084684382 27.5741337288405717 15.0246410590392632
Ar -3.18741711710350684 23.2996469099840624 -16.8034854143018748
Ar -4.20225755039435622 9.36037725943080190 16.5891306154890081
Ar -7.64392908749747946 -9.52432384411045341 -29.8228731471089645
Ar 0.545352525792712428 13.9240554617015260 -0.383786780333776500
Ar -5.27432886808646906 -5.53813781787395865 -20.3014703747109415
Ar 22.9921850152838871 6.78619371666398941 -1.98289905290632484
Ar 19.7720034229251880 -10.2373337687313679 -3.33081818566269172
Ar 0.156776902886395425 6.59630118110908725 8.90749062505743083
Ar 5.57937381862174053 0.233106223140015806 1.02752287819280941
Ar -3.64343561800208793 3.96448881012491006 25.8752124557059595
Ar -0.248491698112870391 20.4489725648023182 -2.51220445353457666
Ar 2.93626708600658270 0.859812213376437984 9.96743307236779508
Ar 3.30384315693043895 -2.92421266591109408 -6.34927042371499883
Ar -6.15490235244551265 -6.84961480075890883 -6.46204144605644260
Ar -23.2388291761596619 -28.1213094673208666 7.13721047187827917
Ar 4.11526291325474780 2.71564143367947342 -0.852030043744060328
Ar 14.6194148692240713 2.80815182256426210 1.93601975975151541
Ar 18.9667954753247869 16.5700888519293095 13.3423444868082761
Ar -28.6124161416877705 2.84353637083477562 -9.23601973326721648
Ar -5.97004594556101331 -16.2230172568109978 -9.22928061840017477
Ar 10.0481077882725955 16.3854819569745231 5.12578711346205651
Ar -7.22508507825336643 6.34615422233080650 -0.680757463730119028
Ar -12.0138912984383506 -10.4653110276797570 -6.43434787584580103
Ar -8.53169926903037457 12.8976589212818862 -0.890361252446473683
Ar -25.3700692950848676 3.33119906434656077 -7.93917685683272722
Ar 2.90163480643285920 -12.9668181360039672 -7.94907759259854707
Ar -13.5963940986222074 11.9896580951935974 -3.55068754869933789
Ar 13.1416029517342476 4.97143783446568488 -3.50841252726170705
Ar -13.3295460955805254 16.0410015777677764 7.05282797577515375
Ar -4.15068335494176122 -19.5111913798076593 21.0255971827539376
Ar -8.42944270819351793 16.3065160593537009 -18.2887817284733885
Ar 0.788636333898691255 9.59016836817029095 22.1772606194495872
Ar -2.92606778628861974 -7.97408054890791007 -21.3519900334304964
Ar -6.39959865978756426 -4.56280461803643256 2.75533571094951402
Ar -4.04423878093174860 -14.9275965394452506 -5.58561473738824965
Ar -27.8524912514281482 0.802052180719123098 -3.02663789713126441
Ar -2.83966529645897792 7.11627121253196915 6.18547332762273783
Ar -8.68327887612401739 -6.67088300493855879 -9.15815450219801264
Ar -11.9620847111198501 -2.20956249614563038 -1.83979975374852245
Ar 22.6848553724304907 12.2047209420099971 1.01238797839832362
Ar 6.29501012040417507 -0.769712471349173866 -6.91454332254278281
Ar 3.49995546789933476 -8.00704920137973453 -0.426526631939732892
Ar 0.385154812289867643 17.8769740351009112 -17.4065226240143041
Ar 21.2288869131365736 10.2327102035561044 -13.0872200859088803
Ar 1.22082587001210396 5.83597435065779457 16.8450099266840283
Ar -7.08754036219628425 6.03412971863339109 -22.3251445579668015
Ar -0.244265849036998037E-01 17.4693605251376454 7.37116730966604194
Ar 15.0981822679441553 9.88940516251130397 -8.49382740142986670
Ar -6.57877688336587152 -15.0484532074656290 14.7230359830473887
Ar -2.22666666633409394 -4.18421900331013674 -2.47007887105670587
Ar 5.20621069851729867 -22.6565181989138011 7.39475674805799521
Ar -8.85828800414884299 -2.47510661993999781 2.35441398531938617
Ar 6.75202354538700167 0.430391383628436597 5.43492495261394382
Ar 11.9263127546080856 8.13267254152258445 2.40081132956567966
Ar -14.5507562394484040 -0.471540677239574602 -13.7058431104765983
Ar 14.1157692422228553 -2.98968593175088149 24.6842798176059546
Ar -3.35107336204723527 -0.681362546744063047 -7.37039916831594510
Ar 7.79269876443546838 3.30687615091469800 -0.732378021069576002E-01
Ar -1.13289059102623746 -17.1672835835708497 -12.9126466371968966
Ar -9.21054349522787241 -10.4846510042527843 -8.38485797788161591
Ar -6.47848777956778044 -3.90736653076878993 -10.6499668409808841
Ar 0.987874979233200667 13.7363585340729077 5.07209659800543733
Ar 8.86097814789463278 9.96103887786039799 1.09373795795780060
Ar -6.58068766844202013 -20.4019345282015756 -7.28935608176262662
Ar -0.448977062720621045 19.5862520159664086 11.0351198968750293
Ar 7.36056937465398153 -2.69594281683156067 7.26081874603436361
Ar 13.8791344546872004 12.1903465249438128 1.24889885444881155
Ar -3.65782753722175302 19.7829061761924159 -11.8161510229542408
Ar 1.49729450944005649 -5.39289977250827679 1.92445849672255198
Ar 18.5861605633917577 3.00868366398259690 2.06440131010935168
Ar 6.20730767975507014 -9.47418398815358032 5.54930507752316249
Ar 3.65054837888884753 3.43181054126032858 -4.31160813615129435
Ar 2.67862616463048520 2.29300605146530545 5.98502962150055051
Ar 24.5113122275150914 4.00733170976478448 13.1412501215423774
Ar -0.600233262008137092 3.62825631372324597 6.38411284716526772
&END COORD
&TOPOLOGY
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
如上面的INPUT所示，在&MOTION的部分，我们使用规定了MD的计算参数，在&FORCE_EVAL的部分，我们规定了使用L-J potential来计算。如果你希望使用AIMD来计算，完全可以将修改&FORCE_EVAL中的具体细节，具体请看CP2K Tutorial - DFT and beyond重的介绍，你可以使用不同泛函，杂化泛函，甚至MP2来计算你的trajectory。
在&MOTION中的&MD部分，我们的ENSEMBLE使用了NVE，总共计算100步，步长5fs，起始温度85K。
总共的step，没有什么好说的，就是算多少步。关于timestep，理论上当然是越小越好，这样可以减少numerical error，因为如果步长太大，很明显，控制MD的newtonian的精度会降低，具体MD的算法请参见书目1中的介绍。
下面简单说一下Ensembles这一部分
在ENSEMBLE中，我们可以使用
1. NVE，这个严格遵循newtonian的规则
2. NVT，这个使用了THERMOSTAT
3. NPT_I, 这个同时使用了THERMOSTAT和BAROSTAT，同时不改变晶胞形状比例
更多的ENSEMBLE请参见manual
为什么要引入THERMOSTAT或者BAROSTAT？
因为如果我们严格遵守Newtonian，我们其实给的温度会在一段时间后收敛于某一个值，这个温度很可能不是我们想sample的温度区间，引入THERMOSTAT可以将温度稳定在某一个位置。至于BAROSTAT也是同样的道理，我们想计算某一压力下的性质，很可能我们的温度，原子多少，晶胞大小和当前的压力并不匹配，我们需要对晶胞的体积进行放缩。
下面给出一段NVT计算中使用Nose-Hoover Thermostat的input
&MD
ENSEMBLE NVT
STEPS 3000
TIMESTEP 5
TEMPERATURE 298
&THERMOSTAT
REGION MASSIVE
&NOSE #Uses the Nose-Hoover thermostat
TIMECON 1000 #timeconstant of the thermostat chain, how often does thermostat adjust your system
&END NOSE
&END
&END MD
关于TIMECON的选择，理论上当然是越大越好，因为TIMECON越大，我们对体系的扰动就越小，但是我们想到equilibrium的时间也就越长，如果使用比较aggressive的thermostat，也很可能过分的adjust我们的系统，使得认为造成的误差增大。因此根据实际需要来选择一个合适的TIMECON，一般建议使用100-1000的范围。
请问REGION MASSIVE与REGION GLOBAL有什么区别呢，您这个例子里可以用REGION GLOBAL吗？
Massive顾名思义，这个thermostat更加aggressive，在算法上相当于在每一个自由度上使用了一个小的thermostat从而组成一个大的thermostat
Global也就是整个系统只有一个thermostat
使用massive thermostat可以加快平衡的过程，当然也引入更多的人为误差。
至于使用哪个根据具体情况测试使用，道理和开始使用较小的TIMECON是一个道理，但是最后建议都使用不太aggressive的thermostat。
对于NPT的ensemble，同样给出一个input
&FORCE_EVAL
...
STRESS_TENSOR ANALYTICAL
...
&END FORCE_EVAL
...
...
&MD
ENSEMBLE NPT_I #constant temperature and pressure using an isotropic cell
STEPS 3000
TIMESTEP 5.
TEMPERATURE 85.0
&BAROSTAT
PRESSURE 0. # PRESSURE, unit[bar]
TIMECON 1000
&END BAROSTAT
&THERMOSTAT
&NOSE
TIMECON 1000
&END NOSE
&END MD
使用不同的ENSEMBLE来平衡你的系统，计算Radial distribution functions (RDF), 和原文链接中进行比较RDF
对于BAROSTAT，我们需要在FORCE_EVAL中计算STRESS_TENSOR，从而评价3D晶胞的压力情况来update你的晶胞大小。
Equilibration
为了使得我们的数据有意义，我们计算的系统必须在平衡状态，而不是某一状态。当然如果我们在transtion path sampling，这个另当别论。
一般来说，我们使用optimized geometry来作为起始状态。
对于液体系统，如果你的密度和压力没有问题，原则上先做NVT，然后NVE
对于压力可能存在问题的体系，先做NVT然后NPT，然后NVT，最后NVE
当然也可以直接做NPT，然后NVE，这个根据具体系统达到平衡的速度来看。
当然我们也可以只使用NVE来平衡我们的系统，在&MD中有一个关键词TEMP_TOL，可以设置一个WINDOW来使系统的温度在某一个区间oscillatie。
还有一个方法就是使用MOTION%MD%ENSEMBLE%LANGEVIN来进行平衡，里面的GAMMA建议使用0.001 [1/fs]。关于Langevin MD更具体的内容，我们会在后文中的second generation Car-Parrinello Molecular Dynamics讲到。
