where the symbol $\mathcal{J}_{ij}$ denotes the exchange interaction matrix between two spins, $S_i$ and $S_j$, the matrix $\mathcal{A}$ represents the single-ion anisotropy. 四态法计算磁耦合，需要计算4种不同的磁基态，allowing the extraction of individual components for the exchange matrix.
合适设置 INCAR, POTCAR and KPOINTS 用于磁性自洽计算, 可以用 Sym4state.jl 产生所有的输入文件用于计算最近邻交换作用和单离子各向异性，过程如下：
pre_and_post
1 2 3 4 5 6 7 8
using Sym4state cd("CrI3") do # hide Sym4state.pre_process( "./POSCAR", [24], # Take Cr element as magnetic 5.0 # There exists an interaction between atoms within a distance of 5 Å. ) end # hide
它将会构建超胞，满足两个任意原子无相互作用。给定的单层 ${CrI3}$ with a cutoff radius of 5 Å, a $2 \times 2 \times 1$ supercell will provide sufficient size. The supercell diagram below labels all the ${Cr}$ atoms:
Within the 5 Å cutoff radius, the monolayer of \ce{CrI3} exhibits two distinct groups of interactions. The first group corresponds to interactions between nearest neighbors, whereas the second group pertains to interactions arising from single-ion anisotropy. It is important to note that all atom pairs within the same group are considered equivalent. This equivalence implies the existence of symmetric operations that can transform one interaction matrix into another, highlighting the underlying symmetry of the system.
这个 shell 脚本创建 Slurm job array 计算所有 11 magnetic configurations, while efficiently managing computational resources by allowing a maximum of 2 jobs to run simultaneously.
Post-process
计算完成后 post_process function 提取不同磁结构的能量，最终创建交换矩阵。
pre_and_post
1 2 3 4 5 6 7 8 9
cd("CrI3") do # hide global pair_mat, coeff_array # hide mv("../oszicar.tar.gz", "./oszicar.tar.gz") # hide run(`tar -xvzf oszicar.tar.gz`) # hide for (idx, dir_name) in enumerate(readlines("cal_dir_list")) # hide cp("oszicar/OSZICAR_$(idx)", dir_name * "OSZICAR") # hide end # hide pair_mat, coeff_array = Sym4state.post_process("./cal.jld2") end # hide
我们可以测试 the dimensions of pair_mat and coeff_array, 它分别存储了不同原子对的起始点和结束点的索引及其对应的相互作用矩阵。
pre_and_post
1 2
size(pair_mat) size(coeff_array)
因此我们看到共有8个相互作用在cutoff radius of 5 Å. 让我们检查pair_mat中的一个特定条目，它包含表示原子对的索引：
In the aforementioned code snippet, we have configured a simulated annealing simulation, commencing at a temperature of 150 K and progressively reducing it to 0 K in steps of 2 K. The simulation operates on a $128 \times 128$ supercell of ${CrI3}$ using the previously computed interaction matrix. To assess the system, we perform a preliminary equilibration phase consisting of 100000 sweeps, followed by a measurement phase comprising 100000 sweeps for acquiring physical quantities. It is worth noting that the magnetic field is absent, rendering the magmom_vector inconsequential.
With the created mcconfig, one can initiate a Monte Carlo simulation as follows:
The parameter backend can be configured to employ CUDABackend() provided by CUDA.jl or any other backends supported by KernelAbstractions.jl to enhance performance utilizing the GPU.
The MCConfig can also be stored into a .toml file by:
pre_and_post
1 2 3
cd("CrI3") do # hide Sym4state.MC.save_config("CrI3.toml", mcconfig) end # hide
or it can also be restored by:
pre_and_post
1 2 3
cd("CrI3") do # hide mcconfig = Sym4state.MC.load_config("CrI3.toml") end # hide
cd("CrI3") do Sym4state.MC.save_config("CrI3.toml", mcconfig) end# 隐藏
或者也可以通过以下方式恢复：
1 2 3
cd("CrI3") do mcconfig = Sym4state.MC.load_config("CrI3.toml") end# 隐藏
函数
reduce_interact_mat_for_a_pair
supercell_check
pre_process
post_process
rm("CrI3", recursive=true)
参考文献
Xiang, H. J., et al. “Predicting the spin-lattice order of frustrated systems from first principles.” Physical Review B 84.22 (2011): 224429. ↩
Šabani, D., C. Bacaksiz, and M. V. Milošević. “Ab initio methodology for magnetic exchange parameters: Generic four-state energy mapping onto a Heisenberg spin Hamiltonian.” Physical Review B 102.1 (2020): 014457. ↩
Xiang, Hongjun, et al. “Magnetic properties and energy-mapping analysis.” Dalton Transactions 42.4 (2013): 823-853. ↩
SAXIS
默认值: SAXIS = (0, 0, 1)
描述
设置全局自旋量子化轴相对于笛卡尔坐标系的方向。
SAXIS 指定由泡利矩阵生成的自旋子空间的相对方向。默认情况下，SAXIS 为 (0, 0, 1)，即自旋量子化轴沿 z 轴方向。
默认方向为 $\sigma_1=\hat x$，$\sigma_2=\hat y$，$\sigma_3 = \hat z
Automag
An automatic workflow software for calculating the ground collinear magnetic state of a given structure and for estimating the critical temperature of the magnetically ordered to paramagnetic phase transition. 用于计算共线磁基态的自动工作流程软件，并用于估计磁相变温度。
Installation安装
Automag is meant to be run on a computing cluster. The first step in the installation is to clone the repository Automag 旨在运行在计算集群上。第一步 installation 是克隆存储库
I assume that you have Python 3 installed on your system, accessible through the python command, with the wheel and the venv packages installed. After cloning, go to the automag directory which has appeared and initialize a Python virtual environment 我假设您的系统上安装了 Python 3，可通过 python 命令，安装了 wheel 和 venv 包。 克隆后，转到已出现的 automag 目录并初始化 Python 虚拟环境
1 2
cd automag python -m venv .venv
After that, activate your virtual environment. If you are working on an Unix-like system, you can do it with the command 之后，激活您的虚拟环境。如果您正在使用类 Unix 系统，您可以使用命令
source .venv/bin/activate
Finally, you need to install all the necessary Python dependencies for Automag with the command 最后，您需要为 Automag 安装所有必要的 Python 依赖项 使用命令
pip install -r requirements.txt
Before using Automag, make sure that enumlib is installed on your system and that that your command line has access to the commands enum.x and makeStr.py. In addition, Automag needs to know how to use VASP, so you need to edit the file automag/ase/run_vasp.py for Automag to correctly load the MKL and MPI libraries (if needed) and call the vasp_std executable on your system. Then add the following lines to your ~/.bashrc file 在使用 Automag 之前，请确保系统上已安装 enumlib 并且 您的命令行可以访问命令 enum.x 和 makeStr.py。 此外，Automag 需要知道如何使用 VASP，因此您需要编辑文件 automag/ase/run_vasp.py 用于 Automag 正确加载 MKL 和 MPI 库 （如果需要）并调用系统上的 vasp_std 可执行文件。然后添加以下几行添加到您的 ~/.bashrc 文件
obviously replacing /PATH/TO with the actual path to these files or directories. The pp folder is the VASP pseudopotential library and it should contain two subfolders named, respectively, potpaw_LDA and potpaw_PBE. 显然将 /PATH/TO 替换为这些文件或目录的实际路径。 pp 文件夹是 VASP 赝势库，它应该包含两个子文件夹分别命名为 potpaw_LDA 和 potpaw_PBE。
Before starting to use Automag, you should also make sure to have the FireWorks library correctly pointing to a MongoDB database and configured for launching jobs through a queue management system (for more detailed information refer to the FireWorks documentation). Last but not least, open the file automag/common/SubmitFirework.py and edit line 23 with the location of your my_launchpad.yaml file. Now you are ready to use Automag. 在开始使用 Automag 之前，您还应该确保拥有 FireWorks 库正确指向 MongoDB 数据库并配置为启动作业（有关更多详细信息，请参阅 添加到 FireWorks 文档）。最后打开文件 automag/common/SubmitFirework.py 并编辑第 23 行编辑 my_launchpad.yaml 文件。现在，您可以开始使用 Automag。
Convergence tests收敛测试
When studying a magnetic structure, you may want to start from convergence tests. Automag allows you to check for convergence of the VASP parameter ENCUT, which is the energy cut-off of the plane wave basis set, and to simultaneously check for convergence of the two parameters SIGMA and kpts, which are the electronic smearing parameter and the k-mesh resolution parameter for Brillouin zone sampling, respectively. 在研究磁性结构时，您可能希望从收敛测试开始。 Automag 允许您检查 VASP 参数 ENCUT 的收敛性，该参数是平面波基组的能量截止，并同时检查 SIGMA 和 kpts 这两个参数的收敛，它们是电子的用于布里渊区采样的涂抹参数和 k-mesh 分辨率参数， 分别。
In order to run convergence tests, go to the folder 0_conv_tests and insert the following input parameters in the file input.py: 要运行收敛测试，请转到文件夹 0_conv_tests 并插入文件 input.py 中的以下输入参数：
mode can be “encut” for convergence tests with respect to ENCUT or “kgrid” for convergence tests with respect to SIGMA and kpts; mode可以是 “ENCUT” 用于关于 ENCUT 或 “kgrid” 的收敛测试 用于 SIGMA 和 kpts 的收敛测试;
poscar_file is the name of the file in POSCAR format which contains the input geometry and which has been put in the folder automag/geometries; poscar_file 是包含输入的 POSCAR 格式文件的名称 geometry 的 ，并且已放入文件夹 automag/geometries;
params is a collection of VASP parameters to be used during single-point energy calculations. params 是要在单点期间使用的 VASP 参数的集合 能量计算。 In addition, the following optional parameters can be specified: 此外，还可以指定以下可选参数：
magnetic_atoms contains the atomic types to be considered magnetic (defaults to transition metals); magnetic_atoms 包含要被视为磁性的原子类型（默认值 过渡到过渡金属）;
configuration is the magnetic configuration to use for convergence tests (defaults to ferromagnetic high-spin); configuration 是用于收敛测试的磁性配置 （默认为铁磁高自旋）;
encut_values contains the trial values for ENCUT (defaults to the interval [500, 1000] eV at steps of 10 eV); encut_values 包含 ENCUT 的 trial 值（默认为区间 [500， 1000] eV，步长为 10 eV）;
sigma_values contains the trial values for SIGMA (defaults to the interval [0.05, 0.20] eV at steps of 0.05 eV); sigma_values包含 SIGMA 的试验值（默认为区间 [0.05， 0.20] eV，步长为 0.05 eV）;
kpts_values contains the trial values for kpts (defaults to the interval [20, 100] A^-1 at steps of 10 A^-1). kpts_values 包含 kpts 的试验值（默认为区间 [20, 100]A^-1 以 10 A^-1 为步长）。 Once you have inserted the input parameters in the file input.py, launch the script 1_submit.py using the python executable of your virtual environment and a number of workflows containing single-point VASP calculations will be saved in your remote database. These calculations need to be run in a separate directory called CalcFold. You can enter it and then submit the jobs with the following commands, being in the automag directory 在文件 input.py 中插入输入参数后，启动使用虚拟环境的 python 可执行文件1_submit.py，以及 许多包含单点 VASP 计算的工作流将保存在您的远程数据库。这些计算需要在单独的目录中运行 称为 CalcFold。您可以输入它，然后使用以下内容提交作业 命令， 位于 automag 目录中
1 2
cd CalcFold nohup qlaunch -r rapidfire -m 10 --nlaunches=infinite &
This will enter the directory CalcFold and it will invoke the qlaunch command in the background, which constantly checks for new calculations in the remote database and submits them to the queue management system of your cluster. The command line option -m 10 means that qlaunch will allow a maximum number of 10 jobs to be in the queue of your system at the same time. If 10 or more jobs are waiting or running in your queue, qlaunch will not submit more jobs until their total number becomes less than 10. If you wish, you can change this parameter to a different number. In the following, I assume that the qlaunch process is always working in the background. 这将进入 CalcFold 目录，并将调用 qlaunch 命令，该命令会不断检查 remote 数据库，并将其提交到集群的队列管理系统。 命令行选项 -m 10 表示 qlaunch 将允许最大数量 的 10 个作业同时位于您的系统队列中。如果 10 或更多作业正在等待或正在您的队列中运行，qlaunch 不会提交更多作业 直到它们的总数小于 10。如果您愿意，您可以更改此 参数设置为不同的数字。在下文中，我假设 qlaunch 进程始终在后台运行。
After all calculations have been completed, you can launch the script named 2_plot_results.py in the 0_conv_tests folder. It will read the output file that Automag wrote in CalcFold and it will produce a plot of the parameters under study versus energy. In addition, the script will also print on screen the values of the parameters under study for which the error in energy is less than 1 meV/atom with respect to the most accurate result. 完成所有计算后，您可以启动名为 2_plot_results.py 0_conv_tests 文件夹中。它将读取输出文件 Automag 在 CalcFold 中写入的，它将生成参数图 正在研究与能源。此外，脚本还将在屏幕上打印能量误差小于 1 meV/atom 相对于最准确的结果。
Calculation of the electronic correlation parameter U by linear response通过线性响应计算电子相关参数 U
The linear response formalism for the calculation of the Hubbard U is based on the application of a series of small perturbations to the first magnetic atom. During the whole process, the perturbed atom must be treated independently from the other atoms of the same species, and we achieve this using a simple trick: we change the chemical identity of the atom to which we want to apply perturbations to a dummy atomic species, but we place the POTCAR file of the original atomic species in the pp folder corresponding to the dummy atom. For example, if we are calculating the value of the Hubbard U for Fe in the system Fe12O18 and we choose Zn as dummy atom, we need to ensure that the same POTCAR file of Fe is used also for Zn in order to have, instead of Zn, a Fe atom that is treated independently from the others. This can be achieved with the following commands 用于计算 Hubbard U 的线性响应形式基于对第一个磁性原子施加一系列小扰动。 在整个过程中，必须独立于同一物种的其他原子，我们用一个简单的技巧来实现这一点：我们更改我们要施加扰动的原子的化学性质添加到虚拟原子种类中，但我们将原始原子的 POTCAR 文件物种。例如，如果我们是 计算系统 Fe12O18 中 Fe 的哈伯德 U 值，我们选择 Zn 作为虚拟原子，我们需要确保也使用相同的 Fe 的 POTCAR 文件 对于 Zn，以便有一个独立处理的 Fe 原子，而不是 Zn 从其他人那里。这可以通过以下命令来实现
1 2 3
cd $VASP_PP_PATH/potpaw_PBE/Zn mv POTCAR _POTCAR cp ../Fe/POTCAR .
In this way the Automag code, and in particular the ASE library, will treat the system as if it was ZnFe11O18, but when they will look for the POTCAR file in the Zn folder, they will find the POTCAR of Fe, so the system will remain Fe12O18. Before launching this calculation, go to the directory 1_lin_response and set the necessary parameters in the file input.py: 通过这种方式，Automag 代码，特别是 ASE 库，将处理系统，就好像它是 ZnFe11O18 一样，但是当他们将在 Zn 文件夹中，他们会找到 Fe 的 POTCAR，因此系统将保持 Fe12O18。 在启动此计算之前，请转到目录 1_lin_response 并设置 文件中的必要参数 input.py：
poscar_file is the name of the file in POSCAR format which contains the input geometry and which has been put in the folder automag/geometries; poscar_file 是包含输入的 POSCAR 格式文件的名称 geometry 的 ，并且已放入文件夹 automag/geometries;
dummy_atom is the name of the dummy atomic species to use for the atom which is subject to perturbations (do not forget to manually put the right POTCAR file in the pp folder for this atom); dummy_atom 是用于原子的虚拟原子种类的名称，其中 会受到干扰（不要忘记手动放置正确的 POTCAR 文件 在此原子的 pp 文件夹中）;
dummy_position is the position of the dummy atom in your POSCAR file (from 0 to N - 1, where N is the number of atoms in the unit cell); dummy_position 是虚拟原子在 POSCAR 文件中的位置（从 0 到 N - 1，其中 N 是晶胞中的原子数）;
perturbations contains the values of the perturbations to apply to the chosen atom in eV; perturbations 包含要应用于所选 原子在 eV 中; params is a collection of VASP parameters to be used during * single-point energy calculations. params 是要在单点期间使用的 VASP 参数的集合能量计算。
In addition, the following optional parameters can be specified: 此外，还可以指定以下可选参数：
magnetic_atoms contains the atomic types to be considered magnetic (defaults to transition metals); magnetic_atoms 包含要被视为磁性的原子类型（默认值 过渡到过渡金属）;
configuration is the magnetic configuration to use for U calculation (defaults to ferromagnetic high-spin). configuration 是用于 U 计算的磁性配置（默认值 到铁磁高自旋）。
Once the input parameters have been inserted in the file input.py, you can launch the script 1_submit.py in order to save the necessary VASP jobs to the remote database. You will see that the instance of qlaunch which is running in the background will immediately send these jobs to the queue management system of your cluster. When all calculations are completed, you will find in CalcFold the file charges.txt, containing the amount of electrons on the partially occupied shell of the chosen atom for each value of the applied perturbation, for both the selfconsistent and the non-selfconsistent runs. Now you can execute the script 2_plot_results.py, which will plot the selfconsistent and the non-selfconsistent responses, it will interpolate them as straight lines to the least squares and it will calculate their slopes. The value of U is obtained from U = 1/X - 1/X0, where X and X0 are the selfconsistent and the non-selfconsistent slopes, respectively. 将输入参数插入文件 input.py 后，您可以 启动脚本1_submit.py，以便将必要的 VASP 作业保存到 远程数据库。您将看到在 后台会立即将这些作业发送到队列管理系统 的集群。完成所有计算后，您将在 CalcFold 中找到 该文件charges.txt，其中包含部分 所选原子的占用壳，用于施加扰动的每个值， 对于自洽和非自洽游程。现在您可以执行 脚本2_plot_results.py，它将绘制 selfconsistent 和 nonself-align 响应，它会将它们作为直线插入到 最小二乘法，它将计算它们的斜率。U 的值由 U = 1/X - 1/X0，其中 X 和 X0 是自洽和非自洽 slopes 的 Slope 值。
Search for the most stable magnetic state寻找最稳定的磁态
The search for the ground collinear magnetic state consists in generating a number of trial configurations and in computing their single-point energy, in order to determine which is the most thermodynamically stable. The trial configurations differ from each other only by the choice of the unit cell and by the initialization of the magnetic moments. Note that the value of the magnetic moment on each atom can change during the single-point energy calculation. Automag generates trial configurations by separately initializing each Wyckoff position occupied by magnetic atoms in a ferromagnetic (FM), antiferromagnetic (AFM) or non magnetic (NM) fashion, taking into account all possible combinations. A completely non-magnetic (NM) configuration is also generated. In this way, overall ferrimagnetic (FiM) states are allowed if the magnetic atoms occupy more than one Wyckoff position. For each magnetic atom, one or two absolute values for the magnetization can be given in input. In the first case, the given value is used for initializing all the spin-up and spin-down states in the configurations generated by Automag. Conversely, if two separate values are given for high-spin (HS) and low-spin (LS) states, then each Wyckoff position occupied by that magnetic atom is separately initialized in a HS or LS fashion, taking into account all possible combinations. In order to run such a calculation, go to the directory 2_coll and set the necessary input parameters in the file input.py: 对共线磁基态的搜索包括生成一个 试验配置的数量和计算它们的单点能量，在 order 来确定哪个是最热力学稳定的。审判 配置彼此之间唯一的区别在于晶胞的选择和 磁矩的初始化。请注意，磁性 在单点能量计算过程中，每个原子上的矩都会发生变化。 Automag 通过单独初始化每个 Wyckoff 来生成试验配置 铁磁 （FM） 中磁性原子占据的位置， 反铁磁 （AFM） 或非磁性 （NM） 方式，同时考虑所有可能的组合。 还会生成完全无磁性 （NM） 配置。这样， 如果磁性原子占据更多，则允许整体亚铁磁性 （FiM） 状态 比一个 Wyckoff 位置。对于每个磁性原子，一个或两个 磁化强度可以在 input 中给出。在第一种情况下，给定值为 用于初始化配置中的所有 spin-up 和 spin-down 状态 由 Automag 生成。相反，如果为 high-spin 给定两个单独的值 （HS） 和低自旋 （LS） 状态，则每个 Wyckoff 位置都由该 磁性原子以 HS 或 LS 方式单独初始化，同时考虑到 所有可能的组合。要运行此类计算，请转到目录 2_coll并设置必要的输入参数：
poscar_file is the name of the file in POSCAR format which contains the input geometry and which has been put in the folder automag/geometries; poscar_file 是包含输入的 POSCAR 格式文件的名称 geometry 的 ，并且已放入文件夹 automag/geometries; supercell_size is the maximum supercell size for generating distinct magnetic configurations, in multiples of the input structure; supercell_size 是产生不同磁性元件的最大超级单元大小 配置，以输入结构的倍数; spin_values a maximum of two values (HS and LS) of the magnetization in Bohr magnetons for each magnetic atom, used to initialize spin-up and spin-down states; spin_values玻尔的磁化强度最多两个值（HS 和 LS） 每个磁性原子的磁子，用于初始化自旋上升和自旋下降状态; params is a collection of VASP parameters to be used during single-point energy calculations. params 是要在单点期间使用的 VASP 参数的集合 能量计算。 In addition, the following optional parameter can be specified: 此外，还可以指定以下可选参数：
lower_cutoff is the minimum value in Bohr magnetons to which a magnetic moment can fall in order for the corresponding configuration to be used for estimating the critical temperature of the material (defaults to zero). lower_cutoff 是磁矩达到的玻尔磁子的最小值 可以进行 材料的临界温度（默认为零）。 Once the input parameters have been inserted in the file input.py, you can launch the script 1_submit.py in order to save the necessary VASP jobs to the remote database. The running instance of qlaunch will send these jobs to the queue management system of your cluster. Once all calculations have been completed, you can launch the script 2_plot_results.py which will produce a number of files containing the histogram plot of the obtained energies for all trial configurations that successfully completed the single-point energy calculation. The script also prints on screen the name of the configuration with lowest energy. 将输入参数插入文件 input.py 后，您可以 启动脚本1_submit.py，以便将必要的 VASP 作业保存到 远程数据库。正在运行的 qlaunch 实例会将这些作业发送到 集群的队列管理系统。一旦所有计算都已完成 complete，您可以启动脚本2_plot_results.py，这将生成一个 包含所有所得能量的直方图的文件数 成功完成单点能量的 Trial 配置 计算。该脚本还会在屏幕上打印配置的名称 以最低的能量。
Calculation of the critical temperature 临界温度的计算 Automag can calculate the critical temperature of the magnetically ordered to paramagnetic phase transition from a Monte Carlo simulation of the effective Hamiltonian which describes the magnetic interaction. Automag computes the coupling constants of the corresponding Heisenberg model and provides all the necessary input files to run the Monte Carlo simulation with the VAMPIRE software package. The simulation itself needs to be run by the user, while Automag can be used to plot and to fit the results. The accuracy of the Heisenberg model is evaluated by computing the Pearson Correlation Coefficient (PCC) between the DFT energies and the predicted energies of a control group of magnetic configurations. It is worth noting that this approach can be applied only if all magnetic atoms have the same absolute value of the magnetization. In order to estimate the critical temperature with Automag, go to the folder 3_monte_carlo and set the necessary parameters in the file input.py: Automag 可以计算出磁力排序的临界温度 顺磁相变来自蒙特卡洛模拟的有效 描述磁相互作用的哈密顿量。Automag 计算 耦合常量，并提供所有 使用 VAMPIRE 软件运行 Monte Carlo 模拟所需的输入文件 包。仿真本身需要由用户运行，而 Automag 可以 用于绘制和拟合结果。海森堡模型的精度为 通过计算 DFT 能量和磁控制组的预测能量 配置。值得注意的是，只有在所有 磁性原子具有相同的磁化绝对值。为了 使用 Automag 估计临界温度，转到 3_monte_carlo 并在 file input.py 中设置必要的参数：
configuration is the name of the magnetic configuration to use for the Monte Carlo simulation (usually you want to put here the name of the configuration with lowest energy, obtained at the previous step); configuration 是用于 Monte 的磁性配置的名称 Carlo 模拟（通常您希望在此处使用 最低能量，在上一步获得）; cutoff_radius is the maximum distance between two magnetic atoms to be considered as interacting neighbors; cutoff_radius 是两个磁性原子之间的最大距离 被视为互动邻居; control_group_size is the relative size of the control group used to evaluate the accuracy of the Heisenberg model; control_group_size 是用于评估的对照组的相对大小 海森堡模型的准确性; append_coupling_constants is a boolean flag which tells Automag whether or not to append the computed values of the coupling constants to the file input.py, which are needed by the script 2_write_vampire_ucf.py. append_coupling_constants 是一个布尔标志，它告诉 Automag 是否 要将耦合常量的计算值附加到文件 input.py， 脚本2_write_vampire_ucf.py需要。 In addition, the following optional parameter can be specified: 此外，还可以指定以下可选参数：
magnetic_atoms contains the atomic types to be considered magnetic (defaults to transition metals). magnetic_atoms 包含要被视为磁性的原子类型（默认值 过渡到过渡金属）。 Once the input parameters have been inserted in the file input.py, you can launch the script 1_coupling_constants.py. It will compute the coupling constants for the given cutoff radius and it will print their values on screen. In addition, it will create a file model.png, which contains a plot of the Heisenberg model energies versus the DFT energies for all the configurations in the control group. The values of the distances between neighbors, the amounts of neighboring pairs of magnetic atoms in the unit cell at each distance and the value of the PCC are also printed on screen. 将输入参数插入文件 input.py 后，您可以 1_coupling_constants.py启动脚本。它将计算耦合 常量，它将在屏幕上打印它们的值。 此外，它还将创建一个文件model.png，其中包含 海森堡模型能量与 DFT 能量的关系 Control 组。相邻要素之间的距离值、 晶胞中每个距离处的相邻磁原子对和 PCC 的值也打印在屏幕上。
We suggest to run the script 1_coupling_constants.py a couple of times with the append_coupling_constants flag set to False and with different values of the cutoff_radius, in order to investigate how many neighbors you need to include for obtaining a well-converged Heisenberg model. Once you are satisfied with the model’s accuracy, run the script a last time with the append_coupling_constants flag set to True and Automag will append the values of the coupling constants and the distances between neighbors to the file input.py. Now you are ready to run the script 2_write_vampire_ucf.py, which will read the file input.py and will produce a VAMPIRE unit cell file vamp.ucf. Now you can run VAMPIRE on your cluster using the unit cell file written by Automag, an input file and a material file specific for your problem. Automag contains a sample input file and a sample material file in the folder 3_monte_carlo/vampire_input, which can be simply edited and adapted to the problem under study. Once the VAMPIRE run is done, you can copy the output file in the folder 3_monte_carlo and run the last script 3_plot_results.py. It will produce a plot of the mean magnetization length (of the spin-up channel for antiferromagnetic materials) versus temperature obtained from the Monte Carlo simulation and it will fit the data using the analytical expression of the mean magnetization length, obtaining the values of the critical temperature and of the critical exponent. The fitted values of these two parameters are printed on screen. 1_coupling_constants.py我们建议使用 append_coupling_constants标志设置为 False，并且 cutoff_radius，为了调查您需要包含多少个邻居 获得收敛良好的 Heisenberg 模型。一旦您对 model 的 accuracy 时，请使用 append_coupling_constants flag 设置为 True 将附加耦合常量的值 以及 input.py 的邻居之间的距离。现在您已准备好 运行脚本 2_write_vampire_ucf.py，它将读取文件 input.py 将生成一个 VAMPIRE 晶胞文件 vamp.ucf。现在，您可以在 使用 Automag 编写的晶胞文件、输入文件和材质的簇 特定于您的问题的文件。Automag 包含一个样本输入文件和一个样本 material 文件，可以是 3_monte_carlo/vampire_input 文件夹中，也可以是简单的 编辑并改编了正在研究的问题。VAMPIRE 运行完成后，您 可以复制 3_monte_carlo 文件夹中output文件并运行最后一个脚本 3_plot_results.py。它将生成平均磁化长度 （的 反铁磁材料的自旋向上通道）与温度的关系 来自 Monte Carlo 模拟，它将使用分析来拟合数据 平均磁化长度的表达式，得到临界 温度和临界指数。这两个的拟合值 参数打印在屏幕上。
磁晶各向异性 magnetic anisotropic energy (MAE) 是指自旋方向在不同方向上的能量差。MAE跟体系对称性有关。 对于简单体二维系，易磁化轴可能在面内或者垂直于面，我们可以选择计算不同方向的某一磁性态能量对比来求MAE 比如CrI3， 我们可以计算沿着x方向的FM能量Ex，然后计算沿着z方向的FM能量Ez。 MAE = (Ex -Ez)/2 对于对称性较低的体系，需要每个面去计算MAE。在vasp中，三维体系研究MAE，可以在高对称轴上算能量差
b. Seeds文件夹，这个文件夹就是USPEX计算时，需要种子文件存放的地方。好了，既然需要种子文件，那么我们来准备种子文件吧。本教程要进行搜索的是 La-H 体系，那么开始准备La、H结构文件吧。首先到日本晶体结构数据NIMS网站(https://crystdb.nims.go.jp/crystdb/search-materials: 免费注册且免费下载，并且里面结构文件多数有其引用的相关文献)上找结构。
c. Specific文件夹,这个文件夹主要存放的是USPEX控制的VASP进行计算的输入文件，由于本教程是La-H体系在0 GPa的变组分结构搜索，主要是用VASP计算USPEX生成结构的能量，因而这个文件内主要是POSCAR进行结构优化的文件： INCAR_1-5: 这样设置的主要原因：考虑有你的初始结构通常远离局部最小点，在这种情况下，
图中绿色的点就是本次计算的结构的，而最下面的那条黑色就是凸包线，位于这条线上的点就是热力学上稳定的结构。在本图中从左到右总共有6个落在线上的点，和OUTPUT.txt里面的6个稳定的成分是相对应的。读者应该知道能够很清楚的知道X轴代表什么意思，怎么计算的。而Y轴形成焓（Enthalpy of formation，eV/atom）是由下面公式计算出来的： Y = (E(AxBy)-xE(A)-yE(B))/(x+y) 在案例12中，x: Mo原子的个数，y:B原子的个数。
案例13：USPEX能很容易地找到最无序的合金结构。这个案例只是为了演示TixCo(1-x)O。您需要在/Seeds/POSCARS中指定初始结构，并且只使用置换突变操作。（在这种情况下，不需要使用任何外部代码。在这个例子中，我们优化（最小化）结构有序度（Oganov and Valle（2009）；Lyakhov Oganov Valle（2010））而不需要进行结构优化（abinitioCode=0）。种子结构（Ti-Co-O结构的超级细胞）中的Ti原子和Co原子被置换，以寻找具有最小的有序度。在这种情况下最小化有序度，我们得到了“特殊准随机结构”的广义形式。）
ML_LOGFILE The log file for all MLFF-related details; training status, current errors and other important quantities can be extracted from here.(search err)
ML_ABN This file contains the collected training structures and a list of selected local reference configurations.
flag = R: interatomic distance between atoms item(1) and item(2).
flag = A: angle defined by atoms item(1), item(2) and item(3) (with the atom item(2) being the apex).
flag = T: torsional angle defined by atoms item(1), item(2), item(3) and item(4).
flag = M: distance between atom item(1) and the center of bond between atoms item(2) and item(3).
flag = B: distance between the center of bond between atoms item(1) and item(2) and the center of bond between atoms item(3) and item(4)
flag = P: ratio of length of the bond between atoms item(1) and item(2) and the length of the bond between atoms item(3) and item(4)
flag = W: function ${\displaystyle {\frac {1-\left(R/c\right)^{M}}{1-\left(R/c\right)^{N}}}}$ wit with 𝑅 being the bond length (in Å) between the atoms item(1) and item(2), 𝑐 is the reference bond length specified as item(3), and the exponents 𝑀 and 𝑁 are defined as item(4) and item(5), respectively.
flag = X, Y, and Z: fractional (direct) coordinates linked with the lattice vectors 𝑎, 𝑏, and 𝑐.
flag = cX, cY, and cZ: Cartesian coordinates 𝑥, 𝑦, and 𝑧.
flag = LR: length of lattice vector item(1)
flag = LA: angle between lattice vectors item(1) and item(2)
flag = LV: cell volume (no item(i) is defined in this case)
complex coordinates
flag = S: simple linear combination of primitive coordinates, i.e., $\left(\xi =\sum { {i=1} }^{ {M} }c{i},q_{i}\right)$.
flag = C: norm of the vector of primitive coordinates, which reads $\left(\xi ={\sqrt {\sum { {i=1} }^{ {M} },(c{i},q_{i})^{2}}}\right)$.
flag = IS: path-based coordinate[2] measuring progress along discretized path represented by ${\displaystyle {\tilde {q}}(j)}$ predefined in file IRCCAR, i.e., ${\displaystyle \xi ={\frac {1}{N-1}}{\frac {\sum {i=1}^{N}(i-1)\exp \left(-\sum {j=1}^{M}c{j}(q{j}-{\tilde {q}}{j}(i))^{2}\right)}{\sum {i=1}^{N}\exp \left(-\sum {j=1}^{M}c{j}(q{j}-{\tilde {q}}{j}(i))^{2}\right)}}}$
flag = IZ: path-based coordinate[2] measuring orthogonal distance from the path ${\displaystyle {\tilde {q}}}$ predefined in file IRCCAR, i.e., ${\displaystyle \xi =-{\frac {1}{c_{1}}}\log \sum _{i=1}^{N}\exp \left(-\sum {j=1}^{M}c{j}(q-{\tilde {q}}(i))^{2}\right)}$ with 𝑁 as defined above The complex coordinates are functions defined in the space spanned by the primitive coordinates.
这种格式的优势在于，可以为每个距离单独设置系数 M 和 N (请参阅上面 W 的定义)。此外，这种格式还允许简单直观地定义多个配位数的按比例总和和/或差异，只需通过与 S 相关联的系数的合适选择即可，这些系数可以是正、负或零。
例子3：约束体积或形状
固定体积变形状
1
LV 0
固定基矢角度变长度
1 2 3
LA 1 2 0 LA 1 3 0 LA 2 3 0
固定形状立方晶系变体积
Note that the S type constraints involving the lengths of the lattice vectors (𝑎𝑖−𝑎𝑗=0) are chosen so as to preserve ratios 𝑎1:𝑎2:𝑎3=1:1:1, as required by the cubic shape.
1 2 3 4 5 6 7 8 9 10 11 12
LA 1 2 0 LA 1 3 0 LA 2 3 0 LR 1 0 LR 2 0 LR 3 0 S 1 0 0 0 0 0 0 S 0 1 0 0 0 0 0 S 0 0 1 0 0 0 0 S 0 0 0 1 -1 0 0 S 0 0 0 1 0 -1 0 S 0 0 0 0 1 -1 0
固定单斜晶系变体积
Here, in order to fix ratios between the lengths of the lattice vectors (𝑎1:𝑎2:𝑎3), we define the constraints of the form 𝑐𝑖∗𝑎𝑖+𝑐𝑗∗𝑎𝑗=0. For instance, if the cell vectors are such that the relative proportions of their lengths are 𝑎1:𝑎2:𝑎3=1:1.5:2, the following ICONST can be used:
1 2 3 4 5 6 7 8 9 10 11 12
LA 1 2 0 LA 1 3 0 LA 2 3 0 LR 1 0 LR 2 0 LR 3 0 S 1 0 0 0 0 0 0 S 0 1 0 0 0 0 0 S 0 0 1 0 0 0 0 S 0 0 0 1.5 -1.0 0.0 0 S 0 0 0 2.0 0.0 -1.0 0 S 0 0 0 0.0 4.0 -3.0 0
从头算分子动力学（ab initio molecular dYnamics, AIMD），又称为第一性原理分子动力学，为研究电子与原子核相互耦合系统的力学进动过程的理论方法。AIMD是在早期经验力场分子动力学基础上发展起来的理论方法。经验力场分子动力学中，针对原子核之间相互作用，均采用经验性的力场，这种力场往往是对相互作用，基于这些相互作用势场好处是计算量小，能够模拟的体系较大，模拟时长也比较长，统计数据也较多。
Now reading vasp MD energies and temperature. Now reading vasp XDATCAR. Total frames 1000, NpT is False Finish reading XDATCAR. Selected time-range:0.0~499.5fs [debug] Now entering function plotfigure.....
Extensions-Analysis-Radial Pair Distribution Function g(r)，选择分析H(type H)在O(type O)周围的概率分布。值得注意的是分析RDF时,横坐标也就是max r不能超过盒子最小边长的一半，也就是得满足最小映像约定。如图4所示，在计算RDF时，如果max r的取值大于盒子最小边长的一半，就有可能重复算到一个粒子和它的映像粒子，这使得程序的周期性判断失准。将生成的dat文件的第一列和第二列作图即可得到RDF图。
使用以下的脚本可以分别得到所有原子，氢原子，氧原子的RMSD，并调用matplotlib绘制RMSD图。网页中有一段话（Note If you use trajectory data from simulations performed under periodic boundary conditions then you must make your molecules whole before performing RMSD calculations so that the centers of mass of the selected and reference structure are properly superimposed.）也就是在计算RMSD的时候选择的分子必须是完整的，不能分处于盒子的两边。这与我们之前的描述是一致的。MD Analysis默认对齐了分子。
u = MDAnalysis.Universe('XDATCAR.pdb', permissive=True) ref = MDAnalysis.Universe('XDATCAR.pdb', permissive=True) # reference (with the default ref_frame=0) ref.trajectory[0] #use first frame as reference R = MDAnalysis.analysis.rms.RMSD(u, ref, select="all", # superimpose on whole backbone of all atoms # align based on all atoms groupselections=["type H","type O"], filename="rmsd_all.dat",center=True)#, # CORE timestep=0.0005 #0.5fs from fs to ps asReader has no dt information, set to 1.0 ps R.run() rmsd = R.rmsd.T # transpose makes it easier for plotting time = rmsd[1]*timestep
import os from pymatgen.core.trajectoryimportTrajectory from pymatgen.io.vasp.outputsimportXdatcar from pymatgen importStructure from pymatgen.analysis.diffusion_analyzerimportDiffusionAnalyzer import numpy as np import pickle
C = diff.conductivity D = diff.diffusivity withopen('result.dat','w') asf: f.write('# AIMD result for Li-ion\n') f.write('temp\conductivity\diffusivity\n') f.write('%d\t%.2f %.10f' %(300,C,D))
Pymatgen 还提供了很多数据库（Materials Project REST API，Crystallography Open Database，and other external data sources）的接口，方便大家从数据库中查询结构和其他数据。
以下是Pymatgen官网提供的后处理的例子：
Top: (left) Phase and (right) Pourbaix diagram from the Materials API. Bottom left: Calculated bandstructure plot using pymatgen’s parsing and plotting utilities. Bottom right: Arrhenius plot using pymatgen’s DiffusionAnalyzer.
[test@ln0%tianhe2 li_sn_s]$ python Python 3.7.3 (default, Mar 27 2019, 22:11:17) [GCC 7.3.0] :: Anaconda, Inc. on linux Type "help", "copyright", "credits" or "license" for more information. >>> import pymatgen >>>
@classmethod def from_structures(cls, structures, specie, temperature, time_step, step_skip, initial_disp=None, initial_structure=None, **kwargs): """ Convenient constructor that takes in a list of Structure objects to perform diffusion analysis. Args: structures ([Structure]): list of Structure objects (must be ordered in sequence of run). E.g., you may have performed sequential VASP runs to obtain sufficient statistics. ... ... """
import os from pymatgen.core.trajectory import Trajectory from pymatgen.io.vasp.outputs import Xdatcar from pymatgen import Structure from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer import numpy as np import pickle
在“Design principles for solid-state lithium superionic conductors”一文中（Wang et al., Nature Materials 2015, 14 , 1026–1031. ），作者用Ab Initio Molecular Dynamic （AIMD）计算了Li 离子在Li$\mathrm{_1}\mathrm{_2}， \mathrm{_7} \mathrm{_3} \mathrm{_1}$$\mathrm{_1}， \mathrm{_2}，和 \mathrm{_4} \mathrm{_4}$ 四种材料中的迁移概率密度（Probability Density），结果如下图所示：
从图中可以看出Li离子在图a所示材料中主要沿c轴方向的通道迁移，而且由于这个通道连通得比较好，Li离子的迁移势垒会比较低（0.220.25 eV）。 Li离子在图b所示的材料的迁移路径形成了一个三维网格，而且由于这个概率密度比图b中的概率密度分布得更加均匀，Li离子的迁移势垒更低（0.180.19 eV）。 图b所示的材料就完全不行了，因为Li离子的概率密度仅分布在特定的位点附近，说明离子不能有效地移动。 Li离子在图d所示材料中也存在迁移局域化的行为。 作者总结说 “A general principle for the design of Li-ion conductors with low activation energy can be distilled from the above findings: all of the sites within the diffusion network should be energetically close to equivalent, with large channels connecting them.” 那么我们如何在自己的计算中画出这样的图呢？Pymatgen 举手说，它可以帮忙！
我们把其中的 pymatgen_diffusion 文件夹放到 Anaconda的site-packages文件夹下，路径是 Windows 系统：……\Anaconda\Lib\site-packages；Linux系统：……/anaconda3/lib/pythonx.x/site-packages，就算安装好了。
接下来我们可以启动python，导入这个模块，如果不报错就没有问题了。
1 2 3 4 5 6
[test@ln0%tianhe2 li_sn_s]$ python Python 3.8.3 (default, Jul 2 2020, 16:21:59) [GCC 7.3.0] :: Anaconda, Inc. on linux Type "help", "copyright", "credits" or "license" for more information. >>> import pymatgen_diffusion >>>
学习用法 我们可以在其github网站上通过例子学习这个模块的用法。
点击打开 probbility_analysis.ipynb 文件。
其内容如下（有所删减）：如果不想看的话可直接查看 开始作图 部分
1 2 3 4 5 6 7 8 9 10 11 12 13 14
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis import json
data = json.load(open("../pymatgen_diffusion/aimd/tests/cNa3PS4_pda.json", "r")) diff_analyzer = DiffusionAnalyzer.from_dict(data) # 初始化DiffusionAnalyzer类
pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diff_analyzer, interval=0.5, species=("Na", "Li")) #可以指定离子 #Save probability distribution to a CHGCAR-like file pda.to_chgcar(filename="CHGCAR_new2.vasp") #保存概率密度文件
开始作图 代码（test.py）如下：
1 2 3 4 5 6 7 8 9
from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis from pymatgen.core.trajectory import Trajectory from pymatgen.io.vasp.outputs import Xdatcar from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
STATIC DIELECTRIC TENSOR (including local field effects in DFT)``` 电子贡献或者
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
```MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))``` ```MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION``` 离子贡献 ```BORN EFFECTIVE CHARGES (in e, cummulative output)``` ```ELASTIC MODULI IONIC CONTR (kBar)``` ```PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)``` ### 频率依赖介电函数：读取WAVECAR IPA： ```javascript {.line-numbers} ALGO = Exact NBANDS = 64 LOPTICS = .TRUE. CSHIFT = 0.100 NEDOS = 2000 #and you might try with the following #LPEAD = .TRUE.
frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects) and frequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)
上面计算完介电常数后，一般正交晶系。 grep "PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)" OUTCAR -A10 & grep "PIEZOELECTRIC TENSOR for field in x, y, z (C/m^2)" OUTCAR -A10
$E = 6062.434883706029P^6 + 2437.756598493882P^4 + -340.5806845162749*P^2 + 10.339600058315137$ 其中参数 A, B, C 分别为
$A = -0.68116 eV *(m^2/C)^2 B = 9.75103 eV *(m^2/C)^4 C = 36.37461 eV *(m^2/C)^6$
绘制能量曲线的脚本
{.line-numbers}
1 2 3 4 5 6 7 8 9 10 11
withopen('energy.dat','r') asf: content = f.readlines() res = [float(i.strip('\n')) for i in content] minres = min(res) res = [1e3*(i-minres) for i in res] import matplotlib.pyplotas plt plt.figure() plt.plot(res,'b.-') plt.xlabel('displacement') plt.ylabel('Energy ($meV$)') plt.show()
### 考虑极化量子 c = 4.0335000000000001 tmp = [] N = 20 for j inrange(N): start = -N/2+j tmp1 = [i+start*c for i in res] # print(tmp1[0]) tmp.append(tmp1) plt.figure() for i intmp: plt.plot(i,'.')
#!/bin/bash path=`pwd` #记录当前文件夹位置并定义到变量 path中，称为初始目录 pwd for cij in`ls -F | grep /$` #ls -F| grep /$ 命令为返回当前文件夹中所有子目录 #将所有子目录添加到变量cij中，并准备进行for循环 do cd ${path}/$cij #进入cij变量所代表的子目录中 pwd for s in strain_* #将cij代表的子目录中所有strain文件夹添加到变量s中 #此部添加变量s的方式可更改为添加变量cij的方式，但不建议这么做 do cd ${path}/$cij/$s #进入初始目录下变量cij和变量s组成的目录，由于for循环的设置，将遍历所有可计算的strain文件夹 pwd mpirun -np 16 vasp_std #执行vasp计算 cd ${path} #返回初始目录 pwd done done
>dfc2 DFC2--a generator of renormalized harmonic FCs from SCPH outputs. XML file containing original Fc2 : STo222.xml Output xml filename with anharmonic correction : ST0222_SCPH2-2_300K.xml FC2 correction file fromSCPH calculation :STo_scph2-2.scph_dfc2 Target temperature : 300
New XML file STo222_SCPH2-2_300K.xml was created successfully.
&cell ## primitive cell is best 6.2507550151332 -0.5 0.5 0.5 0.5 -0.5 0.5 0.5 0.5 -0.5 /
&kpoint 1 # KPMODE = 1: get phonon band 2: get others properties G 0.0 0.0 0.0 H 0.5 -0.5 0.5 201 H 0.5 -0.5 0.5 P 0.25 0.25 0.25 201 P 0.25 0.25 0.25 G 0.0 0.0 0.0 201 G 0.0 0.0 0.0 N 0.0 0.0 0.5 201 /
[scf1355@ln151%nc-l ah]$ ~/soft/alamode/tools/dfc2 DFC2 -- a generator of renormalized harmonic FCs from SCPH outputs. XML file containing original FC2 : crte2-442.xml Output xml filename with anharmonic correction : crte2-new.xml FC2 correction file from SCPH calculation : crte2-442.scph_dfc2 Target temperature : 300
The accuracy of the anharmonic force constants can be checked by increasing the size of the training dataset. If the physical properties, such as thermal conductivity and free energy, change significantly with the size of the training dataset, it indicates that the training is insufficient. If the change is minor, your force constants may be OK.
报错：Warning: The following force constant doesn’t exist in the original file: You can neglect this warning as long as the correction term (the last column) is small. This warning is raised when the absolute value of the correction is larger than 1.0e-10 and the corresponding force constant does not exist in the original FCSXML file. In your case, the correction term is only slightly above the tolerance. So, it should not cause any problems hopefully. When the correction term is much larger than 1.0e-10, this warning indicates a possible error in the inputs. Likely errors are the input structure is not fully symmetrized, or the numerical digits of lattice vectors and fractional coordinates are insufficient. the given KMESH_INTERPOLATE is not commensurate with the supercell size of the original FCSXML.
In the current implementation of alamode, only the former choice, renormalized 2FCs + conventional (old) 3FCs, is supported. Another possible choice may be renormalized 2FCs + renormalized 3FCs, but it is not supported within alamode yet. When MODE = RTA, the anphon code only reads the 2FCs and 3FCs from given XML files. By contrast, when MODE = SCPH, the code additionally reads the 4FCs because they are necessary to compute the renormalized 2FCs based on the self-consistent phonon theory. Therefore, the second combination of renormalized 2FCs + high order FCs using SCPH is not possible.
I performed the SCPH calculation. It is not converged at low temperatures (0K, 100K, 200K), while it achieved convergence at high temperatures (300~1200K).A potential solution is to use a smaller MIXALPHA and a smaller DT.
I would like to use second-order, third-order, and fourth-order force constants to obtain anharmonic potentials and calculate the potential energy as a function of atomic displacement. My goal is to generate a plot similar to the one in the paper “Orbitally driven giant phonon anharmonicity in SnSe” (Nature Phys. 2015).Please use https://github.com/ttadano/alamode/blob/develop/tools/calcpes.py.
now that I have obtained the renormalized force constant file xxx.xml at a finite temperature, I now want to convert the renormalized force constant file directly to FORCE_CONSTANS
How to convert a second-order force constant matrix obtained by almode to FORCE_CONSTANTS in phonopy format？
python convert.py alamode.xml > FORCE_CONSTANTS
Could you give me some suggestion? How shoud I covert this with phonopy?
FULL_FORCE_CONSTANTS = .TRUE. WRITE_FORCE_CONSTANTS = .TRUE. Note: the old FORCE_CONSTANTS will be overwrited
How small the fitting error should be?
It depends on the Taylor expansion potential and displacement magnitude you choose. In the standard harmonic calculation where –mag=0.01 is used in displace.py and the all harmonic interactions are considered (cutoff = None), the fitting error is usually less than 5% (~1–2% in most cases). In the calculation of third-order force constants with –mag=0.04, the fitting error should be small as well. Indeed, in many cases, we obtain much smaller fitting errors (< 1%) than the harmonic case. In the temperature-dependent effective potential method, we try to fit the harmonic potential to the displacement-force datasets sampled by ab initio molecular dynamics at finite temperature. Therefore, the fitting error tends to be much larger (> 10%).
h_ref=`grep " $homo_ref " $j | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'` l_ref=`grep " $lumo_ref " $j | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'` echo "The band gap (in eV) without zero-point vibrations is:" echo "$h_ref $l_ref" |awk '{print ($2-$1)}' echo "The band gap (in eV) including zero-point vibrations is:" echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}' echo "The zero-point renormalization of the band gap (in eV) is:" echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $h_ref $l_ref" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-($11-$10)}'
输出如下：
1 2 3 4 5 6
The band gap (in eV) without zero-point vibrations is: 4.4049 The band gap (in eV) including zero-point vibrations is: 4.05102 The zero-point renormalization of the band gap (in eV) is: -0.353883
bandshift=0 gwrun=-1 dgbd=-1 val=-1 con=-1 gwldadiff=-1 test=-1 while [[ $# -gt 0 ]] do key="$1" case $key in esac shift done
if [ -f "helpscript.perl" ]; then rm helpscript.perl fi
cat > helpscript.perl <<EOF #!/bin/perl
use strict; use warnings;
my \$zahler=0; my @entropy=0; my \$ezp=0; my \$fhelmholtz=0; my \$uenergy; my \$kboltzmann=8.6173303*10**(-5.0); my \$ntemp=8; my \$tmax=700; my \$tmin=0; my \$tstep=100; for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { \$entropy[\$itemp]=0; } while (<>) { chomp; \$_=~s/^/ /; my @help=split(/[\t,\s]+/); \$zahler=\$zahler+1; if (\$zahler == 1) {\$uenergy=\$help[1];} else { my \$homega=\$help[2]/1000; \$ezp=\$ezp+\$homega*0.5; for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { my \$temp=(\$itemp-1)*\$tstep+\$tmin; my \$kbt=\$kboltzmann*\$temp; if (\$temp < 0.0000001) { } else { \$entropy[\$itemp]=\$entropy[\$itemp]-\$kboltzmann*log(1-exp(-\$homega/\$kbt)); \$entropy[\$itemp]=\$entropy[\$itemp]+\$kboltzmann*\$homega/\$kbt*(1/(exp(\$homega/\$kbt)-1)); } } } last if eof; }
\$ezp=\$ezp; #/ \$zahler;
for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++) { my \$temp=(\$itemp-1)*\$tstep+\$tmin; \$entropy[\$itemp]=\$entropy[\$itemp]; # / \$zahler; \$fhelmholtz=\$uenergy+\$ezp-\$temp*\$entropy[\$itemp]; printf ("%15.8e %15.8e\n",\$temp,\$fhelmholtz); } printf ("#NOTEMP %15.8e\n",\$uenergy); EOF
，属于四方晶系。
