使用Multiwfn+VMD绘制分子的表面静电势图(ESP)

May 30, 2025 in 从零开始的科研之旅, first stage

Multiwfn版本：3.8(dev)，Last update: 2025-May-5
VMD版本：1.9.3

VMD版本请尽可能使用1.9.3而不要使用bug较多的1.9.4

目标分子的结构优化

本次的目标分子为二茂铁分子

首先建立二茂铁分子的模型，并且手动设置全重叠的D5d点群，这样可以极大的减少优化用时，关键词如下：

1	#p tpssh/def2tzvp opt freq

体系是过渡金属配合物，因此选择tpssh泛函，由于初猜结构比较合适，算得还是挺快的

                                   -- CHARLES MCCARRY IN "THE GREAT SOUTHWEST"
Job cpu time:       0 days  1 hours 57 minutes 15.2 seconds.
Elapsed time:       0 days  0 hours  5 minutes 58.7 seconds.
File lengths (MBytes):  RWF=    147 Int=      0 D2E=      0 Chk=     16 Scr=      1
Normal termination of Gaussian 16 at Wed May 28 21:43:45 2025.

如下图，优化完成，没有虚频

绘图前的准备工作

如Sob老师所说，绘图之前用户要做的事包括：

把examples.txt和.bat文件拷到Multiwfn可执行文件所在目录
把.bat文件里的VMD路径改成实际路径
把examples.vmd拷到VMD路径下
在vmd.rc里加入上述proc语句

画图

有两种风格的分子表面静电势图(ESP)，分别是表面顶点着色方式和电子密度等值面着色方式的ESP图，我个人更喜欢后者，以后者为例画图吧。

首先运行ESPiso.bat产生电子密度和静电势的cube文件：

PS C:\Multiwfn_3.8_dev_bin_Win64> ls

    Directory: C:\Multiwfn_3.8_dev_bin_Win64

Mode                 LastWriteTime         Length Name
----                 -------------         ------ ----
d----           2025/5/15    12:20                examples
-a---           2025/5/23     9:20        8114156 1.fch
-a---           2024/9/21    17:46           3057 使用Multiwfn发表文章必须在正文里进行引用（包括代算）.txt
-a---           2025/5/28    17:51            105 ESPext.bat
-a---           2018/10/7     5:24             32 ESPext.txt
-a---           2025/5/28    17:51            460 ESPiso.bat
-a---          2021/11/18    10:43             26 ESPiso.txt
-a---           2025/5/28    17:51            344 ESPpt.bat
-a---           2018/9/30    13:53             29 ESPpt.txt
-a---            2025/2/9     9:28         200205 How to cite Multiwfn.pdf
-a---            2023/3/2     2:26        2047000 libiomp5md.dll
-a---           2024/8/27    20:00           1254 LICENSE.txt
-a---           2025/4/13     9:17         416416 Multiwfn quick start.pdf
-a---            2025/5/5     2:50       35666432 Multiwfn.exe
-a---           2025/5/16     1:26          20490 settings.ini

PS C:\Multiwfn_3.8_dev_bin_Win64> .\ESPiso.bat

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 1.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:12 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Please wait...
 Loading various information of the wavefunction
 The highest angular moment basis functions is F
 Loading basis set definition...
 Loading orbitals...
 Converting basis function information to GTF information...
 Back converting basis function information from Cartesian to spherical type...
 Generating density matrix based on SCF orbitals...
 Generating overlap matrix...

 Total/Alpha/Beta electrons:     98.0000     49.0000     49.0000
 Net charge:     0.00000      Expected multiplicity:    1
 Atoms:     23,  Basis functions:    427,  GTFs:    708
 Total energy:   -1651.955643142885 Hartree,   Virial ratio:  2.00246837
 This is a restricted single-determinant wavefunction
 Orbitals from 1 to    49 are occupied
 Title line of this file: Generated by autoGau.

 Loaded 1.fch successfully!

 Formula: H12 C10 Fe1      Total atoms:      23
 Molecule weight:       188.04779 Da
 Point group: C1

 "q": Exit program gracefully          "r": Load a new file
                    ************ Main function menu ************
 0 Show molecular structure and view orbitals
 1 Output all properties at a point       2 Topology analysis
 3 Output and plot specific property in a line
 4 Output and plot specific property in a plane
 5 Output and plot specific property within a spatial region (calc. grid data)
 6 Check & modify wavefunction
 7 Population analysis and calculation of atomic charges
 8 Orbital composition analysis           9 Bond order analysis
 10 Plot total DOS, PDOS, OPDOS, local DOS, COHP and photoelectron spectrum
 11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA/NMR spectrum
 12 Quantitative analysis of molecular surface
 13 Process grid data (No grid data is presented currently)
 14 Adaptive natural density partitioning (AdNDP) analysis
 15 Fuzzy atomic space analysis
 16 Charge decomposition analysis (CDA) and plot orbital interaction diagram
 17 Basin analysis                       18 Electron excitation analysis
 19 Orbital localization analysis        20 Visual study of weak interaction
 21 Energy decomposition analysis        22 Conceptual DFT (CDFT) analysis
 23 ETS-NOCV analysis                    24 (Hyper)polarizability analysis
 25 Electron delocalization and aromaticity analyses
 26 Structure and geometry related analyses
 100 Other functions (Part 1)            200 Other functions (Part 2)
 300 Other functions (Part 3)
 -10 Return to main menu
 -2 Obtain deformation property
 -1 Obtain promolecule property
 0 Set custom operation
             ----------- Available real space functions -----------
 1 Electron density (rho)     2 Gradient norm of rho     3 Laplacian of rho
 4 Value of orbital wavefunction         44 Orbital probability density
 5 Electron spin density
 6 Hamiltonian kinetic energy density K(r)
 7 Lagrangian kinetic energy density G(r)
 8 Electrostatic potential from nuclear charges
 9 Electron localization function (ELF)
 10 Localized orbital locator (LOL)
 11 Local information entropy
 12 Total electrostatic potential (ESP)
 13 Reduced density gradient (RDG)       14 RDG with promolecular approximation
 15 Sign(lambda2)*rho      16 Sign(lambda2)*rho with promolecular approximation
 17 Correlation hole for alpha, ref. point:   0.00000   0.00000   0.00000
 18 Average local ionization energy (ALIE)
 19 Source function, mode: 1, ref. point:   0.00000   0.00000   0.00000
 20 Electron delocal. range func. EDR(r;d)  21 Orbital overlap dist. func. D(r)
 22 Delta-g (promolecular approximation)    23 Delta-g (Hirshfeld partition)
 24 Interaction region indicator (IRI)    25 van der Waals potential (probe=C )
 100 User-defined function (iuserfunc=   -1), see Section 2.7 of manual

 Please select a method to set up grid
 -10 Set extension distance of grid range for mode 1~4, current:  6.000 Bohr
 1 Low quality grid,    covering whole system, about 125000 points in total
 2 Medium quality grid, covering whole system, about 512000 points in total
 3 High quality grid,   covering whole system, about 1728000 points in total
 4 Input the number of points or grid spacing in X,Y,Z, covering whole system
 5 Input original point, grid spacings, and the number of points
 6 Input center coordinate, number of points and extension distance
 7 The same as 6, but input two atoms, the midpoint will be defined as center
 8 Use grid setting of another cube file
 10 Set box of grid data visually using a GUI window
 11 Select a set of atoms, set extension distance around them and grid spacing
 100 Load a set of points from external file
 Coordinate of origin in X,Y,Z is     -12.142516  -10.564851  -10.153107 Bohr
 Coordinate of end point in X,Y,Z is   11.971265   10.285411   10.334542 Bohr
 Grid spacing in X,Y,Z is    0.181307    0.181307    0.181307 Bohr
 Number of points in X,Y,Z is  134  116  114   Total:     1772016
 Note: Virtual orbitals higher than LUMO+9 have been temporarily discarded for saving computational time

 Note: All exponential functions exp(x) with x< -40.000 will be ignored
 Unique GTFs have been constructed. Number of unique GTFs:     708
 Progress: [##################################################]   100.0 %     \
 Calculation of grid data took up wall clock time         4 s
 Note: Previous orbital information has been restored

 Electric dipole moment estimated by integrating electron density
 X component:          0.281807 a.u.      0.716283 Debye
 Y component:         -0.255113 a.u.     -0.648432 Debye
 Z component:          0.140547 a.u.      0.357235 Debye
 Total magnitude:      0.405279 a.u.      1.030117 Debye

 The minimum is  0.60334894E-17 at -12.14252 -10.56485  10.33454 Bohr
 The maximum is  0.95394143E+03 at   0.00503   0.13224   0.00006 Bohr
 Summing up all value and multiply differential element:
   100.725159050150
 Summing up positive value and multiply differential element:
   100.725159050150
 Summing up negative value and multiply differential element:
  0.000000000000000E+000

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
 Exporting cube file, please wait...
 Done! Grid data has been exported to density.cub in current folder
 Hint: If you want to add input file name as prefix of the outputted cube file, you can set "iaddprefix" in settings.ini to 1

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)

 Note: A set of grid data is presented in memory
 "q": Exit program gracefully          "r": Load a new file
                    ************ Main function menu ************
 0 Show molecular structure and view orbitals
 1 Output all properties at a point       2 Topology analysis
 3 Output and plot specific property in a line
 4 Output and plot specific property in a plane
 5 Output and plot specific property within a spatial region (calc. grid data)
 6 Check & modify wavefunction
 7 Population analysis and calculation of atomic charges
 8 Orbital composition analysis           9 Bond order analysis
 10 Plot total DOS, PDOS, OPDOS, local DOS, COHP and photoelectron spectrum
 11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA/NMR spectrum
 12 Quantitative analysis of molecular surface
 13 Process grid data
 14 Adaptive natural density partitioning (AdNDP) analysis
 15 Fuzzy atomic space analysis
 16 Charge decomposition analysis (CDA) and plot orbital interaction diagram
 17 Basin analysis                       18 Electron excitation analysis
 19 Orbital localization analysis        20 Visual study of weak interaction
 21 Energy decomposition analysis        22 Conceptual DFT (CDFT) analysis
 23 ETS-NOCV analysis                    24 (Hyper)polarizability analysis
 25 Electron delocalization and aromaticity analyses
 26 Structure and geometry related analyses
 100 Other functions (Part 1)            200 Other functions (Part 2)
 300 Other functions (Part 3)
 -10 Return to main menu
 -2 Obtain deformation property
 -1 Obtain promolecule property
 0 Set custom operation
             ----------- Available real space functions -----------
 1 Electron density (rho)     2 Gradient norm of rho     3 Laplacian of rho
 4 Value of orbital wavefunction         44 Orbital probability density
 5 Electron spin density
 6 Hamiltonian kinetic energy density K(r)
 7 Lagrangian kinetic energy density G(r)
 8 Electrostatic potential from nuclear charges
 9 Electron localization function (ELF)
 10 Localized orbital locator (LOL)
 11 Local information entropy
 12 Total electrostatic potential (ESP)
 13 Reduced density gradient (RDG)       14 RDG with promolecular approximation
 15 Sign(lambda2)*rho      16 Sign(lambda2)*rho with promolecular approximation
 17 Correlation hole for alpha, ref. point:   0.00000   0.00000   0.00000
 18 Average local ionization energy (ALIE)
 19 Source function, mode: 1, ref. point:   0.00000   0.00000   0.00000
 20 Electron delocal. range func. EDR(r;d)  21 Orbital overlap dist. func. D(r)
 22 Delta-g (promolecular approximation)    23 Delta-g (Hirshfeld partition)
 24 Interaction region indicator (IRI)    25 van der Waals potential (probe=C )
 100 User-defined function (iuserfunc=   -1), see Section 2.7 of manual

 Please select a method to set up grid
 -10 Set extension distance of grid range for mode 1~4, current:  6.000 Bohr
 1 Low quality grid,    covering whole system, about 125000 points in total
 2 Medium quality grid, covering whole system, about 512000 points in total
 3 High quality grid,   covering whole system, about 1728000 points in total
 4 Input the number of points or grid spacing in X,Y,Z, covering whole system
 5 Input original point, grid spacings, and the number of points
 6 Input center coordinate, number of points and extension distance
 7 The same as 6, but input two atoms, the midpoint will be defined as center
 8 Use grid setting of another cube file
 10 Set box of grid data visually using a GUI window
 11 Select a set of atoms, set extension distance around them and grid spacing
 100 Load a set of points from external file
 Coordinate of origin in X,Y,Z is     -12.142516  -10.564851  -10.153107 Bohr
 Coordinate of end point in X,Y,Z is   11.789958   10.321672   10.298280 Bohr
 Grid spacing in X,Y,Z is    0.435136    0.435136    0.435136 Bohr
 Number of points in X,Y,Z is   56   49   48   Total:      131712
 Note: Virtual orbitals higher than LUMO+9 have been temporarily discarded for saving computational time

 Unique GTFs have been constructed. Number of unique GTFs:     708
 Note: ESP will be calculated only for the grids around isosurface of electron density of   0.001000 a.u.
 Detecting the grids for calculating ESP...
 Number of grids to calculate ESP:       10900

 Initializing LIBRETA library (fast version) for ESP evaluation ...
 LIBRETA library has been successfully initialized!

 NOTE: The ESP evaluation code based on LIBRETA library is being used. Please cite Multiwfn original papers (J. Comput. Chem., 33, 580-592 (2012) and J. Chem. Phys., 161, 082503 (2024)) and the paper describing the efficient ESP evaluation algorithm adopted by Multiwfn (Phys. Chem. Chem. Phys., 23, 20323 (2021))
 Progress: [##################################################]   100.0 %     \
 Setting ESP of the grids neighbouring to boundary grids...
 Calculation of grid data took up wall clock time        13 s
 Note: Previous orbital information has been restored

 The minimum is -0.32612104E-01 at  -0.39385  -0.99186   4.64151 Bohr
 The maximum is  0.52214027E-01 at   2.21697  -2.29727  -5.36661 Bohr
 Summing up all value and multiply differential element:
   3.07509251072690
 Summing up positive value and multiply differential element:
   8.67387639395045
 Summing up negative value and multiply differential element:
  -5.59878388322356

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
 Exporting cube file, please wait...
 Done! Grid data has been exported to totesp.cub in current folder
 Hint: If you want to add input file name as prefix of the outputted cube file, you can set "iaddprefix" in settings.ini to 1

                   ---------- Post-processing menu ----------
 -1 Show isosurface graph
 0 Return to main menu
 1 Save graph of isosurface to file in current folder
 2 Export data to a Gaussian-type cube file in current folder
 3 Export data to output.txt in current folder
 4 Set the value of isosurface to be shown, current:   0.05000
 5 Multiply all grid data by a factor
 6 Divide all grid data by a factor
 7 Add a value to all grid data
 8 Substract a value from all grid data
 9 Multiply all grid data by Hirshfeld weights of a fragment (can be used to only make isosurface around interested fragment visible)
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF60135CEB4  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF600F5B65C  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127BDD6  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density1.cub
移动了         1 个文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP1.cub
移动了         1 个文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 2.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density2.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP2.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 3.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density3.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP3.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>Multiwfn 4.fch -ESPrhoiso 0.001  0<ESPiso.txt
 Error: Unable to find the input file you specified in argument

 Multiwfn -- A Multifunctional Wavefunction Analyzer
 Version 3.8(dev), update date: 2025-May-5
 Developer: Tian Lu (Beijing Kein Research Center for Natural Sciences)
 Multiwfn official website: http://sobereva.com/multiwfn
 Multiwfn English forum: http://sobereva.com/wfnbbs
 Multiwfn Chinese forum: http://bbs.keinsci.com/wfn
 ( Number of parallel threads:   4  Current date: 2025-05-30  Time: 20:33:30 )

 Both following papers ***MUST BE CITED IN MAIN TEXT*** if Multiwfn is used:
  Tian Lu, Feiwu Chen, J. Comput. Chem., 33, 580 (2012) DOI: 10.1002/jcc.22885
  Tian Lu, J. Chem. Phys., 161, 082503 (2024) DOI: 10.1063/5.0216272
 See "How to cite Multiwfn.pdf" in Multiwfn binary package for more information

 Now input file path, for example, E:\Planetarian\Yumemi_Hoshino.mwfn
 (.wfn/wfn/wfx/fch/molden/pdb/xyz/mol2/cif/cub... see Section 2.5 of manual)
 Hint: Pressing ENTER button directly can select a file in a GUI window. To reload the past file, inputting "o". Input such as ?miku.fch can open the miku.fch in the same folder as the past file
"5" cannot be found, input again
"1" cannot be found, input again
"3" cannot be found, input again
"2" cannot be found, input again
"0" cannot be found, input again
"5" cannot be found, input again
"12" cannot be found, input again
"1" cannot be found, input again
"2" cannot be found, input again
forrtl: severe (24): end-of-file during read, unit -4, file CONIN$
Image              PC                Routine            Line        Source
Multiwfn.exe       00007FF601370B01  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF60127CBB3  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601CCFB7B  Unknown               Unknown  Unknown
Multiwfn.exe       00007FF601FE53F0  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF9259B7374  Unknown               Unknown  Unknown
ntdll.dll          00007FF92667CC91  Unknown               Unknown  Unknown

C:\Multiwfn_3.8_dev_bin_Win64>move /Y density.cub density4.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y totesp.cub ESP4.cub
系统找不到指定的文件。

C:\Multiwfn_3.8_dev_bin_Win64>move /Y *.cub "C:\VMD"
C:\Multiwfn_3.8_dev_bin_Win64\density1.cub
C:\Multiwfn_3.8_dev_bin_Win64\ESP1.cub
移动了         2 个文件。

然后打开VMD软件，输入iso即可画图

再输入ext即可显示极值点

上图中的材质已被我改为transparent，因此能看到背面的极值点，修改方法可以通过Grapgics-Materials进行修改，然而答主修改时没有任何变化，于是在ESPiso.vmd中直接将EdgyGlass改成了Transparent

怎么添加刻度呢？操作如下图

接下来是自定义图片的导出，File-Render，渲染器选择Tachyon，开始渲染，然后在VMD的根目录下得到了文件vmdscene.dat，接下来调用渲染器进行渲染，渲染的.bat脚本是：

1	tachyon_WIN32.exe vmdscene.dat -aasamples 24 -mediumshade -trans_vmd -res 1024 768 -format BMP -o vmdscene.bmp

其中的1024是长，768是宽，可以改成想要的分辨率

Comment and share

从零开始的电化学阻抗谱原理

May 27, 2025 in 从零开始的科研之旅, first stage

什么是电化学阻抗谱？

有一种叫做黑箱的未知特性系统。可采用一系列测量方法，通过施加输入信号并检测输出响应来探查黑箱。例如，假设将黑箱置于暗室中，然后施加特定波长的光照。若观察到响应信号（如电流），则可判定箱内物质具有光活性。输入与输出之间的关系称为“传递函数”，阻抗谱是传递函数的一种特例。

阻抗：电路中电阻、电感、电容对交流电个阻碍作用个统称。阻抗是一个复数，实部称为电阻，虚部称为电抗。电化学过程主要涉及电荷在电极和电解质界面上的传输和存储。例如，电极表面的双电层效应会形成双电层电容，电解质中的离子扩散也会导致电容行为。这些机制使电容成为电化学系统中主要的电抗元件。相比之下，电感通常与磁场的变化相关，比如电流通过线圈时产生的磁场。然而，在常规的电化学系统中，磁场变化通常不显著，因此电感的影响非常小，通常电化学系统中考虑电阻、容抗而不考虑感抗。

电化学阻抗谱（EIS）通过施加交流电压（或电流）并测量产生的电流（或电压），能够研究电极和电化学池的电化学特性，利用等效电路进行分析，可以确定与电极电化学反应相关的重要参数，如电荷转移电阻、双电层电容和Warburg阻抗

阻抗的复数表达形式

交流电（AC）电路中的电流、电压以及阻抗通常使用复数来表示，主要原因有两个：

交流信号（以及许多其他正弦波现象）具有一个幅值和一个相位，它们分别与复数的模和幅角非常相似
复数的基本运算，如加法、减法、乘法和除法，更容易执行和编程

注意：

因为符号i 在交流电路中用于表示电流，这里我们用j作为虚数单位，即j² = −1
符号 ℜe代表复数的实部

预备知识

复数的标准形式 Z = a + bj 复数的模 $$ \lvert Z \rvert =\sqrt {a^2+b^2} $$

欧拉公式 e^ix = cos x + isin x 复数的指数形式 Z = re^jθ = rcos θ + rjsin θ 取实数部分： ℜe(re^jθ) = rcos θ

单变量复函数f(t)的导数的实部等于f(t)的实部的导数，证明如下： $$ $$ **复函数$f(t)$的不定积分的实部等于$f(t)$的实部的不定积分，证明如下** $$ $$

复数形式的阻抗

电压源

v(t) = V₀cos (ωt) = ℜe(V₀cos (ωt) + V₀jsin (ωt)) = ℜe(V₀e^jωt)

电阻

对于如下图的交变电路

定义Z_R = R，那么电流$i=\Re e(\dfrac{V_R}{R})$或$I=\dfrac{V_R}{Z_R}$

电容

电容上电流的定义为：$I=C(\dfrac{d v(t)}{dt})$，则有如下推导： $$ \begin{align} v(t)=&\dfrac{1}{C}\int i dt\\ v(t)=&V_0\cos (\omega t)=\Re e(V_0e^{j\omega t})\\ \Re e(V_0e^{j\omega t})=&\dfrac{1}{C}\int idt\\ \dfrac{d}{dt}(\Re e(V_0e^{j\omega t}))=&\dfrac{d}{dt}(\dfrac{1}{C}\int idt)\\ \Re e(j\omega V_0e^{j\omega t})=&\dfrac{1}{C}i\\ i=&\Re e(j\omega CV_0e^{j\omega t}) \end{align} $$ 电容器C两端的电压为V_C = V₀e^jωt，那么电容对电流的阻碍作用$Z_C=\dfrac{1}{j\omega C}=-\dfrac{j}{\omega C}$

与电阻的阻抗（电阻）Z_R不同，电容的阻抗（容抗）Z_C是一个纯虚数

电感

在电化学系统中，电感并不常见

电感L上的电流i和电压v(t)之间的关系为：$v(t)=L\dfrac{di}{dt}$，则有： $$ \begin{align} \Re e(V_0e^{j\omega t})=&L\dfrac{di}{dt}\\ \int \Re e(V_0e^{j\omega t})=&\int L\dfrac{di}{dt}dt\\ \Re e(\dfrac{1}{j\omega}V_0e^{j\omega t})=&Li\\ i=&\Re e(\dfrac{1}{j\omega L}V_0e^{j\omega t}) \end{align} $$ 电感的对电流的阻碍作用为：Z_L = jωL，类似的，电感的阻抗（感抗）也是一个纯虚数

交变电路中的电阻和电容

电压源为：E(t) = E₀sin (ωt)

电容表现出与电阻相似的行为，即电压与电流的幅值比恒定，但电流相位相对于外加电压存在-90°偏移:$I=E_0\omega C\sin(\omega t-\dfrac{\pi}{2})$

电阻和电容的电流响应可用两个参数表征：

阻抗模量|Z|
- 电阻：|Z_R|恒定
- 电容：$\lvert Z_C \rvert=\dfrac{1}\omega C{}$
电压与电流之间的相位差θ
- 电阻：θ = 0^∘
- 电容：θ = −90^∘

使用阻抗表达式Z = |Z|(cos θ + sin θ)，在交流测量中可以统一处理这两种电路元件

电阻和电容的串联电路比较简单，而并联的情况稍复杂，讨论一下：

施加正弦电压E(t)与并联RC电路响应电流(t)的关系。总电流是流经电阻R(t)和电容C(t)的电流之和

总阻抗： $$ \begin{align} \dfrac{1}{Z}=&\dfrac{1}{Z'}+\dfrac{1}{Z''}\\ =&\dfrac{1}{R}+j\omega C\\ Z=&\dfrac{1}{\dfrac{1}{R}+j\omega C}\\ =&\dfrac{R}{1+(\omega CR)^2}-\dfrac{\omega CR^2 j}{1+(\omega CR)^2}\\ \lvert Z \rvert=&\sqrt{Z'^2+Z''^2}\\ =&\dfrac{R}{\sqrt{1+(\omega CR)^2}} \end{align} $$ 总结一下串联和并联 RC 电路的复阻抗特性：

	Z = Z^′ + jZ^″	Z^′ = \|Z\|cos θ	Z^″ = \|Z\|sin θ	$\lvert Z \rvert =\sqrt{Z'^2+Z''^2}$	$\theta=\tan^{-1}\dfrac{Z''}{Z'}$
串联RC电路	$R-\dfrac{j}{\omega C}$	R	$-\dfrac{1}{\omega C}$	$\sqrt{R^2+\dfrac{1}{\omega^2C^2}}$	$\tan^{-1}(-\dfrac{1}{\omega RC})$
并联RC电路	$(\dfrac{1}{R}-\dfrac{\omega C}{j})^{-1}$	$\dfrac{R}{(\omega RC)^2+1}$	$\dfrac{-\omega R^2C}{(\omega RC)^2+1}$	$\dfrac{1}{\sqrt{\dfrac{1}{R^2}+\omega^2C^2}}$	tan⁻¹(ωRC)

阻抗的图解

如上所述，阻抗Z既可用阻抗模值|Z(ω)|和相位角θ(ω)表示，也可用阻抗的实部Z^′(ω)和虚部Z^″(ω)表示，分别对应两种图：Bode图和Nyquist图

(a)串联与(b)并联 RC 电路的波特图。虚线标示频率 f = 1/(2πRC)处相位角θ = −45°的位置。

(a)串联与(b)并联 RC 电路的Nyquist图。虚线标示频率 f=1/(2πRC)对应的相位角θ=-45°位置。

串联RC电路的阻抗是Z_R与Z_C之和。因此在高频区，由于Z_R ≫ Z_C ，电路阻抗等于电阻R；而在低频区，由于Z_R ≪ Z_C ，阻抗等于Z_C 。如Bode图所示，串联RC电路在高频区表现为电阻特性（|Z|与频率无关，相位角θ为 0°），在低频区表现为电容特性（|Z|与频率成反比，相位角θ为−90°）。Nyquist图中，串联 RC 电路的阻抗轨迹为一条从实轴垂直延伸的直线，与实轴交点的值等于 R。该电路的时间常数为RC，意味着当频率$ω=\dfrac{1}{RC}$时，|Z_R| = |Z_C|且相位角θ = −45°

并联RC电路的阻抗Z通过公式$\frac{1}{Z}=\frac{1}{Z_R}+\frac{1}{Z_C}$计算。因此，并联 RC 电路，阻抗较小的电路元件特性会显现。在高频区，并联RC电路的Z 等于Z_C ；而在低频区则等于Z_R，这与串联RC电路的情况相反。在并联 RC 电路的Nyquist图中，阻抗轨迹 Z 呈现为半圆形，圆弧直径等于 R。该电路的时间常数τ = RC，与串联电路相同。当频率$\omega =\frac{1}{RC}$时，|Z_R| = |Z_C|且相位角θ = −45°，此时对应Nyquist图中圆弧的顶点。

与Bode图相比，Nyquist图更适合通过轨迹形状来辨识电路结构。但需指出的是，当两个并联 RC 电路的时间常数相近时，即便其电阻和电容值不同，它们的半圆也难以区分。如下图所示，当两个并联 RC 电路的时间常数仅相差10倍时，可观察到由两段圆弧重叠形成的畸变半圆。只有当时间常数差异达到 100 倍以上时，两段圆弧才能被清晰区分。因此在构建等效电路时，仅依Nyquist图并不可取，必须验证该等效电路能否合理解释目标电化学现象中的每个基本过程

不同电容值 C 下的电路Nyquist图：(a) C=0.2 毫法，(b) 1 毫法，(c) 5 毫法，(d) 10 毫法。两个并联 RC 电路的时间常数比分别为：(a) 250，(b) 50，(c) 10，(d) 5。

下面列举一些常见的电路模型以及相应的Nyquist、Bode图吧

电感的Nyquist图在第四象限
“电阻等价于点”
“电容等价于垂直线”

电容和电阻的并联反映到Nyquist图上即为一个半圆弧
有几个电阻-电容并联电路就有可能出现几个半圆弧
电阻$\rm{R}_0$+电阻-电容并联电路就等价于将半圆弧向右平移$\rm{R_0}$个单位

定义法拉第阻抗Z_f = R_ct + Z_W

上图左侧为理想高频部分，此时电荷转移电阻R_ct很大，阻抗由电荷转移主导
上图右侧为理想低频部分，此时“电容断路（电容特性：通高频阻低频）”，阻抗由物质扩散主导（溶液离子扩散、固体电解质中离子扩散等等）
如果是一个频率非常广泛的、既包含高频又包含低频的电化学阻抗测试，那么Nyquist图中可能同时出现高频部分的半圆弧和低频部分的直线，如下图所示

电极的等效电路

在电化学阻抗谱（EIS）中，电极的电化学行为通常由包含电阻和电容的等效电路表示。其中，Randles 电路常被用作电极的等效电路模型

定义法拉第阻抗Z_f = R_ct + Z_W
溶液电阻R_s：即使使用参比电极，IR降也不可忽略。电解质传导过程中涉及的电阻称为溶液电阻（R_s ），可通过Nyquist奎斯特图中高频极限处与实轴的交点获得
电荷转移电阻R_ct：当过电位较大时，电极发生极化，过电位与电流之间呈现指数依赖，而当电极电位接近平衡电位时，过电位与电流的指数关系可简化为线性关系。因此，在极低过电位区域（<10 mV）内，电极表现为电阻特性，其电阻值称为电荷转移电阻，对应于Nyquist奎斯特图中半圆的直径。由于电荷转移电阻的概念仅适用于获得线性关系的极小过电位区域，因此电化学阻抗谱中的电压振幅通常设置为≤10 mV。
双电层电容C_dl：电极与电解质界面处的电双层类似于电容器，虽然无法直接从Nyquist图中确定C_dl 的数值，但如果确定了R_ct ，则可通过以下关系式利用圆弧顶点频率对应并联电路的时间常数计算出C_dl ：τ = R_ct × C_dl
Wurburg阻抗：源于电解液中活性物质的扩散过程。在半无限扩散条件下，Z_W可通过以下方程表示：

$$ Z_W=\dfrac{RT}{n^2F^2Ac\sqrt{D}}\cdot\dfrac{1}{\sqrt{\omega}}\cdot\dfrac{1-j}{\sqrt{2}} $$

其中，n为反应中转移的电子数，A 为电极的几何面积，c为活性物质的浓度，D为扩散系数。Wurburg阻抗可用一个电容器与电阻器无限级联的电路表示。在Wurburg阻抗中，Z’和-Z’‘相等，因此在Nyquist图中可观察到一条斜率为-45°的直线。通过绘制 Z’和 Z’’随ω变化的关系曲线，可从两条平行直线的斜率确定扩散系数。值得注意的是，该式源自能斯特方程中浓度电势依赖项($\dfrac{\partial E}{\partial c}$)。对于不符合能斯特方程的电极（如发生固态扩散的电极），必须考虑($\dfrac{\partial E}{\partial c}$)项才能计算扩散系数

如下图所示，当扩散系数 D 增大时，活性物质能快速响应电极表面浓度变化而迅速补充至电极界面。因此较大的 D 值会导致Z_W值减小。反之，Z_W随D值减小而增大，从而使直线向低频方向移动。最终在Bode图中，Z_W的线性区会与并联RC电路引起的|Z|变化区域重叠。这一现象在下下图（😛）所示的Nyquist图中更为明显：当 D 值较大时，-45°斜率的直线与圆弧可清晰区分；但随着 D 值减小，圆弧与直线开始重叠。这种重叠现象表明活性物质的扩散速度低于电荷转移反应的时间常数，本质上意味着扩散过程成为速率控制步骤。

不同 RT/(n2F2AcD^0.5 )值下 Randles电路的 Bode 图：(a) 1, (b) 10, (c) 100, (d) 300 Ω s −0.5 — 不同 RT/(n^2F2AcD^0.5 )值下 Randles电路的 Bode 图：(a) 1, (b) 10, (c) 100, (d) 300 Ω s −0.5

不同 RT/(n2F2AcD^0.5 )值下 Randles电路的Nyquist图：(a) 1, (b) 10, (c) 100, (d) 300 Ω·s −0.5 — 不同 RT/(n^2F2AcD^0.5 )值下 Randles电路的Nyquist图：(a) 1, (b) 10, (c) 100, (d) 300 Ω·s −0.5

一个性能良好的固态电解质材料，其离子电导率应当较大，扩散系数D应当较大，那么Nyquist图应当接近于a或b的形貌

从理想电极到真实电极——常相位元件CPE

CPE被定义为： $$ Z_{CPE}=\dfrac{1}{Y_O(j\omega)^n} $$

θ = 90^∘(1 − n)

其中θ是与理想情况（φ=90°）的相位偏差，n 是从 0 到 1 的常数，定义了与理想行为的偏差

n	相位角φ（理想行为）	与理想行为的偏差θ	阻抗方程
1	90	0	$Z(\omega)=\dfrac{1}{j\omega \rm{C_{dl}}}$
0.9	81	9	$Z_{\omega}=\dfrac{1}{Y_O(j\omega)^{0.9}}$
0.8	72	18	$Z_{\omega}=\dfrac{1}{Y_O(j\omega)^{0.8}}$
0	0	90	$Z_{\omega}=\dfrac{1}{Y_O(j\omega)^{0}}=\dfrac{1}{Y_O}=R$
0.5	45	45	$Z(\omega)=\dfrac{1}{Y_O}\sqrt{j\omega}$

两个极端情况：

n=1，CPE表现为理想电容器
n=0，CPE 表现为电阻器
n=0.5，$Z(\omega)=\dfrac{1}{Y_O}\sqrt{j\omega}$，即等价于 Warburg 阻抗

用传输线模型来直观理解扩散阻抗

Warburg阻抗是一个复杂的元素，它表示氧化还原物质向电极表面的质量传递，并被描绘为奈奎斯特图的低频范围上的-45°线。这种行为是指化学物质的时间相关（非稳态）半无限扩散，如下图A所示，在电极/电解质界面处存在单个边界，距离 x=0。朝向本体溶液，并且在静止条件下，扩散层扩展到由电化学电池的尺寸限定的无限长度（x→∞），同时浓度梯度随时间下降。

将扩散层视为由无数个 “电阻（R_d）- 电容（C_d）” 单元串联，其中：

R_d：单位长度的扩散电阻（反映扩散阻力）
C_d：单位长度的 “化学电容”（反映浓度存储能力）
半无限扩散——无限长传输线，阻抗为Warburg直线
有限扩散——有限长传输线，末端接电阻（投射）或电容（反射）
- 当有限扩散区对扩散物质可渗透时，随着时间和电流流过电化学电池，建立了稳态浓度梯度（dC/dx=常数）。这是透射或可渗透边界的情况。
- 当有限扩散区在一段时间后对扩散物质不渗透时，电荷转移是不可能的，并且扩散物质的浓度梯度变为零（dC/dx=0）。这是反射或不渗透边界的情况

A）电化学电池的半无限区，（B）半无限区的传输线（TL）描述，（C）t<td 跨度处的有限边界扩散，（D）透射，和（E）t>td 处的反射边界

==具有透射边界的有限长度线性扩散的典型例子包括固体氧化物或聚合物膜燃料电池，和旋转盘电极实验中的扩散层==

固态电解质的阻抗谱分析

以陶瓷固态电解质为例，假设其仅存在单离子传输而不具有电子导电性

在极高频率区域（>1 MHz）的阻抗数据对于通过 EIS 测量分析固体电解质至关重要，这与其它材料的 EIS 分析存在显著差异。若某固体电解质的体相离子电导率为10⁻⁴S ⋅ cm⁻¹ ，且离子迁移相关电容为10⁻¹¹F ⋅ cm⁻² ，则对应1/(2π × τ)（τ：时间常数）的特征频率约为∼10⁷Hz（10 MHz）。若未对该频率区域采取特殊防护措施进行EIS测量，高频区的阻抗数据将会缺失，仅能观测到部分圆弧。在此情况下，EIS数据中包含的来自电缆和电池的电子电阻分量无法消除。此外，由晶粒内部电阻和晶界电阻构成的离子电阻也无法分离，从而导致固体电解质重要信息的丢失。

由于大多数陶瓷固体电解质是由不同尺寸和取向的晶粒聚集而成的多晶材料，除了晶粒内部电阻外，离子在晶粒间界面处的传输电阻对离子电导率具有重要影响。晶粒内部电阻$R_{\rm{gi}}$ （又称体电阻）是材料固有的离子传输阻力。因此，若测试样品不含孔隙，通过晶粒内部电阻倒数计算得到的体离子电导率等于该固体电解质单晶的离子电导率。晶界电阻$R_{\rm{gb}}$ 表示离子穿越晶粒间界面（晶界）及/或位于晶界的化学组成不同微相时所受的传输阻力。下图展示了陶瓷固体电解质的等效电路及Nyquist图，该电路由两个串联的并联 RC 回路构成：一个电阻对应$R_{\rm{gi}}$ ，另一个对应$R_{\rm{gb}}$ 。由于两个并联RC电路具有不同的时间常数（τ = R × C），电化学阻抗谱能够分别测定这两个电阻值。

由两个并联RC电路组成的固体电解质等效电路，分别对应晶粒内部与晶界的串联连接。(b)是根据(a)所示等效电路计算得出的固体电解质典型Nyquist图。需注意，为简化起见图中省略了电极响应部分。

由于开发和改进陶瓷固体电解质最有效的方法是通过降低$\rm{R_{gb}}$ 使离子电导率接近单晶水平，因此分离$\rm{R_{gi}}$ 和$\rm{R_{gb}}$ 并分别评估它们至关重要。要通过高频区半圆直径确定$\rm{R_{gi}}$ ，应通过增加电解质厚度、减小电极面积或降低测量温度来降低特征频率。然而这些方法存在局限性，例如当电极面积小于样品几何面积时会干扰样品中的电流分布，因此需要测量直至高频区的阻抗数据。本文介绍了我们测量固体电解质高达 100 MHz 高频阻抗的方法。

在固体电解质的电化学阻抗谱（EIS）测量中，通常采用电极/电解质/电极的三明治结构组件。这种配置的优势在于能够实现简单的两端测量（注：设备与电极之间需连接两根电压线和两根电流线，共计四根电缆）。采用电极/电解质/电极结构进行EIS测量时，电极材料的选择至关重要。电极材料需满足以下必要条件：能够构建稳定界面且不与电解质发生化学反应，同时具备足够的电子导电性。金或铂等贵金属常被用作阻塞电极，某些情况下也会采用能与移动离子发生可逆反应的可逆电极。圆盘状固体电解质样品通常通过粉末压片成型，必要时还需进行煅烧处理。样品的相对密度最好大于 90%，理想情况下应超过 95%。

在电极涂覆工艺中，有两个关键因素会影响固体电解质的阻抗测量结果。其一是电极厚度：过薄的电极可能导致阻抗数据偏离原始值，进而造成对固体电解质离子传输行为的误判。其二是固体电解质的表面状态，该因素会显著影响阻抗测量。本节将详细阐述这两个因素如何影响阻抗数据。

下图展示了不同电极厚度及固体电解质表面是否抛光处理的氟离子导体Nyquist图。所有情况下均采用溅射法沉积直径为7mm 的铂电极。阻抗测试使用 Keysight E4990A 阻抗分析仪在室温空气环境下完成。所有Nyquist图均呈现两个特征：由晶粒内部与晶界电阻（具有相近时间常数）形成的畸变弧线（两个半圆弧融合），以及低频区对应阻塞电极与固体电解质界面电容特性的直线段。不同电极厚度的Nyquist图存在显著差异：对于 200 nm 厚电极（图 2a），离子传输电阻（晶粒内部与晶界电阻）形成的阻抗弧起始于坐标原点；而 50 nm 薄电极（图 2b）的阻抗弧则发生水平偏移，且由阻抗弧外推得到的高频极限阻抗值偏离原点。这种因电极厚度不同而产生的差异可归因于样品中的电流分布；对于较薄的电极，其横向电阻不容忽视

采用铂电极的氟离子传导固体电解质的Nyquist图：(a) 200 纳米厚度电极，(b) 50 纳米厚度电极。(c) 样品表面抛光后 50 纳米铂电极的Nyquist图。右侧为放大视图以便观察。圆圈表示测量值，虚线表示拟合曲线。

如上图，当样品经 2000 目砂纸抛光并沉积 50 纳米铂电极后，所得阻抗弧严重偏离原点。这表明抛光过程引入了除离子传输过程之外的额外电阻。此外，抛光样品的离子传导过程对应阻抗弧略有扩大，直线斜率趋于平缓。这种斜率变化反映了抛光前后样品阻塞行为的差异：抛光样品在<1 MHz 频率区间出现与阻塞行为重叠的寄生电阻。固态电解质抛光前后的阻抗行为变化归因于抛光过程中（为去除煅烧样品表面偏析层而进行的）样品表面劣化。

为消除电缆电感对高频测量的干扰，设备与电解池（样品）之间的连接电缆应尽可能缩短。建议采用特性阻抗为 50Ω的同轴电缆，并将其紧贴样品安装。此外，必须合理设置电流回流路径。对于高达 100MHz 的高频测量，需确保电缆不会引起相位差畸变。在进行高频测量前，应对系统进行开路/短路/负载校准

离子电导率的计算

$$ \sigma'(\omega)=\dfrac{L}{S}\times\dfrac{Z'(\omega)}{\vert Z'(\omega) \rvert ^2} $$

对于使用固态电解质模具进行电化学阻抗谱测试的情况，离子电导率的计算公式为： $$ \sigma =\dfrac{L}{S\times R} $$ 其中L为固态电解质pellet圆片厚度，S为圆片pellet的面积，R为电阻（Nyquist图拐点处的横坐标）

实战演练——以$\rm{Li_{5.5}PS_{4.5}Cl_{1.5}}$的电化学阻抗谱图为例

这是之前测试的图，简单画了一下，2Tons（249MPa）压力下的测试结果。可以看到，这与常见的阻抗谱差别非常大，并没有出现超高频的半圆弧，只出现了低频的扩散部分和中高频的拐点。这是收阻抗仪的最高频率所限，使用的阻抗仪最高频只能达到约1MHz左右，对于10mS/cm的固态电解质，想要完整的测出半圆弧可能需要100MHz甚至更高的频率

由于设备频率上限较低，只能得到粒子抵抗和粒界阻抗的总值
活化能分析中无法判断随温度变化特性值发生变化的是粒子阻抗还是粒界阻抗
在通过变温方式尝试分离粒子阻抗和粒界阻抗时，可控温度范围有限，且无法对变温引起的参数条件变化进行校正

==此外，交变电流源的最高频率不等于最高有效频率==，有效的最高频率常常受限于测试端子类型和连线长度，在进行变温控制的时候，==夹具及导线材料的介电常数因为温度变化而变化，如果不及时矫正阻抗会使结果产生较大偏离==

下图展示了超高有效频率下的LLTO型固态电解质的阻抗谱

参考资料

Lazanas, A. Ch.; Prodromidis, M. I. Electrochemical Impedance Spectroscopy─a Tutorial. ACS Meas. Sci. Au 2023, 3 (3), 162–193. https://doi.org/10.1021/acsmeasuresciau.2c00070.
Lazanas, A. Ch.; Prodromidis, M. I. Correction to “Electrochemical Impedance Spectroscopy─a Tutorial.” ACS Meas. Sci. Au 2025, 5 (1), 156–156. https://doi.org/10.1021/acsmeasuresciau.5c00007.
Ariyoshi, K.; Mineshige, A.; Takeno, M.; Fukutsuka, T.; Abe, T.; Uchida, S.; Siroma, Z. Electrochemical Impedance Spectroscopy Part 2: Applications. Electrochemistry 2022, 90 (10), 102008–102008. https://doi.org/10.5796/electrochemistry.22-66080.
Ariyoshi, K.; Siroma, Z.; Mineshige, A.; Takeno, M.; Fukutsuka, T.; Abe, T.; Uchida, S. Electrochemical Impedance Spectroscopy Part 1: Fundamentals. Electrochemistry 2022, 90 (10), 102007–102007. https://doi.org/10.5796/electrochemistry.22-66071.
Hu, W.; Peng, Y.; Wei, Y.; Yang, Y. Application of Electrochemical Impedance Spectroscopy to Degradation and Aging Research of Lithium-Ion Batteries. J. Phys. Chem. C 2023, 127 (9), 4465–4495. https://doi.org/10.1021/acs.jpcc.3c00033.

Comment and share

科研中的小问题

May 25, 2025 in 从零开始的科研之旅, first stage

实验部分

球磨

球磨的原理是什么？
球磨的作用是什么？
200rpm合不出来，300rpm就能合出来？

X-ray diffraction

X射线是什么？
X射线衍射的原理和用途是什么？
XRD谱图的横坐标和纵坐标分别是什么意义？
XRD为什么会产生峰？峰是什么意义？
横坐标2theta角度为什么通常选择20^∘~80^∘
单晶衍射和粉末衍射分别是什么？
关于消光？

近年来锂离子电池大规模应用的快速发展要求系统地了解锂离子电池的退化机制，这在很大程度上取决于先进表征方法的利用和发展。基于 X 射线的技术，如 X 射线吸收光谱（XAS）、X 射线衍射（XRD）和 X 射线光电子能谱（XPS），可以确定化学成分和元素的价态和晶体结构的变化。电子和扫描探针显微技术，（3）如扫描电子显微镜（SEM）、透射电子显微镜（TEM）和低温透射电镜，可以观察电极的形貌变化和固体电解质界面（SEI）的厚度。基于滴定的技术，如滴定气相色谱（TGC）（4）和质谱滴定（MST），（5,6）可以定量死 Li 金属和 SEI 组分。此外，固态核磁共振（NMR）（7）和电子顺磁共振（EPR）（8）能够实现对 Li 镀层的时间分辨和定量检测。然而，这些方法大多是侵入性技术，需要额外昂贵以及具有特殊电池设计的精密仪器进行测试，在实际应用中是不可行的。除了这些物理特性之外，电化学方法，如差分电压分析（DVA）、（9,10）差分充电电压（DCV）、（11）增量容量分析（ICA）、（12）和库仑效率（CE）、（13,14）也在电池老化研究中发挥了重要作用。然而，它们是相当复杂的确定方法，需要很长的监测时间。

电化学阻抗谱测试

为什么测出来的Nyquist图没有半圆部分呢？

理论知识

记得看Sob老师的博客，多逛逛计算化学公社呀

Comment and share

从零开始的X射线原理

May 24, 2025 in 从零开始的科研之旅, first stage

前言

1917 年，赫尔首次描述了简单粉末衍射仪的构造，这是在 1895 年威廉·康拉德·伦琴发现 X 射线后不久。衍射仪测量 X 射线反射的角度，从而获得其包含的结构信息。如今，这项技术的分辨率有了显著提高，它被广泛用作分析相信息和解决固态材料晶体结构的工具。X 射线衍射技术用于确定样品的组成或晶体结构。对于大分子和无机化合物等较大的晶体，它可用于确定样品中原子的结构。如果晶体尺寸太小，它可以确定样品的组成、结晶度和相纯度。

什么是X射线？

X射线是一种电磁波，由内壳层电子跃迁所得：

对于原子序数为Z、仅含一个电子的类氢原子，内壳层电子跃迁所辐射的能量可由Rydberg方程计算： $$ \bar{v}=(\dfrac{1}{n_i^2}-\dfrac{1}{n_f^2})RZ^2 $$ 可以明显看出，与电子跃迁相关的能量差随着原子序数的增加而大幅增大，并且在这种跃迁过程中发射的辐射波长随着原子序数的增加从 10⁻⁷ m (1000 Å) 范围移至 10⁻¹⁰ m (1 Å)范围（此时辐射的电磁波被定义为X射线，此辐射类型被称为特征辐射）、

如下图，要引发这种内壳层跃迁，需要产生一个电子空位，必须从K层中移除一个电子。由高压电场加速过的阴极电子书轰击靶材，将其部分能量传递给靶材料的电子，从而导致电子激发。如果入射电子的能量足够高，一些电子可能会击出靶中 K 壳层的一个电子，从而产生一个空位

不同能级跃迁回K层对应着不同的谱线，如$\rm{K_{\alpha}}$线和$\rm{K_{\beta}}$线。对于铜的X射线光谱，在低能分辨率下只有2条特征线被观察到，但是，在更高的分辨率下，很容易看出Kα线是双峰，它被标记为$\rm{K_{α1}}$和$\rm{K_{α2}}$。因为铜的2p轨道，能级L分裂为$\rm{L_{II}}$和$\rm{L_{III}}$，间距非常小（0.020kev），波长Kα1 (= 1.54056 Å) 和Kα2 (= 1.54439 Å) 也非常接近。因此，通常说的$\rm{K_{\alpha}}$其实是$\rm{K_{α1}}$和$\rm{K_{α2}}$的加权平均值

阳极	$\rm{K_{\alpha}}$	$\rm{K_{\beta}}$
铜	1.54184Å	1.39222Å
钼	0.71073Å	0.63229Å

以铜靶为例 $$ E=\dfrac{hc}{\lambda}=\dfrac{6.626\times10^{-34}\times2.998\times 10^8}{1.542\times 10^{-10}\times1.6\times10^{-19}}=8.05~\rm{kV} $$ 8.05 kV为铜的$\rm{K_{\alpha 1}}$线的最低激发电压，实际操作时通常选择40kV。激发电压选择过低会造成谱线不明显、电流过大热损耗大等问题，选择40kV是长期实验优化的结果。

由于种种原因，如量子力学的选择定则：电子跃迁需要满足角动量变化Δl=±1、散热和成本等问题，通常X射线光谱中不会出现$\rm{K_{\gamma}}$甚至更高的线，即使能够发生N→K的跃迁（比如高Z元素W)，也会因为强度极弱而湮灭在噪声中，因此X射线光谱中通常见到的是$\rm{K_{\alpha}}$和$\rm{K_{\beta}}$

此外，在高速电子轰击靶材的过程中除了由内壳层电子跃迁所辐射的特征X射线，还存在着轫致辐射：高速电子被靶材原子核电场偏转，动能转化为连续电磁辐射。

如下图，当加速电压（单位：kV）较小、高速电子无法引起内壳层跃迁时（加速电压小于临界电压），那么此时原子的辐射均为轫致辐射

钼靶 X 射线光谱随施加电压的变化关系 $$ V_{临界}=\dfrac{hc}{e\lambda_{K-edge}} $$

在所得的一系列X射线中，只有满足布拉格条件（波长与晶面间距d匹配）的射线才可用于样品检测，因此需要通过滤光片等手段只保留窄带波长，去除了除 K_α X 射线之外的所有电磁波

XRD谱图上的峰是什么？

XRD——X-ray diffraction（X射线衍射），其本质是波的干涉与叠加。

当波（如光波、声波、X射线）遇到障碍物或穿过周期性结构时，波前会发生弯曲和重新分布，这种现象称为衍射。障碍物或结构的尺寸与波长相当，这是发生衍射的关键条件，正因如此，与原子间距匹配的X射线适合用作发生衍射的电磁波。X射线是电磁波，照射到晶体时，每个原子中的电子受X射线的作用，成为次级波源，向四周散射X射线（这种散射也被称为相干散射），这些不同路径的X射线相遇并发生干涉（衍射发生的本质是弹性散射波，（即此时的次级X射线）的干涉）

干涉分为两类：

相长干涉（Constructive Interference）：满足布拉格条件时，所有晶面的散射波相位一致，两列波波峰对齐，叠加后振幅增强 → 形成峰。
相消干涉（Destructive Interference）：波峰与波谷对齐，叠加后振幅抵消 → 形成暗区。

布拉格定律： nλ = 2dsin θ 其中λ是外加波长，θ是衍射角，d 是原子平面之间的距离。在大多数情况下，我们假定n=1，即一级衍射，因为总可以将 n=2,3,…的衍射峰解释为来自(nh nk nl) 晶面的衍射——即来自晶面间距为$\rm{d_{hkl}}$的$\dfrac{1}{n}$倍的晶面$\rm{d_{nhnknl}}$。然后，原子平面之间的距离可用于确定成分或晶体结构。

当路径 ABC 与 A’B’C’之间的距离相差整数个波长（λ）时，衍射的 X 射线表现出相长干涉

X 射线衍射的结果绘制出了在各自的 2θ位置处不同衍射角的信号强度。2θ位置对应于样品中晶体或原子之间的特定间距，该间距由入射到样品中的 X 射线束的衍射角决定。峰的强度与该相或具有该间距的分子数量有关。峰的强度越大，具有该特定间距的晶体或分子数量就越多。

消光是什么？为什么只有特定的晶面处才会出现峰？

衍射峰强度取决于晶胞内所有原子散射辐射之间的相位关系，但经常出现这种情况：虽然布拉格定律预测某个峰应该存在，但其实际强度却为零（这是因为布拉格定律不涉及原子位置，仅与晶胞的大小和形状有关）。

例如，对于具有体心立方BCC晶胞的(100)晶面衍射峰强度：相位关系表明，在晶胞顶面和底面（即(100)晶面）散射的 X 射线虽然发生相长干涉，但与晶胞中心原子散射的X射线存在180^∘的相位差，因此最终强度为零。下表给出了不同立方晶格中特定衍射峰出现与否的选择定则。

Bravais Lattice	Reflections Present	Reflections Absent
Simple Cubic	All	None
Body-Centered Cubic	(h+k+l)= even	(h+k+l)= odd
Face-Centered Cubic	h,k,l unmixed (either all odd or all even)	h,k,l mixed
Base-Centered Cubic	h,k unmixed (either all even or all odd)	h,k mixed

为什么谱图的横坐标是2θ？

实验中，样品和探测器的转动是围绕同一轴进行的，因此 总偏转角的计算需要结合样品和探测器的相对运动，如下图所示

峰宽度的意义是什么？

峰的宽度与晶体尺寸成反比。较窄的峰对应较大的晶体。较宽的峰意味着可能存在较小的晶体、晶体结构缺陷，或者样品在本质上可能是非晶态的，即一种缺乏完美结晶度的固体。对于较小的样品，将由X 射线衍射（XRD）所得的图谱与数据库中标准样品的谱图对比可确定样品的组成。此结论可由谢乐公式（Debye-Scherrer）直接给出 $$ D=\dfrac{K\lambda}{B \cos\theta} $$

K为Scherrer常数，也称形状因子，若B为衍射峰的半峰高宽，则K=0.89；
D为晶粒垂直于晶面方向的平均厚度（Å）；
B为实测样品衍射峰半高宽度
θ为布拉格衍射角，单位为角度
λ为X射线]波长，对于Cu kα来说，一般为1.54056 Å

为什么晶相能相长干涉形成峰，而非晶相不能？

晶体是长程有序的周期性结构，满足布拉格条件，即所有晶面(hkl)的的间距d恒定，当入射角θ满足布拉格公式时，所有晶面的散射波相位一致 → 相长干涉 → 尖锐峰。

非晶体可能短程有序，而长程无序，d值不固定，散射波来自不同局部区域，相位杂乱，无法形成尖锐的峰

举个例子

例子仅供参考

晶体硅的晶胞参数a=0.543nm，X射线波长为0.154nm，那么我们可以预测某个晶面会在何处出现峰

对于立方晶系： $$ d_{hkl}=\dfrac{a}{\sqrt{h^2+k^2+l^2}} $$ 以(111)晶面为例 $$ d_{111}=\dfrac{0.543}{\sqrt{3}}=0.314nm $$ 根据布拉格公式： $$ \theta =\arcsin (\dfrac{\lambda}{2d_{111}})=14.3^\circ $$

如何确定2theta角度？

虽然理论上2θ范围是0°~180°，但实际操作中极少覆盖全部范围，原因包括：

物理限制：大多数晶面间距d有一个最小值，那么$\theta_{\rm{max}} \le \arcsin(\dfrac{\lambda}{2d_{\rm{min}}})$
信号强度：高角度区（如2θ>120°）的衍射强度极弱，信噪比差。
时间成本：全范围扫描耗时过长，而关键峰通常集中在20°~80°

一般选择20~80°，一般情况

常规晶体材料的晶面间距d大多在1~10Å之间

Cu靶X射线波长λ≈1.54Å，代入公式：
- 当d=10Å时，( $\sinθ = λ/(2d) ≈0.077$ → θ≈4.4° → 2θ≈8.8°)
- 当d=1Å时，( $\sinθ = 0.77$ → θ≈50° → 2θ≈100°)

参考资料

Sadoway, D. Introduction to Solid State Chemistry.
Reig, A. CHEM322: Inorganic Chemistry.

Comment and share

关于电子能谱

May 18, 2025 in 从零开始的科研之旅, first stage

什么是电子能谱？

电子能谱就是用一定能量的==电子束或光子与样品表面相互作用==，使样品表面原子中不同能级的电子激发成自由电子，并研究这些自由电子的电子强度（电子数目）按其能量变化的分布曲线。根据使用的激发源的不同，电子能谱又分为==X射线光电子能谱（XPS）、紫外光电子能谱（UPS）、俄歇电子能谱（AES）==

电子能谱的基础是光电效应：当能量为hν的单色光照射到样品上时，样品中的原子或分子吸收能量从而使电子脱离样品成为自由电子（即光电子）逸出，与此同时还能产生==俄歇电子和荧光X射线等== hν = E_b + E_k + E_r 其中E_b为原子的始态和终态能量差，可以看成发射的光电子的结合能；E_r为原子的反冲能量，可忽略

X射线光电子能谱

X射线光电子能谱（XPS），也称为化学分析电子能谱（ESCA）。价电子对X光子的光电效应截面远小于内层电子，所以XPS主要用于研究内层电子的结合能。由于内层电子不参与成键，保留了较多的原子性质，所以常用于表面组分的化学分析，也可以进一步应用于确定这些==元素的化学或电子状态==。它允许以非破坏性方式测定样品的原子组成，由于X射线仅能穿透样品 5 – 20 Å，因此可以进行表面特异性分析。

电子在暴露于足够能量的电磁辐射（波长在X射线范围的高能光子）时会从材料表面射出，即光电效应： KE = hν − E_b − φ 其中E_b是电离能（也叫电子结合能），φ 是光谱仪的功函数。

不同的原子具有不同的电子壳层和不同的结合能，因此XPS可用于识别样品表面的元素组成。计算发射的具有各种E_k的光电子的数量，并生成XPS光谱，这有助于识别除H和He之外的所有状态的所有元素。

Instrumentation of X-ray photoelectron spectroscopy

如上图所示，X射线光电子能谱仪器的主要部件包括以下部分：

X射线源
超高真空（UHV）
电子透镜
XPS能量分析仪
检测器
监视器

首先，高速电子撞击金属阳极使其发生内壳层跃迁，辐射出特征X射线。特征X射线脱离射线源部分后进入到UHV超高真空腔，在这里，约10⁻⁹~10⁻¹⁰ torr甚至更低的压强可以减少XPS光谱中的噪声量并避免样品表面污染产生任何额外峰。接着，X射线在UHV环境中照射到样品表面，样品释放出光电子，光电子被电子透镜所汇聚并减速。大量减速后的光电子接着进入XPS能量分析仪，这是一个两端施加有电压的半球腔体。多通道检测器可以检测出电子的动能，类似于胶片。高性能计算机用作监视器，用于控制镜头和分析仪的电压，它还会产生 XPS 光谱——计算机自动绘制E_k强度与结合能的关系生成的。

定性分析

原子的电子能级可分为两种类型——与原子核结合的核能级和仅弱结合的价能级。原子的价能级是与其他原子的价能级相互作用以在分子和化合物中形成化学键的能级。它们的特性（价键特性）和能量通过这个过程发生变化，成为形成的新物种的特性。

原子的核心能级电子的能量几乎与原子结合的化学物质无关，因为它们不参与键合过程。因此，核心水平结合能的识别因此提供了元素的独特特征。元素周期表中的所有元素都可以通过这种方式识别（除了H和He，它们没有核心能级）。

定量分析

定量分析需要测量相对峰强度，峰的强度取决于表面膜中存在的元素的原子浓度（N）。

化学状态分析

化学位移是同一分子或两种不同化合物中两种化学不同状态的同一原子之间的结合能差异。它是样品电离原子周围的原子性质和原子数的特征。它可以用来查看样品表面原子的化学键是什么样的。

它涉及化学键性质、官能团的性质和原子的氧化态。
XPS 的一个重要功能是它能够区分不同化学状态的相同元素，因为电子的 BE 会因氧化引起的元素周围电荷密度的变化而不同

深度分析

对薄表面薄膜中元素深度分布的估计称为深度表面分析。这样做是为了知道薄膜是否均匀。

优势与局限性

优势：

可以分析所有样品，包括固体、气体和液体样品，但大多数是固体样品。
定量效果好
出色的化学状态测定能力
它适用于从生物材料到金属的各种材料。
这是一种非破坏性技术。

局限性：

缺乏良好的空间分辨率，即空间分辨率差
仅中等绝对灵敏度
它无法检测氢和氦
它不适用于痕量分析。

为什么常常用Al作为靶材

X射线衍射（XRD）实验中选择的是铜靶或者钼靶，为什么这里用的是铝靶呢？

XPS对X射线源的光子能量和能量分辨率（单色性/线宽）有非常特殊的要求。XPS的核心是精确测量样品表面原子发射的光电子的动能，进而计算其结合能。结合能的微小偏移（< 1 eV）往往对应着元素的化学态变化（如氧化态、化学键）。==如果X射线本身的线宽很宽，就会导致光电子谱峰变宽，使得相邻的、化学位移很小的峰无法分辨，严重影响化学态分析的精度==

XRD测量的是晶面间距引起的衍射峰位置（角度）。衍射峰的自然宽度通常比X射线源的线宽要宽得多（受晶体尺寸、应力等因素影响）。因此，X射线源本身的线宽（如Cu Kα约2.5 eV）对衍射峰位置测量的影响相对较小，强度更重要。

Al Kα的能量是1486.6 eV，这个能量足够高，能够有效地激发周期表中绝大多数元素（从轻元素如Li到重元素）的核心能级电子（如C 1s, O 1s, N 1s, Fe 2p, Au 4f等），这些能级是XPS分析化学态的关键。能量不是越高越好。能量过高（如Cu Kα的8048 eV）会产生更强的韧致辐射本底和更多的二次电子，增加谱图的背景噪音。同时，高能光子的光电离截面（电离效率）在特定能量范围外会降低。1486.6 eV是一个在激发效率和背景噪音之间取得良好平衡的能量点。

为什么光电方程中有功函数这一项？

什么是功函数呢？功函数 Φ 是指将一个电子从==固体材料==内部（==费米能级处==）移动到真空中静止状态（真空能级）所需的最小能量。想象电子在固体内部处于一个“能量坑”的底部（费米能级）。要把这个电子从坑里拉出来扔到坑外平坦的地面上（真空能级），你需要克服坑壁的阻挡做功。这个“坑壁的高度”就是功函数 Φ（单位：电子伏特）。

为什么光电方程中有功函数这一项？理解的关键在于区分 ==“原子内结合能 (Eb)”==和==“电子逸出样品进入真空并被检测器测到”== 这两个过程所涉及的能量参考点。

Einstein的原始光电方程是针对于自由原子或者气体的，Eb 是电子相对于原子真空能级的结合能（即把电子从原子轨道移到离原子无穷远处所需能量）：

对于孤立原子，真空能级是明确的参考点
测得的光电子动能直接反应原子内电子的结合能E_b

但是对于==固体==，情况变得更加复杂了。固体中，原子不是孤立的，电子处于==周期性势场中==。固体中的电子能量通常参考费米能级 (Fermi Level, EF)——绝对零度下电子占据的最高能级。总结一下固体中的光电效应：内壳层电子吸收光子能量跃迁到无穷远轨道处（依然处在固体体相中），这一步需要克服结合能，从无穷远轨道处跃迁到真空层（脱离表面）还需要克服一部分能量，这部分能量被称为功函数φ。

所以功函数补偿了谱仪测量参考点（其费米能级）与理论计算参考点（真空能级）之间的差异。它确保了测量动能 KE 与样品内结合能 Eb（相对于样品费米能级）之间的定量关系成立。

XPS怎么确定元素的氧化态？

graph TD
A[氧化态升高] --> B[原子周围电荷密度降低]
B --> C[原子实有效正电荷增加]
C --> D[内壳层电子受核束缚增强]
D --> E[结合能E_b升高]
E --> F[XPS峰向高结合能方向移动]

XPS的测试项目

全谱
- 得出元素的成分信息，也可以对所测元素进行半定量计算
精细谱
- 判断化学态，确定不同化学态的百分比含量
刻蚀/溅射
- 去除样品表面杂志后进行分析
- 深度剖析，探测不同深度下的样品成分
价带谱
- 主要测量价带顶，也可以分析价带结构，探测外层价电子信息，也可以用来区分化学环境价态
俄歇谱
- 主要测量价带顶，也可以分析价带结构，探测外层价电子信息，也可以用来区分化学环境价态
Mapping
- 表征不同元素成分在特定区域的分布情况，分辨率〈EDS，XPS深度10nm左右，EDS深度1um左右

样品制备

粉末
- 样品量20mg以上，压片制样，粘到双面胶带上分析；颗粒越细越好，且分散均匀；在压片样品表面分析时尽量选用平整测试区域，越平整信号越强

Comment and share

使用Multiwfn加载ELFCAR文件绘制等值面图

May 18, 2025 in 从零开始的科研之旅, VASP

本文仅为学习记录，并非教程

Multiwfn是一个波函数分析程序，功能非常强大，也能绘制ELF等值面图，它支持多种类型的输入文件，既有量子化学计算程序的输出文件，又有第一性原理计算程序的输出文件

以下是Multiwfn加载ELFCAR绘制等值面图的过程

step1：用ELFCAR绘图其实是基于格点数据文件的绘图，需要修改根目录下的配置文件settings.ini，将参数iuserfunc改为-1或者-3
step2：加载ELFCAR文件
step3：4 Output and plot specific property in a plane
step4：100 User-defined function (iuserfunc=
step5：5 Shaded relief map with projection
step6：1: XY plane
step7：1a

如果将Multiwfn的根目录添加到PATH环境变量以方便调用，那么画出来的图可能会遇到如下如题

如果直接在根目录下双击.exe文件运行的话就能正常绘图

Comment and share

掺杂质量计算的python脚本

May 17, 2025 in 从零开始的科研之旅, VASP

起因是我想合成各式各样的材料，需要用到各种药品，计算起来比较麻烦，尤其是非化学计量比的化合物比如Li_5.3PS_4.3Cl_0.7Br，遂用AI写了一个python脚本，个人觉得还是比较实用，脚本文件和相对原子质量的文本文件需要在同一文件夹内

from decimal import Decimal, getcontext
import re
from sympy import symbols, Matrix, linsolve, Rational
from collections import defaultdict

# Improve Decimal precision
getcontext().prec = 100

class FormulaParser:
    @staticmethod
    def parse(formula: str) -> dict:
        elements = {}
        pattern = re.compile(r"([A-Z][a-z]?)([\d.]*)?")
        idx = 0
        while idx < len(formula):
            m = pattern.match(formula, idx)
            if not m:
                raise ValueError(f"Invalid formula at: {formula[idx:]}")
            elem, num = m.group(1), m.group(2)
            coeff = Decimal(num) if num else Decimal(1)
            elements[elem] = elements.get(elem, Decimal(0)) + coeff
            idx = m.end()
        return elements

class AtomicMassReader:
    @staticmethod
    def read(filename="atomic_masses.txt") -> dict:
        masses = {}
        try:
            with open(filename, 'r') as f:
                for line in f:
                    line = line.strip()
                    if not line or line.startswith('#'):
                        continue
                    parts = line.split()
                    if len(parts) != 2:
                        raise ValueError(f"Line format error: {line}")
                    masses[parts[0]] = Decimal(parts[1])
        except FileNotFoundError:
            raise FileNotFoundError(f"File {filename} not found")
        if not masses:
            raise ValueError("Atomic masses file empty")
        return masses

class SynthesisCalculator:
    def __init__(self, atomic_masses: dict):
        self.atomic_masses = atomic_masses
        self.target = {}
        self.target_mm = Decimal(0)
        self.reagents = []
        self.reagent_formulas = []
        self.molar_masses = []

    def input_target(self):
        f = input("Enter target formula: ").strip()
        self.target = FormulaParser.parse(f)
        for e, c in self.target.items():
            if e not in self.atomic_masses:
                raise ValueError(f"Mass of element {e} not defined")
            self.target_mm += c * self.atomic_masses[e]
        print(f"Parsed: {self.target}")
        print(f"Target molar mass = {self.target_mm.normalize()} g/mol")

    def input_reagents(self):
        n = int(input("Enter number of reagents: "))
        for i in range(n):
            f = input(f"Reagent {i+1} formula: ").strip()
            self.reagent_formulas.append(f)
            parsed = FormulaParser.parse(f)
            mm = sum(parsed[e] * self.atomic_masses[e] for e in parsed)
            self.reagents.append(parsed)
            self.molar_masses.append(mm)
            print(f"{f} molar mass = {mm.normalize()} g/mol")

    def validate(self):
        tgt = set(self.target)
        prov = set().union(*(r.keys() for r in self.reagents))
        miss = tgt - prov
        extra = prov - tgt
        if miss:
            print(f"Error: Missing elements: {miss}")
            exit()
        if extra:
            print(f"Warning: Extra elements: {extra}")

    def solve_basis(self):
        # Construct A, b corresponding to 1 mol 
        elems = sorted(set(self.target) | set().union(*(r.keys() for r in self.reagents)))
        A = [[r.get(e, Decimal(0)) for r in self.reagents] for e in elems]
        b = [self.target.get(e, Decimal(0)) for e in elems]
        # Convert to rational number to avoid floating point error
        A_rat = [[Rational(str(val)) for val in row] for row in A]
        b_rat = [Rational(str(val)) for val in b]
        vars = symbols(f'x0:{len(self.reagents)}')
        sol = linsolve((Matrix(A_rat), Matrix(b_rat)), vars)
        if not sol:
            print("Error: No exact solution for 1 mol target.")
            exit()
        tup = next(iter(sol))
        # Convert rational to decimal exact representation
        base = {}
        for i, r in enumerate(tup):
            base[vars[i]] = Decimal(r.p) / Decimal(r.q)
        return base, elems

    def calculate(self):
        t = input("Input type (0 mass g, 1 moles): ").strip()
        if t not in ['0','1']:
            print("Invalid type")
            exit()
        v = Decimal(input("Amount: ").strip())
        n = v / self.target_mm if t=='0' else v
        print(f"Target moles = {n.normalize()}")

        base, elems = self.solve_basis()
        print("\nReagent requirements:\n" + "="*40)
        total_mass = Decimal(0)
        for i, var in enumerate(sorted(base.keys(), key=lambda x: int(str(x)[1:]))):
            moles = base[var] * n
            mass = moles * self.molar_masses[i]
            total_mass += mass
            print(f"Reagent {i+1} ({self.reagent_formulas[i]}):")
            print(f"  Moles: {moles.normalize():.10f} mol")
            print(f"  Mass:  {mass.normalize():.10f} g\n")
        print(f"Total mass = {total_mass.normalize():.10f} g")

        #Verify element balance
        print("Verification:")
        for e in elems:
            actual = sum((base[symbols(f'x{i}')] * n) * r.get(e, Decimal(0))
                         for i, r in enumerate(self.reagents))
            target_amt = self.target.get(e, Decimal(0)) * n
            print(f"{e}: target {target_amt.normalize():.10f}, actual {actual.normalize():.10f}")


def main():
    print("High-Precision Compound Synthesis Calculator")
    print("="*60)
    try:
        masses = AtomicMassReader.read()
        print(f"Loaded elements: {', '.join(masses.keys())}")
    except Exception as e:
        print(e)
        exit()
    calc = SynthesisCalculator(masses)
    calc.input_target()
    calc.input_reagents()
    calc.validate()
    calc.calculate()

if __name__ == "__main__":
    main()

附一张使用截图，后面也会放到Github上面

Comment and share

初学ELF电子定域化函数的计算与分析

May 14, 2025 in 从零开始的科研之旅, VASP

本文仅为学习记录，并非教程

这次记录文主要侧重于应用，参考了Sob老师的博文和网络上的各大教程

以FCC Cu为例进行ELF计算与分析

采用FCC面心立方最密堆积的铜单胞，使用vasp计算，先进行结构优化，确认收敛后进行ELF计算

Global Parameters
ISTART =  1            (Read existing wavefunction, if there)
ISPIN  =  2            (Non-Spin polarised DFT)
MAGMOM = 4*0.6
# ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
 ENCUT  =  530        (Cut-off energy for plane wave basis set, in eV)
# PREC   =  Accurate   (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid, helps GGA convergence)
LASPH  = .TRUE.        (Give more accurate total energies and band structure calculations)
PREC   = Accurate      (Accurate strictly avoids any aliasing or wrap around errors)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells, be careful)
# LPLANE = .TRUE.      (Real space distribution, supercells)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGXF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGYF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGZF    = 300        (FFT grid mesh density for nice charge/potential plots)

Electronic Relaxation
ISMEAR =  1            (Gaussian smearing, metals:1)
SIGMA  =  0.2         (Smearing value in eV, metals:0.2)
NELM   =  90           (Max electronic SCF steps)
NELMIN =  6            (Min electronic SCF steps)
EDIFF  =  1E-08        (SCF energy convergence, in eV)
# GGA  =  PS           (PBEsol exchange-correlation)

Ionic Relaxation
NSW    =  0         (Max ionic steps)
IBRION =  -1            (Algorithm: 0-MD, 1-Quasi-New, 2-CG)
ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-02        (Ionic convergence, eV/AA)
# ISYM =  2            (Symmetry: 0=none, 2=GGA, 3=hybrids)
LELF = .TRUE.
NPAR = 1

在优化完的结构上进行自洽计算，添加LELF = .TRUE.这个参数，确保ELFCAR文件输出

官方wiki还提到必须显示开启NPAR = 1，然而网上的各大教程都没有提到这一点

If LELF is set, NPAR=1 has to be set explicitely in the INCAR file in addition

且有一篇讨论帖：https://wiki.vasp.at/forum/viewtopic.php?t=20020

不过我自己尝试的时候发现超算平台上添加NPAR = 1会报错，而在自己电脑上则不会报错，有点迷、、、

用Vesta载入ELFCAR文件

先设置一下等值面数值

Isosurface Level（等值面水平）是一个关键参数，用于控制可视化电子局域化程度的阈值，VESTA会绘制出ELF等于该值的等值面（isosurface），直观展示电子局域化程度高于该阈值的区

Properties-Isosurfaces-Isosurfaces level，设置完后可以按以下Tab键就会自动保存应用

这里如果设置0.75的话啥也看不到，因为铜单胞中存在着金属键，并没有很强的局域化电子，而是均匀分布的电子气

高定域性通常出现在共价作用区域、孤对电子区域、原子内核及壳层结构区域

然后就可以查看ELF图了

首先是3D的图

左上角-Edit-Lattice planes，新建一个图，根据想查看的面的来设置米勒指数hlk，以及设置距离d

然后是2D的图

左上角-Utilities-2D Data Display-Slice即可

可以看到，顶点Cu与顶点Cu原子之间、顶点Cu和面心Cu之间的电子的局域化程度比较高并且没有明显的偏向偏离

接下来看一下一维时的电子定域化情况

左上角-Utilities-Line profile即可

首先棱上的情况

Cu和Cu两原子中间存在均匀的非零的电子局域函数，虽然不足0.5但分布平坦，这是金属键合的特征

两个铜原子的ELF(r)值相近

面对角线上的情况

三个铜原子的ELF(r)值相近

离子晶体CsF

再看一下离子键的情况，以CsF为例

氟和铯之间出现深蓝色，这说明在这个区域电子的分布较为离域化，即氟和铯之间不存在共用电子，键的类型偏向离子键

再看看对角线上的电荷分布

两端为铯中间为氟，可以看到氟和铯ELF(r)值差别较大

共价晶体——硅

如箭头所指，硅原子和硅原子之间有较强的电子定域化作用，键的类型为共价键

如下图，体对角线上的ELF图也能很好证明，Si-Si之间的区域电子化定域程度极高，呈现明显的尖峰

共价晶体——灰锡

我仍记得高中时课本讲的一个稍微有些特殊的例子——灰锡，灰锡是一个共价晶体

可以看到，灰锡与晶体硅结构上类似，Sn($\frac{1}{4},\frac{1}{4},\frac{1}{4}$)与Sn(0,0,0)与以及Sn(0.5,0.5,0)成键，Sn($\frac{3}{4},\frac{3}{4},\frac{1}{4}$)与Sn(0.5,0.5,0)和Sn(1,1,0)成键。成键电子局域化程度较高，为共价电子

Comment and share

2025年5月Intel平台VASP.6.4.3的编译安装

May 12, 2025 in 从零开始的科研之旅, VASP

2025年5月Intel平台VASP.6.4.3的编译安装

WSL的安装

选择Arch Linux，个人比较熟悉，仅供学习使用不作为稳定的生产力环境，配置如下，两颗E5 2676 v3

打开powershell或者终端，输入wsl --install archlinux以自动安装，这一步需要使用代理(MAGIC)，主机上使用Clash-Verge进行代理（运气好也能连上），终端中输入：

1 2	set http_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，http代理 set https_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，https代理

如果没有条件使用代理，也可以手动安装，参考Arch Wiki，https://wiki.archlinuxcn.org/wiki/%E5%9C%A8_WSL_%E4%B8%8A%E5%AE%89%E8%A3%85_Arch_Linux，手动下载.wsl镜像，然后安装

笔者计算机C盘性能太差，于是想要将wsl2迁移到另一个分区

安装成功

[root@NAS ~]# fastfetch
                  -`                     root@NAS
                 .o+`                    --------
                `ooo/                    OS: Arch Linux x86_64
               `+oooo:                   Host: Windows Subsystem for Linux - archlinux (2.4.12)
              `+oooooo:                  Kernel: Linux 5.15.167.4-microsoft-standard-WSL2
              -+oooooo+:                 Uptime: 5 hours, 16 mins
            `/:-:++oooo+:                Packages: 155 (pacman)
           `/++++/+++++++:               Shell: bash 5.2.37
          `/++++++++++++++:              WM: WSLg 1.0.65 (Wayland)
         `/+++ooooooooooooo/`            Terminal: xterm-256color
        ooosssso++osssssso+`           CPU: Intel(R) Xeon(R) E5-2676 v3 (48) @ 2.39 GHz
       .oossssso-````/ossssss+`          Memory: 1.09 GiB / 15.57 GiB (7%)
      -osssssso.      :ssssssso.         Swap: 304.95 MiB / 4.00 GiB (7%)
     :osssssss/        osssso+++.        Disk (/): 18.67 GiB / 1006.85 GiB (2%) - ext4
    /ossssssss/        +ssssooo/-        Disk (/mnt/c): 159.05 GiB / 237.01 GiB (67%) - 9p
  `/ossssso+/:-        -:/+osssso+-      Disk (/mnt/d): 13.16 TiB / 14.55 TiB (90%) - 9p
 `+sso+:-`                 `.-/+oso:     Disk (/mnt/e): 1.70 TiB / 3.64 TiB (47%) - 9p
`++:.                           `-/+/    Disk (/mnt/f): 49.45 GiB / 476.94 GiB (10%) - 9p
.`                                 `/    Disk (/mnt/g): 25.16 TiB / 29.10 TiB (86%) - 9p
                                         Local IP (eth0): 172.31.36.99/20
                                         Locale: en_US.UTF-8

编译器的安装

pacman -S make
pacman -S cmake
pacman -S gcc
pacman -S gcc-fortran
####发现gcc和gfortran编译器太新，会编报错(C23)，故降级软件包
#添加中科大源来安装yay
pacman -Syy yay
yay -S downgrade #安装降级工具
downgrade -S gcc-libs #选择12.1.0版本，安装gcc的依赖
downgrade -S gcc #安装gcc，选择12.1.0版本
downgrade -S gcc-fortran #安装gforran，选择12.1.0版本

前往Intel官网下载oneAPI BaseTool kithttps://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html?packages=oneapi-toolkit&oneapi-toolkit-os=linux&oneapi-lin=online

1 2	chmod +x xxxxxxxxxx.sh #赋予下载下来的sh文件可执行权限 xxxx.sh #运行sh文件

oneAPI BaseTool kit不需要全部安装，只需安装MKL部分，

HPC也只需安装一部分

此处博主踩了较多的坑：

编译器采用2022版，不可采用过新编译器
- 2024年10月后，Intel删除了ifort编译器
gcc编译器不可太新，否则编译某些fortran模块会报错，博主从14降级到了12.1.0

VASP.6.4.3的编译安装

源码包可以去各大私密性较强的论坛或者微信公众号获取，仅供个人学习使用，不用于学术、商业用途

选择arch/目录中的模板makefile.include.intel，复制到.(上一层目录)，并修改名称为makefile.include，修改目录如下

# Default precompiler options
CPP_OPTIONS = -DHOST=\"LinuxIFC\" \
              -DMPI -DMPI_BLOCK=8000 -Duse_collective \
              -DscaLAPACK \
              -DCACHE_SIZE=4000 \
              -Davoidalloc \
              -Dvasp6 \
              -Duse_bse_te \
              -Dtbdyn \
              -Dfock_dblbuf

CPP         = fpp -f_com=no -free -w0  $*$(FUFFIX) $*$(SUFFIX) $(CPP_OPTIONS)

FC          = mpiifort
FCL         = mpiifort

FREE        = -free -names lowercase

FFLAGS      = -assume byterecl -w

OFLAG       = -O2 -xhost  #添加-xhost标记
OFLAG_IN    = $(OFLAG)
DEBUG       = -O0

OBJECTS     = fftmpiw.o fftmpi_map.o fftw3d.o fft3dlib.o
OBJECTS_O1 += fftw3d.o fftmpi.o fftmpiw.o
OBJECTS_O2 += fft3dlib.o

# For what used to be vasp.5.lib
CPP_LIB     = $(CPP)
FC_LIB      = $(FC)
CC_LIB      = icc
CFLAGS_LIB  = -O
FFLAGS_LIB  = -O1
FREE_LIB    = $(FREE)

OBJECTS_LIB = linpack_double.o

# For the parser library
CXX_PARS    = icpc
LLIBS       = -lstdc++

##
## Customize as of this point! Of course you may change the preceding
## part of this file as well if you like, but it should rarely be
## necessary ...
##

# When compiling on the target machine itself, change this to the
# relevant target when cross-compiling for another architecture
VASP_TARGET_CPU ?= -xHOST
FFLAGS     += $(VASP_TARGET_CPU)

# Intel MKL (FFTW, BLAS, LAPACK, and scaLAPACK)
# (Note: for Intel Parallel Studio's MKL use -mkl instead of -qmkl)
FCL        += -qmkl=sequential
MKLROOT    ?= /opt/intel/oneapi/mkl/2022.1.0 #一定要添加MKLROOT的目录，否则会报错


/home/storm/intel/oneapi/mkl/2022.1.0
LLIBS      += -L$(MKLROOT)/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64
INCS        =-I$(MKLROOT)/include/fftw

# HDF5-support (optional but strongly recommended)
#CPP_OPTIONS+= -DVASP_HDF5
#HDF5_ROOT  ?= /path/to/your/hdf5/installation
#LLIBS      += -L$(HDF5_ROOT)/lib -lhdf5_fortran
#INCS       += -I$(HDF5_ROOT)/include

# For the VASP-2-Wannier90 interface (optional)
#CPP_OPTIONS    += -DVASP2WANNIER90
#WANNIER90_ROOT ?= /path/to/your/wannier90/installation
#LLIBS          += -L$(WANNIER90_ROOT)/lib -lwannier

然后开始编译

发现使用命令make all -j48会编译报错

于是改为官方wiki的命令make DEPS=1 -j48，编译成功

测试

make test进行测试，有一些报错，懒得管了、、、能跑就行

Comment and share

Vaspkit的使用

May 12, 2025 in 从零开始的科研之旅, VASP

安装、配置好vaspkit之后输入vaspkit可启动软件包

(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/trash-vaspkit$ vaspkit
            \\\///
           / _  _ \         Hey, you must know what you are doing.
         (| (o)(o) |)       Otherwise you might get wrong results.
 o-----.OOOo--()--oOOO.------------------------------------------o
 |         VASPKIT Standard Edition 1.5.1 (27 Jan. 2024)         |
 |         Lead Developer: Vei WANG (wangvei@icloud.com)         |
 |      Main Contributors: Gang TANG, Nan XU & Jin-Cheng LIU     |
 |  Online Tutorials Available on Website: https://vaspkit.com   |
 o-----.oooO-----------------------------------------------------o
        (   )   Oooo.                          VASPKIT Made Simple
         \ (    (   )
          \_)    ) /
                (_/
 ===================== Structural Utilities ======================
 01) VASP Input-Files Generator    02) Mechanical Properties
 03) K-Path for Band-Structure     04) Structure Editor
 05) Catalysis-ElectroChem Kit     06) Symmetry Analysis
 07) Materials Databases           08) Advanced Structure Models
 ===================== Electronic Utilities ======================
 11) Density-of-States             21) Band-Structure
 23) 3D Band-Structure             25) Hybrid-DFT Band-Structure
 26) Fermi-Surface                 28) Band-Structure Unfolding
 31) Charge-Density Analysis       42) Potential Analysis
 44) Piezoelectric Properties      51) Wave-Function Analysis
 62) Magnetic Analysis             65) Spin-Texture
 68) Transport Properties
 ======================== Misc Utilities =========================
 71) Optical Properties            72) Molecular-Dynamics Kit
 74) User Interface                78) VASP2other Interface
 84) ABACUS Interface              91) Semiconductor Kit
 92) 2D-Material Kit               95) Phonon Analysis
 0)  Quit
 ------------>>

以下是中文翻译对照：

===================== 结构工具 ======================
01) VASP 输入文件生成器
02) 力学性质
03) 能带结构的 K 路径
04) 结构编辑器
05) 催化 - 电化学套件
06) 对称性分析
07) 材料数据库
08) 高级结构模型
==================== 电子工具 ======================
11) 密度态
21) 能带结构
23) 三维能带结构
25) 混合 DFT 能带结构
26) 费米面
28) 能带结构展开
31) 电荷密度分析
42) 势分析
44) 压电性质
51) 波函数分析
62) 磁性分析
65) 自旋纹理
68) 输运性质
====================== 杂项工具 ========================
71) 光学性质
72) 分子动力学套件
74) 用户界面
78) VASP 转其他格式接口
84) ABACUS 接口
91) 半导体套件
92) 二维材料套件
95) 声子分析
0) 退出

我目前还是初学者小白，使用得最多的功能还是生成输入文件功能尤其是生成INCAR文件，

 ------------>>
01
 ==================== VASP Input Files Options ===================
 101) Customize INCAR File
 102) Generate KPOINTS File for SCF Calculation
 103) Generate POTCAR File with Default Setting
 104) Generate POTCAR File with User Specified Potential
 105) Generate POSCAR File from cif (no fractional occupations)
 106) Generate POSCAR File from Material Studio xsd (retain fixes)
 107) Reformat POSCAR File in Specified Order of Elements
 108) Successive Procedure to Generate VASP Files and Check
 109) Submit Job Queue

 0)   Quit
 9)   Back
 ------------>>
  ------------>>
101
 +---------------------------- Tip ------------------------------+
 |          WARNNING: You MUST know what wou are doing!          |
 |Some Parameters in INCAR file need to be set/adjusted manually.|
 +---------------------------------------------------------------+
 ======================== INCAR Options ==========================
 ST) Static-Calculation            SR) Standard Relaxation
 MG) Magnetic Properties           SO) Spin-Orbit Coupling
 D3) DFT-D3 no-damping Correction  H6) HSE06 Calculation
 PU) DFT+U Calculation             MD) Molecular Dynamics
 GW) GW0 Calculation               BS) BSE Calculation
 DC) Elastic Constant              EL) ELF Calculation
 BD) Bader Charge Analysis         OP) Optical Properties
 EC) Static Dielectric Constant    PC) Decomposed Charge Density
 PH) Phonon-Calculation            PY) Phonon with Phononpy
 NE) Nudged Elastic Band (NEB)     DM) The Dimer Method
 FQ) Frequence Calculation         LR) Lattice Relaxation
 MT) Meta-GGA Calculation          PZ) Piezoelectric Calculation

 0)   Quit
 9)   Back
 ------------>>
 Input Key-Parameters (STH6D3 means HSE06-D3 Static-Calcualtion)

以下是中文对照：

+---------------------------- 提示 ------------------------------+
|          警告：你必须清楚自己在做什么！          |
|INCAR 文件中的一些参数需要手动设置/调整。|
+---------------------------------------------------------------+
======================= INCAR 选项 ==========================
ST) 静态计算            SR) 标准松弛
MG) 磁性性质            SO) 自旋轨道耦合
D3) DFT-D3 无阻尼校正  H6) HSE06 计算
PU) DFT+U 计算          MD) 分子动力学
GW) GW0 计算            BS) BSE 计算
DC) 弹性常数            EL) ELF 计算
BD) Bader 电荷分析      OP) 光学性质
EC) 静态介电常数        PC) 分解电荷密度
PH) 声子计算            PY) 使用 Phonopy 的声子计算
NE) 强制弹性带（NEB）   DM) 二聚体方法
FQ) 频率计算            LR) 晶格松弛
MT) Meta-GGA 计算       PZ) 压电计算

LR晶格弛豫就对应着结构优化，当然生成的INCAR文件可能仍需要进行一些修改，比如ENCUT可能需要调整为赝势中ENMAX的1.3倍

生成KPOINTS文件

对于非能带的计算，选择102) Generate KPOINTS File for SCF Calculation使用程序自动撒点即可，但是需要用户选择撒点方式和撒点密度，选择Gamma Scheme，使用0.04精度

关于生成 KPOINTS 文件

在 VASP 计算中，KPOINTS 文件用于指定 k 点网格的设置。对于非能带计算，比如进行能量计算、优化结构等，通常不需要像计算能带结构那样沿特定的高对称路径采样 k 点，而是可以在整个布里渊区均匀地分布 k 点，这就是所说的 “自动撒点”。

撒点方式

Monkhorst-Pack 方式 ：这是最常用的自动撒点方式。它通过指定一组三个整数（nx, ny, nz），在三个倒易晶格方向上均匀地分布 k 点，形成的网格可以系统地逼近布里渊区的积分。例如，对于简单晶格，若指定 k 点网格为 3×3×3，程序就会按照 Monkhorst-Pack 方式在这三个方向各取 3 个点，生成一个均匀分布的 k 点网格。

Gamma 方式 ：这种方式会在 k 点网格中包含 Γ 点（布里渊区的中心点）。对于金属等具有高对称性的材料，在 Γ 点附近采样很重要，这时候 Gamma 方式撒点更合适。比如指定 k 点网格为 4×4×4 并采用 Gamma 方式撒点，生成的网格会以 Γ 点为中心，向周围均匀分布其他 k 点。

K 点密度

K 点密度指的是在布里渊区单位体积内分布的 k 点数量，它决定了 k 点网格的疏密程度。较高的 K 点密度意味着更精细地采样布里渊区，能够更精确地计算系统的物理性质，但也会增加计算量和资源消耗。

例如，对于一个计算资源充足的大型超算任务，为了获得高精度结果，可以选择较大的 K 点密度，如 8×8×8 甚至更高；而对于初步探索或计算资源有限的情况，可以适当降低 K 点密度，如 4×4×4。

生成KPOINTS的同时，POTCAR也会被自动生成，前提是在vaspkit的配置文件中正确设置了赝势的路径和赝势稳健的种类（默认type为PBE）

Comment and share

NEWER POSTS
OLDER POSTS
page 2 of 8

目标分子的结构优化

绘图前的准备工作

画图

什么是电化学阻抗谱？

阻抗的复数表达形式

预备知识

复数形式的阻抗

交变电路中的电阻和电容

阻抗的图解

下面列举一些常见的电路模型以及相应的Nyquist、Bode图吧

电极的等效电路

从理想电极到真实电极——常相位元件CPE

用传输线模型来直观理解扩散阻抗

固态电解质的阻抗谱分析

离子电导率的计算

实战演练——以$\rm{Li_{5.5}PS_{4.5}Cl_{1.5}}$的电化学阻抗谱图为例

参考资料

实验部分

球磨

X-ray diffraction

电化学阻抗谱测试

理论知识

前言

什么是X射线？

XRD谱图上的峰是什么？

消光是什么？为什么只有特定的晶面处才会出现峰？

为什么谱图的横坐标是2θ？

峰宽度的意义是什么？

为什么晶相能相长干涉形成峰，而非晶相不能？

举个例子

如何确定2theta角度？

参考资料

什么是电子能谱？

X射线光电子能谱

定性分析

定量分析

化学状态分析

深度分析

优势与局限性

为什么常常用Al作为靶材

为什么光电方程中有功函数这一项？

XPS怎么确定元素的氧化态？

XPS的测试项目

样品制备

以FCC Cu为例进行ELF计算与分析

先设置一下等值面数值

然后就可以查看ELF图了

首先是3D的图

然后是2D的图

接下来看一下一维时的电子定域化情况

离子晶体CsF

共价晶体——硅

共价晶体——灰锡

2025年5月Intel平台VASP.6.4.3的编译安装

WSL的安装

编译器的安装

VASP.6.4.3的编译安装

测试

生成KPOINTS文件

John Doe