博文
Phonopy 热力学性质及ZPE的计算
|||
关注：
1） 热力学性质计算的原理；
2） 简谐计算与准简谐计算的差别；即Phonopy 与Phonopy-qha得到的热力学性质的差别在哪
3） Phonnopy 计算热力学性质的具体步骤。 把关键参数写入mesh.conf 文件，执行phonopy -t mesh.conf即可
4) Phonopy 扩胞时，没找到正确的对称性应如何处理，如不处理对计算结果有什么影响？
phonopy -d --dim="3 3 1" --tolerance=0.001
check了一下，发现 Phonopy 扩胞时，没找到正确的对称性得到的SPOSCAR与找到正确对称性时扩胞后的SPOSCAR完全一样，请问，是否可据此为没找到正确对称性对计算结果应该没有什么影响？
题记：A friend in need is a friend indeed. Thanks must be given to Prasad and Huayung
Prasad, my real friend, he is so glad to share his knowledge and experience with others. He should become a good teacher. However, not everyone like him, one guy in my group is so caring her knowledge and experience in case that will be learned by others. Which kind of person will we like? Which kind of person should we be? It is underthinking. And I have my answer as you guys see this blog.
So, though Prasad left here and was far away from us. I still looked for help from him and got this hints how to calculate thermodynamic propeties.
Phonnopy 计算热力学性质的具体步骤：
About phonopy: In order to get ZPE we have to run thermal properties. The command to run thermal properties in phonopy is:
phonopy -t mesh.conf
this command exactly similar to the command to obtain the phonon DOS [ phonopy -p mesh.conf], the only difference is the -t option.
All the four files ( that we use to obtain the phonon DOS ) are necessary to obtain thermal properties.
1. mesh.conf,
2. vasprun.xml,
3. POSCAR,
4. FORCE_CONSTANT
once we run the command: phonopy -t mesh.conf
it generates a file called thermal_properties.yaml, this is a simple txt file.
Please open this file, almost at the top of the file you see zero_point_energy:
the units are kJ/mol/unit cell [ the unit cell is POSCAR, not the SUPERCELL called SPOSCAR]
When we ran the command: phonopy -t mesh.conf
along with thermal_properties.yaml the program also writes the text file on the screen. This text is same as the text written in thermal_properties.yaml. So don't worry about this. Use the thermal_properties.yaml file to obtain the ZPE.
The ZPE values are in kJ/mol/unit cell. Prof. Neil may ask to convert into eV. So please convert to eV.
This mesh.conf from Huayun
ATOM_NAME = H Sc
DIM = 3 3 1
FORCE_CONSTANTS = READ
MP = 8 8 8
GAMMA_CENTER = .TRUE.
TPROP = .TRUE.
英文手册解读：
Thermal properties, free energy, heat capacity, and entropy, are calculated from their statistical thermodynamic expressions (see Thermodynamic properties).
Thermal properties are calculated from phonon frequencies on a sampling mesh in the reciprocal space. Therefore These tags are used with MP tag and their convergence with respect to the sampling mesh has to be checked. Usually this calculation is not computationally demanding, so the convergence is easily achieved with increasing the density of the sampling mesh.
-p option can be used together to plot the thermal propreties.
Phonon frequencies have to be calculated in THz. Therefore unit conversion factor to THz may be specified using --factor option.
The calculated values are written into thermal_properties.yaml.
The unit systems of free energy, heat capacity, and entropy are kJ/mol, J/K/mol, and J/K/mol, respectively, where 1 mol means 
your input unit cell (not formula unit), i.e. you have to divide the value by number of formula unit in your unit cell by yourself. For example, in MgO (conventional) unit cell, if you want to compare with experimental results in kJ/mol, you have to divide the phonopy output by four.
TMIN, TMAX, and TSTEP tags are used to specify the temperature range to be calculated. The default values of them are 0, 1000, and 10, respectively.
Experimental
Mean square displacements projected to Cartesian axes as a function of temperature are calculated from the number of phonon excitations. The usages of TMAX, TMIN, TSTEP tags are same as those in thermal properties tags. The result is writen into thermal_displacements.yaml. See the detail of the method,Thermal displacement.
Eigenvectors are projected along the direction specified by this tag. Projection direction is specified in reduced coordinates, i.e., with respect to a, b, c axes.
Thermodynamic properties
![E = .sum_{.mathbf{q}.nu}.hbar.omega(.mathbf{q}.nu).left[.frac{1}{2} +
.frac{1}{.exp(.hbar.omega(.mathbf{q}.nu)/k_.mathrm{B} T)-1}.right]](http://phonopy.sourceforge.net/_images/math/0a5af4ec198d887291070bcb21b5e16d0eef91bf.png)
![C_V &= .left(.frac{.partial E}{.partial T} .right)_V ..
&= .sum_{.mathbf{q}.nu} k_.mathrm{B}
.left(.frac{.hbar.omega(.mathbf{q}.nu)}{k_.mathrm{B} T} .right)^2
.frac{.exp(.hbar.omega(.mathbf{q}.nu)/k_.mathrm{B}
T)}{[.exp(.hbar.omega(.mathbf{q}.nu)/k_.mathrm{B} T)-1]^2}](http://phonopy.sourceforge.net/_images/math/ff44399310b797d866c5fe8622260e59350bbe44.png)
![F &= -k_.mathrm{B} T .ln Z ..
&= .varphi + .frac{1}{2} .sum_{.mathbf{q}.nu}
.hbar.omega(.mathbf{q}.nu) + k_.mathrm{B} T .sum_{.mathbf{q}.nu} .ln
.bigl[1 -.exp(-.hbar.omega(.mathbf{q}.nu)/k_.mathrm{B} T) .bigr]](http://phonopy.sourceforge.net/_images/math/6c0581b8bb4f1c33e2873fa14c575ce3f87c0416.png)
![S &= -.frac{.partial F}{.partial T} ..
&= .frac{1}{2T}.sum_{.mathbf{q}.nu}.hbar.omega(.mathbf{q}.nu).coth(.hbar.omega(.mathbf{q}.nu)/2k_.mathrm{B}T)-k_.mathrm{B} .sum_{.mathbf{q}.nu}.ln.left[2.sinh(.hbar.omega(.mathbf{q}.nu)/2k_.mathrm{B}T).right]](http://phonopy.sourceforge.net/_images/math/6cf55cf9080958655c884d0a005a3057868aed50.png)
MP numbers give uniform meshes in each axis. As the default behavior, the center of mesh is determined by the Monkhorst-Pack scheme, i.e., for odd number, a point comes to the center, and for even number, the center is shifted half in the distance between neighboring mesh points.
Examples of an even mesh with 
center in two ways,
MP_SHIFT gives the shifts in direction along the corresponding reciprocal axes (
, 
, 
). 0 or 1/2 (0.5) can be used as these values. 1/2 means the half mesh shift with respect to neighboring grid points in each direction.
Instead of employing the Monkhorst-Pack scheme for the mesh sampling, center mesh is used. The default value is .FALSE..
Total and partial density of states are drawn with some parameters. The example makes DOS be calculated from frequency=0 to 40 with 0.1 pitch.
By setting this tag, EIGENVECTORS = .TRUE. is automatically set. PDOS tag controls how elements of eigenvectors are added. Each value gives the atom index in primitive cell. , separates the atom sets. Therefore in the example, atom 1 and 2 are summarized as one curve and atom 3, 4, 5, and, 6 are summarized as the other curve.
This tag specifies the deviation of a smearing function. The unit is same as that of final result of DOS, i.e., for VASP without --factor option, it is THz. The default value is the value given by the difference of maximum and minimum frequencies divided by 100.
By setting .TRUE., DOS at lower phonon frequencies are fit to a Debye model. By default, the DOS from 0 to 1/4 of the maximum phonon frequencies are used for the fitting. The function used to the fitting is 
where 
is the parameter and the Debye frequency is 
where 
is the number of atoms in unit cell. Users have to unserstand that this is not a unique way to determine Debye frequency. Debye frequency is dependent on how to parameterize it.
DEBYE_MODEL = .TRUE.
对称性问题：
phonopy -d --dim="3 3 1"
_
_ __ | |__ ___ _ __ ___ _ __ _ _
| '_ | '_ / _ | '_ / _ | '_ | | | |
| |_) | | | | (_) | | | | (_) || |_) | |_| |
| .__/|_| |_|___/|_| |_|___(_) .__/ __, |
|_| |_| |___/
1.6.1
Creating displacements
Supercell: [3 3 1]
Spacegroup: P4mm (99) 对称性实际应为p4/mmm (123)
disp.yaml and supercells have been created.
_
___ _ __ __| |
/ _ '_ / _` |
| __/ | | | (_| |
___|_| |_|__,_|
http://phonopy.sourceforge.net/command-options.html?highlight=tolerance
The specified value is used as allowed tolerance to find symmetry of crystal structure. The default value is 1e-5.
Using this option, various crystal symmetry information is just printed out and phonopy stops without going to phonon analysis.
This tag can be used together with the --cell, --wien2k, or --primitive_axis option
资料搜索：
1. phonopy-qha.使用说明
This example is to explain how to use phonopy-qha.
Following the instruction, thermal properties at constant pressure (P=0) for silicon may be calculated.
The shell script can be used only for the crystal structure whose atoms are pinned at the reduced positions by the crystal symmetry.
For general crystals, the unit cells have to be relaxed before starting phonon calculations.
Create VASP input files. POTCAR is not included in this example.
% phonopy-qha e-v.dat QHA-*/thermal_properties.yaml -p --sparse=50
calc.sh内容
#$ -S /bin/zsh #$ -cwd #$ -N Sivol-num #$ -pe mpi* 4 #$ -e err.log #$ -o std.log cp $TMPDIR/machines machines mpirun vasp5212mpi >& ../progress-num |
scripqha.sh内容（不用修改）
phonopy-qha e-v.dat thermal_properties.yaml-{-{5..1},{0..5}} --sparse=50
# Vinet EOS
# T E_0 B_0 B’_0 V_0
0.000000 -14.796263 75.231724 4.758283 66.697923
2.000000 -14.796263 75.231723 4.758283 66.697923
4.000000 -14.796263 75.231718 4.758284 66.697923
6.000000 -14.796263 75.231695 4.758286 66.697924
8.000000 -14.796263 75.231634 4.758294 66.697928
10.000000 -14.796264 75.231510 4.758308 66.697934
9.8.1 -t, --tmax, --tmin, --tstep
These correspond to TPROP, TMAX, TMIN, and TSTEP tags, respectively (Thermal properties related tags).
LLL:
thermal_properties.yaml；一般都在force文件夹里。
A::::好的，我找找看....
LLL:::
原子数目;单位是KJ/mol
乘以0.010364 ；是电子伏特natom: 8
zero_point_energy: 217.3900395
high_T_entropy: 374.5454014
Using phonopy results of thermal properties, thermal expansion and heat capacity at constant pressure can be calculated under the quasi-harmonic approximation.
phonopy-qha is the script to calculate them. An example of the usage is as follows:
phonopy-qha e-v.dat thermal_properties-{1..10}.yaml
1st argument is the filename of volume-energy data (in the above expample, e-v.dat).
The volume and energy of the cell (default units are in 3 and eV, respectively). An example of the volume-energy file is:
# cell volume energy of cell other than phonon
156.7387309525 -104.5290025375
154.4138492700 -104.6868148175
.....
Lines starting with # are ignored.
The other arguments are the filenames of thermal_properties.yaml calculated
at the respective volumes given in the 1st argument.
The thermal_properties.yaml at volume points have to be calculated with the same temperature ranges and same temperature steps.
thermal_properties.yaml can be calculated by following Thermal properties related tags, where the physical unit of the Helmholtz free energy is kJ/mol as the default, i.e., no need to convert the physical unit in usual cases.
The example for Aluminum is found in the example directory.
If the condition under puressure is expected, PV terms may be included in the energies.
https://blog.sciencenet.cn/blog-567091-726895.html
上一篇：PWSCF计算小技巧
下一篇：PWSCF计算声子谱-gamma点的测试很重要
