|
Skelton's tips: Phonons & Phonopy: - Atomic structure optimization : to obtain optimized structures (update: 2013-08-24)
- In many cases, the atomic positions are drifted if one use default option for structure relaxations. In other words, atomic forces are hardly reaching small values due to the uncertainty in calculation. - To obtain optimized geometry, it 9is suggested to use the following tag in INCAR. - Or you can increase calculation precision by lowering "EDIFF" or by increasing "NELMIN" It will increase the precision of the force and thereby you wil obtain the optimized structure. NELMIN = 6 - Importance of k-space projection (LREAL=FALSE)
- Not same energies for same structures, LREAL=Auto or LREAL=TRUE - Sometimes, the same structures can give different total energies if either LREAL=TRUE or LREAL=Auto are used. For example, one makes same structure but the atomic basis are equally displaced by constant vector. Many cases it is ok, because the energy difference is small, less than few meV/atom. However, when one calculates the formation energy, the error can be large if the supercell size is very big. If it happens, you should re-calculate the total energy using k-space projection scheme. (LREAL = False) - When you have to calculate very precise calculation for phonon mode or for small barrier calculations, it is strongly recommended to use the "LREAL=F" tag, with other high precision settings even though the supercell size is very large. - Unit of density-of-states (DOS) in DOSCAR file (update: 2013-08-24)
- The unit of DOS is the states per electrons. If you run the single interation, you will always obtain "states/eV" - However, sometimes you will get different unit in DOSCAR. In many cases, it happens becasue you run the molecular dynamics simulation or structure relaxation. Then the DOSCAR contatins the sum or average value of DOSs. Note that the NSW value is important while the number of interation is not. - One can project the electronic wavefunction on arbitrary positions. It is very helpful to find the integrated charge for given spheres without postprocessing of charge density files (CHG or CHGCAR) - To do that, you should add following lines in POSCAR file, with the use of proper RWIGS, NPAR, LORBIT tags in the INCAR file. Empty (spheres) 1 0.0000000000000000 0.000000000000000 0.0000000000000000 - Drawing local density of states (LDOS) (update: 2013-10-27)
- Local density of states is the real space distribution of the wavefunction square, i.e. charge density for given energy range. Following below steps to obtain the LDOS - Obtain wavefunction -> obtain all partial charge density, rho(n,k,x,y,z) -> You can obtain energy eigenvalue from EIGENVAL -> Now you know rho(E(n,k),x,y,z) -> Now average over rho for y,z. to obtain rho(E,x) -> Now draw rho(E,x) with smearing. - Note that you need very many k-point. I know and you know that it is the time consuming job. - Frequently, one should check the difference between LOCPOT files - In this case, it is a good idea to fix the number of charge-density or local-potential grids. - Related tags are NGX, NGY, NGZ. Maybe sometimes with NGXF, NGYF, NGZF - Phonopy code is the very useful tool for phonon calculations. - Good texts to read: - IR-active mode calculations
- You can calculate the IR-active mode intensities using the following code. - Dielectric function for metallic system
- Dielectric functions are composed of interband-transition-term (bound electron term) plus intraband-transition-term (free electron term: Drude term) - Although you can obtain inerband-transition term from VASP, you should calculate intraband-transition-term by yourself especially for metallic system. - It means that you need plasma frequency from the band structure. There are two possible method to calculate the plasma freqeuncy using vasp code. - Dielectric constant for dielectrics with clamped-ion (update: 2016-01-28)
- Using the LEPSILON-tag, one can calculate the ion-clamped dielectric constant as written in VASP manual [LINK]. - Optimizing structure - Calculating wavefunction - Turning on LEPSILON-tag (LEPSILON = TRUE) and commenting out NPAR-tag or no NPAR tag in the INCAR-file. - Dielectric constant for dielectrics with relaxed-ion (update: 2016-01-28)
- Using the LEPSILON-tag & IBRION-tag, one can calculate the ion-clamped dielectric constant as written in VASP manual [LINK]. - Optimizing structure. Note that ionic forces should be very small because of the force constant calculations. - Calculating wavefunction - LEPSILON = TRUE / no NPAR tag / IBRION = 5,6,7,8 (one of them) - Prediction of plasma frequency (update: 2013-10-27)
- The plasma frequency I mention here is not the 'Screened plasma frequency'. You shoud distinguish between just plasma frequency and screened plasma frequency. The screened plasma frequency is the frequency where the real part of the dielctric function (inter-band term + intra-band term) becomes zero. And the screened plasma frequency is smaller than the plasma frequency. - Using BoltzTraP code, you can obtain electrical conductivity over electon scattering time (SIGMA/TAU) within constant relaxation time approximation. Then, SIGMA/TAU can be transformed to plasma frequency. You shoud be careful when calculating plasma frequency. As the plasma frequency is obtaiend from the intra-band transition, I mean the group velocity of the electron, the k-point sampling is very important. You need sufficiently many k-points to obtain correct band dispersion. It is recommended to test the energy convergence with increasing number of k-points. - You may also obtain the plasma frequency within vasp code. With the input tag of "LOPTICS=TRUE" you can find the "plasma frequency" in the OUTCAR file. - Electrical conductivity using VASP code (update: 2013-11-13, 2015-04-17)
- Within semiclassical theory, one may obtain the electrical conductivity from the band structure by calculating the intraband structure. However, when you use the semiclassical theory, you should get the relaxation time from outside or you should calculate the relaxation time. If you already know the electron scattering time or relaxation time, you can get the electrical conductivity using VASP. - There are two ways to get the electrical conductivity. One is using "BoltzTraP" code. And the other thing is to read the OUTCAR file from VASP code. Using LOPTICS = TRUE tag and RTIME = x.xx (femtotsec) in INCAR with sufficiently many kpts, you can finally found the reliable value of electrical conductivity tensor in OUTCAR from VASP. Unit of electrical conductivity per relaxation time is 10^22 [S/m/sec]. For example, if the tensor is [ [0.020 0.000 0.000] [0.000, 0.020, 0.000] [0.000, 0.000, 0.020]], then Sigma/Tau = 0.020 x 10^22 S/m/sec = 2 x 10^20 S/m/sec For recent VASP version of 5.3.3, you can set the relaxation time ( RTIME in unit of femtosec ) and obtain the electrical conductivity tensor in unit of mega S/m written in OUTCAR. // 2015-04-17 And the Energy dependent conductivity is also written in the "vasprun.xml" file. // 2015-04-17 - You may reduce the sentence length of MAGMOM tag as shown below, MAGMOM = 1 1 1 1 1 ==> MAGMOM = 5*1 MAGMOM = 1 1 -3 -3 4 ==> MAGMOM = 2*1 2*-3 4 It really saves your time to make INCAR file for spin poloarized calculation of large supercell - Many body perturbation theory (MBPT) : GW and BSE in VASP (I tested MBPT calculations with vasp version of 5.2.xx)
- MBPT calculations are very time consuming and very memory consuming. - Calculation flow: (1) DFT(Exc=LDA or GGA) calculations (or DFT+U, Hybrid-DFT) for wavefunctions (WAVECAR) and their derivatives (WAVEDER). You need sufficient many unoccupied bands. (2) GW calculations (ALGO=GW0, scGW0, scGW) (3) BSE calcualtoin (ALGO=BSE) : although you can not find the BSE tag in your manual, it works. And check the "vasprun.xml" which contains frequency dependent dielectric functions. - Related tags: LOPTICS = T for derivative of wavefunction (see manual) NBANDS number of total bands including occupied and unoccupied, convergence test is needed (see manual) ENCUT convergence test is needed (see manual) CSHIFT smearing value for optical response results (freq. dependent dielectric function) (see manual) ALGO (=GW0, scGW0, scGW, BSE) one shot G0W0 tag (GW0), W fixed GW tag (scGW0), fully self consistent GW (scGW), BSE calculation tag (see manual) NOMEGA (see manual) NBANDSO number of occupied bands you considered for BSE calculations NBANDSV number of unoccipied (virtual) bands you considered for BSE calculations - Charge analysis (update: 2015-04-20)
- Bader Charge The concept of Bader Charge is based on the space separation based on the minimum charge density positions. - CHGCAR or CHG file You can directly analyze the atomic charge density using CHGCAR file. Or as mentioned in above (#6), you can get the atomic charge using DOSCAR or OUTCAR. - Unit conversion for VASP phonon frequency from density functional perturbation theory (DFPT) results, Unit of phonopy (update: 2015-06-24)
- [VASP] When one run DFPT calculations or BSE calculations, one obtain the vibration frequency in OUTCAR or vasprun.xml files.VASP. freq (THz) : ang. moment w (2pi THz) : 1/wavelength (1/cm) : energy (meV) = 1 : 2 pi : 10^10/c(speed of light: c) : 10^15xh(planck const: h) = 1 : 6.238186 : 33.35641 : 4.135669 - [phonopy] the default unit for frequency is THz. Therefore you can obtain phonon energy by producting planckconstant (h) 1 THz = 4.135669 meV - example of PbTe phonon (VASP + phonopy, cubic lattice, lattice parameter of a=6.567843. - First prin file:////Users/wangvei/Library/Group%20Containers/UBF8T346G9.Office/TemporaryItems/msohtmlclip/clip_image001.pngciples calculations of Lattice thermal conductivity
- Thermal conductivity (kappa_tot) is a sum of electrical term and lattice term. - Electrical term (kappa_elec) = (Lorenz Number) x ( Elec. conductivitiy ) x (Temp.), where you can obtain Lorenz number and Elec. conductivity from the BoltzTrrap (electron boltzmann transport equation) - - Lattice therm (kappa_latt) = (kappa_phonon) can be obtained by performing phono3py calculations, which calculates the transition rate of phonon-phonon scattering. (phono3py LINK) Using phono3py code with VASP (2015-12-21)
(0) prepare the optimized supercell with the VASP parameter to be used in force calculations. The force should be less than 1.0e-08 eV/angstrom. Of course, you should optimized the structure considering the k-point mesh you will use in (1). WAVECAR file is needed to reduce computational resources for step (2).(1) generate POSCAR files with pair-distance configurations. There might be one to ten thousands POSCAR files. phono3py -d --dim="4 3 3" -c POSCAR-unitcell --cutoff_pair_distance="4.4" supercell size, unitcell file, cutoff pair-distance in angstrom(2) calculate force for pair-distance configurations using VASP. Make directories for configurations of POSCAR-00365.Copy INCAR, KPOINTS, POTCAR POSCAR, and WAVECAR. Using WAVECAR preapared from (0) will be useful to reduce the computational cost. In my case, I used the gamma point only but with the large supercell containing 200-400 atoms. Energy cutoff is also known to be sensitive to the lattice thermal conductivity tensor. For example, when the monoclinic unitcell volume is optimized within 400 eV and internal coordinate is optimized within 300 eV. then there could be non-zero off-digonal term in lattice thermal conductivity tensor.(3) collect vasprun.xml file phono3py --cf3 disp-{00001..00999}/vasprun.xml (or) phono3py --cf3 disp-*/vasprun.xml disp_fc3.yaml, POSCAR-unitcell files should be in same directory.(4) create fc2.hdf and fc3.hdf phono3py --dim="4 3 3" -c POSCAR-unitcell "This step is not madatory, but you can avoid calculating fc2 and fc3 at every run tim." from "http://phonopy.sourceforge.net/phono3py/workflow.html" The file size of fc3.hdf is very large. It is about 1 to 10 GB.(5) calculate thermal conductivity (5.1) split thermal 3phonon process configurations phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm ---wgp q-mesh grid, write grid point for split calculationThe computational time for 3phonon-transition rates are very time consuming when q-mesh grid is fine. Thermal conductivity, especially for low temperature region, q-mesh grid affects the lattice thermal conductivity. Therefore it is recommended to split the calculation set and calculate the transition rate separately using many cpu cores. After run the commands in (5.1), you will get "grid_address-m111111" and "ir_grid_points.yaml"(5.2) calculate 3phonon transition rate for each q grid point phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm --write_gamma --gp="3" Check q grid point number (gp) from ir_grid_points.yaml and write correct gp. The calculation time is dependent on q-mesh size. Each gp-point calculation time is varying from 10 minutes to 10 days depending supercell size and q-mesh size. Note that there are many gp if q-mesh size is fine.(5.3) sum up thermal conductivity phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm --read_gamma > z_TC_def.txt Now you will get the lattice thermal conductivitiy file, "z_TC_def.txt" Type "tail -n 150 z_TC_def.txt" and find the lattice conductivity tensor as a function of temperature.(*) Useful tag: --isotope, --mv, --bmfp for isotope effect, mass variance effect, boundary scattering effect. (**) Auxilary tools kaccum and gaccum commands for accumulated lattice thermal conductivity and 3ph transition rate, respectively. - Band structure using VASP+Wannier90 (2016-04-05)
We can draw band structure by calculatinng Wannier Interpolation using Wannier90. You need wannier90 compiled VASP. Although I don't know well about wannier90, I want to share my knowledge to prevent other's failures Simple method (1-1) Preconverge wavefunction (1-2) Turn on LWANNIER90=T tag and run wan90 compiled VASP. Then Maximally Localized Wannier Functions (MLWFs) are generated with wannier90.win, wannier90.mmn, wannier90.eig, wannier90.out files. Note that, in my case, the vasp run is failed when I use NPAR tag. If you have same problem, comment out the NPAR tag in the INCAR file. (2017-03-24) (1-3) Re-Run wannier90. wannier90.x wannier90 (1-4) Interpolate wannier functionals: first modify the wannier90.win file. For details, see the tutorial file. wannier90.x wannier90 GW band structure from wannier90 I will explain how to make GW band structure. You considered diamond Si and there are two Si atoms in primitive cell. Then there are 8 electrons. For GW calculations you may need many number of bands, about 96 or larger than it. Please note that NPAR should be tagged out or equal to number of nodes(cores). And you may know that NBANDS is dependent on NPAR tag. First check NBANDS without NPAR tag. Also you don't need too many wannier functions. It is a good idea to reduce the number of bands. Actually the number of bands for wannier90 is very sensitive to the band structure. So I recommend to check the suitable number of bands by just running simple DFT. Here I assume that the optimal number of wannier bands is 4 per Si atom (8 per primitive cell). In shortly, you need to determine "NBANDS" for VASP and "NUM_WAN" for wannier90. Here I will use NBANDS=96 (I have 12 core nodes, so NBANDS should be divived by 12) Here I will use NUM_WAN = 8 (the number of wannier functions are related to the symmetries) (1-1) DFT calculations with default NBANDS (1-2) Optical calculations (LOPTICS=TRUE, NBANDS=96 ) (1-3) GW calculations with wannier90 (ALGO=GW0, LWANNIER90=TRUE, no NPAR related tag) If you do not set wannier90.win file, then wannier90.win will be generated. But it is too heavy to play with it since there are too many BANDs. Make wannier90.win file (or modify wannier90.win file) as follows. (The default wannier90.win can be obtained by running without wannier90.win file. GW calculation is too time consuming. Therefore just run DFT with LWANNIER90=T tag) num_wann = 8 num_bands = 8 exclude_bands : 9-96 (1-4) Run wannier90 After getting wannier90.win, wannier90.eig, wannier90.mmn, rerun the wannier90.x cmd> wannier90.x wannier90 Now you wannierize and obtain MLWFs. (1-5) Intepolate First modify or generate or make suitable wannier90.win or NAME.win. And run by typing "wannier90.x wannier90" or "wannier90.x NAME". Now you can interpolate for DOS, BANDS, or Fermi-Surface. IMPORTANT TIP for wannier90 with VASP (2016-11-02): you should comment out the NPAR tag when the calculation is crashed. When I ran the vasp calculation for Mg2Si primitive cell, my calculation is crashed with following error: internal error in GENERATE_KPOINTS_TRANS: G vector not found. - By removing NPAR tag, I did succeed to calculate the wannier90 with VASP.Phonon thermal conductivity of low dimensional structure (2016/04/25)
Please be careful that if there is vacuum, I mean your structure is 1d or 2d structures with vacuum, you should normalize the thermal conductivity tensors. Also be careful that there will be dipole interactions between adjacent supercell and the force from displacement supercell will be dependent on the lateral size of supercell. For example, I tested 1D carbon (C) atomic chain along z direction and the C-C distance is set to 1.28 angstrom. I considered 1x1x10 supercells containing 10 C atoms with various a=b values, 10, 14, 20, and 25 angstroms. I find that the thermal conductivity value of zz component times a square seems to be converged as a becoming larger from 10 to 25 (k_zz x a^2 = 2836, 2968, 3136, 3215 Wxm/K at 300K, 4560, 5181, 5281, 4685 Wxm/K at 500 K. It implies that the cross section area should be checked carefully because the phono3py code assumes that one's structure is bulk. The reason why the k_zz is not exactly same is that there is dipole when single C atom is displaced and it makes large cross sectional interactions between adjacent supercells. - I think that if we consider the neutral atomic chain which has negligible dipole interaction between supercell, the k_zz will be converged more fastly.How to compute defect formation energy of charged defect D with q state.
Eform[Dq] = TE[Dq] - E[with out defect, supercell] - atomic chemical potential + q (electron chemical potential) = TE[Dq] - E[with out defect, supercell] - atomic chemical potential + q (Fermi level + Evbm). Here, the Evbm position should be used for charged defect supercell. However, due to the charge distribution change, the Evbm or Ecbm of defective system is not same to the Evbm of Ecbm of bulk without defect. So we should find the proper position of vbm by aligning the local potentials. However, due to the small supercell size, the aligning supercell is size dependent and the unwanted charged defect interaction is generated. So, one should remove these two effects. For local potential ailgnment, the local potential change is estimaed with a functoin of a distance from defect. Because in the relaxed supercell the local potential average is not easy, for the easy way, I recommend to see the atomic local potential around the atomic site with a distance of 1 angstrom. This atomic site local potential is already computed and written in the VASP code. It is relatively convenient to use. VASP-MD Simulations http://blog.163.com/xiaowei_090513/blog/static/1177183592010099379680/ This file determines the kind of job which VASP will perform; single point energy calculation (SPE), geometry optimisation (GO - coarse/fine), molecular dynamics (MD - nve/nvt), spin polarised calculation (mag). Examples can be found in /home/cs/model/vasp_util. Example; INCAR.spe $system = single point energy calc NELMIN = 4 minimum number of electronic SCF cycles EDIFF = 1E-6 stooping criterion for electronic convergence NSW = 0 number of ionic shifts ISMEAR = 0 treatment of partial occupancies of electronic levels Example; INCAR.coarse $system = coarse geom optimisation NELMIN = 4 EDIFF = 1E-2 EDIFFG = -1E-2 stopping criterion for forces Fmax < 0.01 eV/A IBRION = 2 minimisation method, good away from minimum ISIF = 3 optimise coords and cell pars LREAL =.TRUE. do calc in real space - quicker ISTART = 0 start with a random wavefunction NSW = 20 maximum of 20 ionic shifts ISMEAR = 0 LCHARG =.FALSE. don't write CHG and CHGCAR files Example; INCAR.fine $system = geom optimisation NELMIN = 4 EDIFF = 1E-6 EDIFFG = -1E-4 PREC = high increase energy cut-off by 25% IBRION = 1 minimisation method, good close to minimum ISIF = 3 NSW = 50 ISMEAR = 0 LCHARG=.FALSE. Example; INCAR.mag $system = collinear mag structure calc IBRION = 1 ISIF = 3 NPAR = 1 forces mag structure to be written in output file EDIFF = 1E-6 EDIFFG = -1E-3 PREC = high RWIGS = 1.376 0.900 1.233 1.302 radii for spherical integration of spin density, 1 per atom ISPIN = 2 do spin polarised calc MAGMOM = 24*0 5 -5 -5 5 initial mag moments for 28 atomsNSW = 20 Example; INCAR.nve $system = molecular dynamics ALGO = V MAXMIX = 40 IBRION = 0 do molecular dynamics NSW = 6000 number of time steps NBLOCK = 1 store structure every time step POTIM = 3.0 time step 3fs TEBEG = 673 target temperature ISYM = 0 turn off symmetry SMASS = -3 NVE ensemble LREAL =.TRUE. LCHARG =.FALSE. NELMIN = 4 PREC = LOW reduce energy cut-off by 25% for MD ISTART = 0 ISMEAR = 0; SIGMA=0.1 Example; INCAR.nvt $system = molecular dynamics ALGO = V MAXMIX = 40 IBRION = 0 NSW = 6000 NBLOCK = 1 POTIM = 3.0 TEBEG = 673 ISYM = 0 SMASS = 2 NVT ensemble, value determines frequency of coupling to heat bath LREAL =.TRUE. LCHARG =.FALSE. NELMIN = 4 PREC = LOW ISTART = 0 ISMEAR = 0; SIGMA=0.1 Example; INCAR.scale $system = molecular dynamics quench ALGO = V MAXMIX = 40 IBRION = 0 NSW = 50 NBLOCK = 5 rescale temperature every 5 steps POTIM = 3.0 TEBEG = 683 in |
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=50 NBLOCK*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 -3 NVE ensemble;-1 用来做模拟退火。大于0 NVT 系综。 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 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. 正则系综（canonical ensemble），全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数（N）、体积（V）、温度（T）。正则系综是蒙特卡罗方法模拟处理 的典型代表。假定N个粒子处在体积为V的盒子内，将其埋入温度恒为T的热浴中。此时，总能量（E）和系统压强（P）可能在某一平均值附近起伏变化。平衡体 系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能F（N,V,T）。 2. 微正则系综（micro-canonical ensemble），简写为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. 巨正则系综（grand canonical ensemble），简写为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）dp dq表示在t时刻出现在（p，q）点附近相体积元dp dq内的相点数在全部相点数中所占的比值，即表示任一系统在t时刻其运动状态处于（p，q）附近的相体积元dp dq内的概率。显然 ，概率密度函数ρ（p，q，t）满足归一化条件∫ρ（p,q,t）dpdq=1。 统计物理学的认为系统的任意宏观量I（t）是相应微观量L（p，q）在一定宏观条件下对系统一切可能的微观运动状态的统计平均值，即 I（t）＝∫L（p，q）ρ（p，q，t）dp dq。由此可见，经典统计物理的基本课题是确定各种条件下系综的概率密度函数ρ（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=50 NBLOCK*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 -3 NVE ensemble;-1 用来做模拟退火。大于0 NVT 系综。 SMASS=1,2,3 是没有区别的。都是NVT ensemble。 SMASS只要是大于0就是NVT系综。 学习了，但是要指出你的一点小错误。 NBLOCK= 1 指的是多少个离子步（时间步长只是对于MD而言，对于其它的计算，NBLOCK也是起控制作用的）写一次XDATCAR KBLOCK=50 NBLOCK*KBLOCK 个离子步写一次PCDAT. CONTCAR是每个离子步之后都会写出来的，但是会用新的把老的覆盖 CHG是在每10个离子步写一次，不会覆盖 CHGCAR是在任务正常结束之后才写的。 如何用VASP计算铁磁、反铁磁和顺磁 顺磁，意味进行non-spin polarized的计算，也就是ISPIN=1。 铁磁，意味进行spin-polarized的计算，ISPIN=2，而且每个磁性原子的初始磁矩设置为一样的值，也就是磁性原子的MAGMOM设置为一 样的值。对非磁性原子也可以设置成一样的非零值（与磁性原子的一样）或零，最后收敛的结果，非磁性原子的local磁矩很小，快接近0，很小的情况，很可 能意味着真的是非磁性原子也会被极化而出现很小的local磁矩。 反铁磁，也意味着要进行spin-polarized的计算，ISPIN=2，这是需采用反铁磁的磁胞来进行计算，意味着此时计算所采用的晶胞不再是铁磁 计算时的最小原胞。比如对铁晶体的铁磁状态，你可以采用bcc的原胞来计算，但是在进行反铁磁的Fe计算，这是你需要采用sc的结构来计算，计算的晶胞中 包括两个原子，你要设置一个原子的MAGMOM为正的，另一个原子的MAGMOM设置为负，但是它们的绝对值一样。因此在进行反铁磁的计算时，应该确定好 反铁磁的磁胞，以及磁序，要判断哪种磁序和磁胞是最可能的反铁磁状态，那只能是先做好各种可能的排列组合，然后分别计算这些可能组合的情况，最后比较它们 的总能，总能最低的就是可能的磁序。同样也可以与它们同铁磁或顺磁的进行比较。了解到该材料究竟是铁磁的、还是顺磁或反铁磁的。 亚铁磁，也意味要进行spin-polarized的计算，ISPIN=2，与反铁磁的计算类似，不同的是原子正负磁矩的绝对值不是样大。非共线的磁性， 那需采用专门的non-collinear的来进行计算，除了要设置ISPIN，MAGMOM的设置还需要指定每个原子在x,y,z方向上的大小。这种情 况会复杂一些。 举个例子来说，对于Mn-Cu(001)c(2x2)这种体系，原胞里面有2个Mn原子，那么你直接让两个Mn原子的MAGMOM的绝对值一样，符号相反就可以了，再加上ISPIN=2。这样就可以实现进行反铁磁的计算了。 | 
发表于 2021-12-27 10:24:34
发表于 2021-12-28 08:15:33
