计算化学公社

标题: 使用ElaStic和Quantum Espresso计算弹性常数与模量 [打印本页]
作者
Author: biogon    时间: 2021-5-24 17:46
标题: 使用ElaStic和Quantum Espresso计算弹性常数与模量
本帖最后由 biogon 于 2022-7-14 13:41 编辑

      关于如何使用ElaStic和exciting计算弹性常数和模量，在http://exciting.wikidot.com/carbon-ElaStic已经有了教程，但是此教程在使用Quantum Espresso计算弹性模量时很多东西是不适用的，故因此写本教程。
1.   计算前准备（默认已经安装好Quantum Espresso）
1.1 ElaStic版本与依赖库的选择
ElaStic最初的版本是基于python2.*环境下运行的，最新版本为ElaStic 1.1，后来有人基于python 3.*环境重写了该程序，称为ElaStic 2020 其可在python 3.* 环境下运行。二者均需要numpy才可运行。已安装anaconda python 3的使用者可以同时安装python 2.* 版本的运行环境，使用时激活即可，在此不再赘述。
1.2 ElaStic 的安装
ElaStic 1.1的下载地址：http://exciting.wdfiles.com/loca ... /ElaStic_1.1.tar.gz
ElaStic 2020 的github： https://github.com/ponychen123/Elastic2020
ElaStic 1.1的安装相当简单，下载压缩包，解压后，设置如下环境变量
  1. export ElaSticROOT=/home/$USER/software/ElaStic_1.1
  2. export PATH=/home/$USER /software/ElaStic_1.1:$PATH
复制代码
环境变量根据自己情况可自行修改
ElaStic2020与此基本相同，不再赘述
1.3 garce 的安装
ElaStic1.1的可视化分析所需xmgrace，这个包一般系统是不自带的，且官方网站目前是无法下载到的，但是可以利用系统的包管理器一键安装。
  1. $ sudo apt-get install grace  #Ubuntu
  2. $ sudo yum install grace  #cent os
  3. $ sudo zypper install xmgrace  #sles
复制代码
即可成功安装，如无法安装请更新源。
1.4  Sgroup
Sgroup是判断空间群的小工具，其在ElaStic文件夹中有一压缩包adon_v1_0.tar.gz，执行以下命令。
  1. $ tar -zxvf adon_v1_0.tar.gz
  2. $ cd SpaceGroups
  3. $ make
  4. $ cp sgroup ..
复制代码
即可
至此所有的准备都已经完成，即可开始准备计算。
2.   开始计算
2.1 输入文件的准备
由于我们是使用QuantumEspresso和ElaStic计算，所以计算前需要新建文件夹叫***_2nd，在里面准备***.in这样的pw.x输入文件，内容如下：
  1. &CONTROL
  2.    calculation      = 'relax'
  3.    prefix           = '****'
  4.    outdir           = './'
  5.    pseudo_dir       = '/……/pp'
  6.    tprnfor          = .true.
  7.    forc_conv_thr    = 1d-4
  8.    tstress          = .true.
  9. /
  10. &SYSTEM
  11.    ibrav            = 0
  12.    nat              = 2
  13.    ntyp             = 2
  14.    ecutwfc          = 40
  15.    ecutrho          = 200
  16.    occupations      = 'smearing'
  17.    smearing         = 'gaussian'
  18.    degauss          = 0.02
  19.    celldm(1)         = ****
  20. &ELECTRONS
  21.    conv_thr         = 1d-8
  22. /
  23. &IONS
  24.    ion_dynamics     = 'bfgs'
  25. /
  26. ATOMIC_SPECIES
  27.    ** 1.0 **_ **.uspp.F.UPF
  28.    **  1.0 **_**.uspp.F.UPF
  29. ATOMIC_POSITIONS (crystal)
  30. **            0.0000000000        0.0000000000       -0.0000000000
  31. **            0.2500000000        0.2500000000        0.2500000000
  32. K_POINTS automatic
  33. ** ** ** ** ** **
  34. CELL_PARAMETERS (alat)
  35. 0.707107   0.707107   0.000000
  36. 0.000000   0.707107   0.707107
  37. 0.707107   0.000000   0.707107
复制代码
需要注意一下，晶格常数必须使用alat表示的作为输入，celldm(1)是a.u.，否则最后结果是完全错误的。其他的内容根据自己的实际修改输入文件即可。
2.2 生成计算弹性常数的输入文件
使用ElaStic 1.1和Quantum Espresso的计算流程如下：
  1. $ cd /$USER/……/***_2nd
  2. $ ElaStic_Setup

  4.      Which DFT code would you like to apply for the calculations?
  5.      exciting  ---------=>  1                                    
  6.      WIEN2k    ---------=>  2                                    
  7.      Quantum ESPRESSO --=>  3

  9. >>>> Please choose (1, 2, or 3): 3

  11.      Energy  ---=>  1           
  12.      Stress  ---=>  2   
  13. >>>> Please choose the method of the calculation (choose 1 or 2): 1

  15.      2nd  ---=>  2           
  16.      3rd  ---=>  3   
  17. >>>> Please choose the order of the ElaStic constant (choose 2 or 3): 2

  19. >>>> Please enter the name of the Quantum-ESPRESSO input file: ***.in         

  21.      Number and name of space group: 227 (F d -3 m) [origin choice 2]      
  22.      Cubic I structure in the Laue classification.      
  23.      This structure has 3 independent second-order ElaStic constants.

  25. >>>> Please enter the maximum Lagrangian strain
  26.      The suggested value is between 0.030 and 0.150: 0.05
  27.      The maximum Lagrangian strain is 0.09

  29. >>>> Please enter the number of the distorted structures [odd number > 4]: 15
  30.      The number of the distorted structures is 15

  32. $
复制代码
1.    3 选择你所用的计算程序
2.    1 选择使用能量还是应力计算弹性常数
3.    2 选择要计算感兴趣的弹性常数的阶数
4.    ***.in，你的原始输入文件的名称
5.    0.09，我们要对其进行计算的最大应变的绝对值
6.    15 在最大负应变和最大正应变之间产生的，等距分布的变形结构的数量。
最后该脚本会生成一些新目录和文件，可以使用ls命令查看。
ElaStic 2020 的使用流程基本与此相同，除了一开始不用选择计算程序，如果你只装了Quantum Espresso。
3.   进行计算
在原始输入文件***.in所在目录***_2nd下新建一个shell脚本命名为诸如elastic_pwx_run.sh之类的名字，写入如下内容
  1. #!/bin/bash
  2. label=`ls -d Dst??`
  3. for Dstn in $label ; do
  4.     cd $Dstn
  5.     Dstn_num_list=`ls -d ${Dstn}_??`
  6.     for Dstn_num in $Dstn_num_list ; do
  7.         cd $Dstn_num/
  8.           for inf in *.in
  9.           do
  10.           echo Running ${inf} ...
  11.           time  pw.x < ${inf} > ${inf//in/out}
  12.           echo ${inf} has finished
  13.           echo
  14.           done
  15.         cd ../
  16.     done
  17.     cd ../
  18. done
复制代码
如在超算等等特殊环境下运行，自行在#!/bin/bash后添加所需内容即可。
chmod+x对其赋予可执行权限，运行之，等待其运行完成即可进行分析。
4.   数据分析
本流程需要在有GUI机器或者GUI机器进行X11远程转发才可以进行，否则需要手动作图分析。
这里以http://exciting.wikidot.com/carbon-elastic的图作为示例进行讲解。
计算弹性常数的方法在Computer Physics Communications 184 (2013) 1861–1873 一文中有详细描述，在此不再复述。
执行
  1. $ ElaStic_Analyze
复制代码
此时，屏幕上将自动出现xmgrace图
(, 下载次数 Times of downloads: 86) (, 下载次数 Times of downloads: 95) (, 下载次数 Times of downloads: 110)
三图分别为Dst01、Dst02、Dst03的d2E-η 关系图
如果发现任何曲线都无收敛趋势，则应继续增大变形结构的数目，否则最后的结果误差会较大。
通过分析Dst01的第一个图，我们注意到4阶多项式的曲线显示出在4060-4080 GPa处收敛的趋势。从拟合的角度来看，可以将其假定为二阶导数的收敛值。对于这种变形类型，此值等于体积模量的9倍。因此，体积模量的值为451-453GPa。
ElaStic_Analyze将生成文件ElaStic_2nd.in，在接下来的部分中用于计算模量用。
如何手动读取数据分析？运行完ElaStic_Analyze会在/....../***_2nd目录下生成Energy-vs-Strain文件夹一个，里面有Dst**_d2e.dat，这个就是上述图的数据的原始文件，可以提取里面的数据自行作图或直接肉眼观察是否收敛
5.   数据后处理
为了获得二阶弹性常数的数值，需要进行如下步骤进行分析。
第一步是编辑文件ElaStic_2nd.in，其格式为：
  1. Dst01     eta_max     Fit_order
  2. Dst02     eta_max     Fit_order
  3. Dst03     eta_max     Fit_order
复制代码
在此文件的每一行中，应为eta_max和Fit_order输入一个数值。根据前面各图的可视化分析结果，通过观察相应图的曲线收敛区域可得到这两个值：eta_max是位于曲线的收敛区域中的最大变形值，Fit_order代表收敛曲线的阶次。从上面可知，我们可以将ElaStic_2nd中的参数设置为如下的数值：
  1. Dst01     0.05     4
  2. Dst02     0.05     4
  3. Dst03     0.05     6
复制代码
现在，为了计算弹性常数和模量，执行ElaStic_Result
  1. $ ElaStic_Result
复制代码
弹性常数将输出到文件ElaStic_2nd.out中。从中可以读取得到弹性常数和模量的数值。
  1.     The output of ElaStic code                                             
  2.     Today is Mon May 10 10:38:45 2021

  4.     Symmetry of the second-order elastic constant matrix in Voigt notation.
  5.     for, space group-number between 207 and 230, Cubic I structure.        

  7.                C11     C12     C12      0       0       0                  
  8.                C12     C11     C12      0       0       0                  
  9.                C12     C12     C11      0       0       0                  
  10.                 0       0       0      C44      0       0                  
  11.                 0       0       0       0      C44      0                  
  12.                 0       0       0       0       0      C44                 

  14.     Elastic constant (stiffness) matrix in GPa:                             

  16.         -0.0         0.0         0.0         0.0         0.0         0.0
  17.          0.0        -0.0         0.0         0.0         0.0         0.0
  18.          0.0         0.0        -0.0         0.0         0.0         0.0
  19.          0.0         0.0         0.0         0.0         0.0         0.0
  20.          0.0         0.0         0.0         0.0         0.0         0.0
  21.          0.0         0.0         0.0         0.0         0.0         0.0


  24.     Elastic compliance matrix in 1/GPa:

  26.    -42.97254    -5.76901    -5.76901    -0.00000    -0.00000    -0.00000
  27.     -5.76901   -42.97254    -5.76901    -0.00000    -0.00000    -0.00000
  28.     -5.76901    -5.76901   -42.97254    -0.00000    -0.00000    -0.00000
  29.      0.00000     0.00000     0.00000   204.08372     0.00000     0.00000
  30.      0.00000     0.00000     0.00000     0.00000   204.08372     0.00000
  31.      0.00000     0.00000     0.00000     0.00000     0.00000   204.08372

  33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  35.     Voigt bulk  modulus, B_V =    -0.01  GPa
  36.     Voigt shear modulus, G_V =    -0.00  GPa

  38.     Reuss bulk  modulus, B_R =    -0.01  GPa
  39.     Reuss shear modulus, G_R =     0.01  GPa

  41.     Hill bulk  modulus,  B_H =    -0.01  GPa
  42.     Hill shear modulus,  G_H =     0.00  GPa

  44. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  46.     Voigt Young modulus,  E_V =    -0.01  GPa
  47.     Voigt Poisson ratio, nu_V =     0.32

  49.     Reuss Young modulus,  E_R =     0.08  GPa
  50.     Reuss Poisson ratio, nu_R =     2.64

  52.     Hill Young modulus,   E_H =     0.02  GPa
  53.     Hill Poisson ratio,  nu_H =     0.94

  55. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  57.     Elastic Anisotropy in polycrystalline, AVR = -158.322 %

  59. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  61.     Eigenvalues of elastic constant (stiffness) matrix:   

  63.             -0.0
  64.             -0.0
  65.             -0.0
  66.              0.0
  67.              0.0
  68.              0.0

  70.     ... Have a G00D Day, Week, Month, Year, and Century (if you are lucky) ...   
  71.                Bye-Bye! Tschuess! Ciao! Poka! Zia Jian! KhodaHafez!
复制代码

如果得到的弹性常数和模量有负值存在，这说明计算可能出现了问题，应检查原始输入文件中的***.in中的各项参数是否正确。
大家有什么问题以及建议都可以在下面回复。


作者
Author: yuxu1994    时间: 2021-5-25 13:27
@kyuu 来来来
作者
Author: ZZG    时间: 2021-10-16 13:53
请问下，我按照教程，安装成功之后。执行ElaStic_Setup 脚本发现生成的输入文件中出现重复的原子坐标和k点，这是为什么呢？我删除软件又装还是这样。
作者
Author: biogon    时间: 2021-10-16 20:59
ZZG 发表于 2021-10-16 13:53
请问下，我按照教程，安装成功之后。执行ElaStic_Setup 脚本发现生成的输入文件中出现重复的原子坐标和k点 ...

原始的输入文件咋写的
作者
Author: ZZG    时间: 2021-10-17 17:28
这是我的原始输入文件。现在是celldm（1）出现重复，其他的好了、
作者
Author: ZZG    时间: 2021-10-17 17:29
biogon 发表于 2021-10-16 20:59
原始的输入文件咋写的

您好，我的文件上传了，请看一下
作者
Author: ZZG    时间: 2021-10-20 20:43
请问下，我在执行ElaStic_Analyze的时候一直出现这种报错是为什么呢？
作者
Author: biogon    时间: 2021-10-21 08:44
ZZG 发表于 2021-10-20 20:43
请问下，我在执行ElaStic_Analyze的时候一直出现这种报错是为什么呢？

你用的是基于python3的elastic2020么？
作者
Author: ZZG    时间: 2021-10-21 09:26
biogon 发表于 2021-10-21 08:44
你用的是基于python3的elastic2020么？

是的，问了一些人，有的人说下载的ElaStic脚本有问题，我就是在楼主的链接里下载的
作者
Author: biogon    时间: 2021-10-21 10:53
ZZG 发表于 2021-10-21 09:26
是的，问了一些人，有的人说下载的ElaStic脚本有问题，我就是在楼主的链接里下载的

这个玩意我后来发现有bug，还没空检查bug在哪，改用基于python2的elastic1.1就没有问题
作者
Author: ZZG    时间: 2021-10-21 15:22
biogon 发表于 2021-10-21 10:53
这个玩意我后来发现有bug，还没空检查bug在哪，改用基于python2的elastic1.1就没有问题

嗯嗯，我又在重新看作者写的文章，看能不能用elastic1.1.请问下，在使用之前是不是要进行这样的配置环境呀？
作者
Author: ZZG    时间: 2021-10-21 15:23
就是这样路径的添加
作者
Author: biogon    时间: 2021-10-21 15:28
ZZG 发表于 2021-10-21 15:23
就是这样路径的添加

你仔细看下我写的就能明白elastic2020和1.1的配置几乎是无差别的，这个是exciting程序的环境变量的配置，不是elastic1.1的
作者
Author: ZZG    时间: 2021-10-25 18:56
请问下，我设置扭曲程度到最大0.15，为什么之后的Dst-3最后不收敛呢?我是应该去前面的最大扭曲吗？还是Dst-3扭曲方式就无法最后收敛呢？我设置0.05扭曲最后基本都收敛了。
作者
Author: biogon    时间: 2021-10-26 09:37
ZZG 发表于 2021-10-25 18:56
请问下，我设置扭曲程度到最大0.15，为什么之后的Dst-3最后不收敛呢?我是应该去前面的最大扭曲吗？还是Dst- ...

感觉变形量设的太高了，一般0.09就差不多了，如果你不想重算的话就Dst03把2视作是收敛的，看结果误差如何
作者
Author: ZZG    时间: 2021-10-26 21:01
biogon 发表于 2021-10-26 09:37
感觉变形量设的太高了，一般0.09就差不多了，如果你不想重算的话就Dst03把2视作是收敛的，看结果误差如何

误差有些大，20 GPa左右。请问这个扭曲程度的大小对结果有什么影响吗？我看变化挺大，或者对计算上有什么意义吗？
作者
Author: biogon    时间: 2021-10-27 09:15
ZZG 发表于 2021-10-26 21:01
误差有些大，20 GPa左右。请问这个扭曲程度的大小对结果有什么影响吗？我看变化挺大，或者对计算上有什么 ...

变形量设大了变形结构不够多的话会放大误差
作者
Author: ZZG    时间: 2021-10-27 11:54
biogon 发表于 2021-10-27 09:15
变形量设大了变形结构不够多的话会放大误差

变形结构我设到了41，是否是扭曲程度和添加的bianxing结构足够大，结果越精确呢？扭曲程度上限是0.15，我设置51和变形结构是否足够呢？
作者
Author: biogon    时间: 2021-10-28 09:12
ZZG 发表于 2021-10-27 11:54
变形结构我设到了41，是否是扭曲程度和添加的bianxing结构足够大，结果越精确呢？扭曲程度上限是0.15，我 ...

理论上是越多越准的
作者
Author: ZZG    时间: 2021-10-31 09:31
biogon 发表于 2021-10-28 09:12
理论上是越多越准的

嗯嗯  我计算一个16原子的体系，设置0.05扭曲 21个结构。算Dst01的所有文件将近一天。用了72个核数。这个合理吗？而且有的文件计算过程出现了stop，但是job done了
作者
Author: biogon    时间: 2021-10-31 11:12
ZZG 发表于 2021-10-31 09:31
嗯嗯  我计算一个16原子的体系，设置0.05扭曲 21个结构。算Dst01的所有文件将近一天。用了72个核数。这个 ...

那些非平衡结构一般都是要迭代很久的，变形量越大耗时越久
stop的你要看计算过程和输出的能量有没有问题
作者
Author: ZZG    时间: 2021-11-1 17:27
biogon 发表于 2021-10-31 11:12
那些非平衡结构一般都是要迭代很久的，变形量越大耗时越久
stop的你要看计算过程和输出的能量有没有问题

好的，别人告诉我，Smearing错了，金属不能用    occupations = 'smearing', smearing = 'gauss', degauss = 1.0d-5，半导体的可以。是这样吗？
作者
Author: mxh    时间: 2021-12-21 18:22
请问大家，我执行elastic_pwx_run.sh出现了图片上的错误，什么原因导致，怎样去解决
作者
Author: biogon    时间: 2021-12-21 20:15
ZZG 发表于 2021-11-1 17:27
好的，别人告诉我，Smearing错了，金属不能用    occupations = 'smearing', smearing = 'gauss', degaus ...

应该可以
作者
Author: biogon    时间: 2021-12-21 20:15
mxh 发表于 2021-12-21 18:22
请问大家，我执行elastic_pwx_run.sh出现了图片上的错误，什么原因导致，怎样去解决

检查输入文件，看是不是有问题
作者
Author: mxh    时间: 2021-12-23 10:43
请问大家，ElaStic和thermo_pw程序哪个能算二维材料的弹性常数。如果都不能的话，对接QE的哪个程序可以算二维材料的弹性常数
作者
Author: mxh    时间: 2021-12-23 10:51
请问大家，ElaStic和thermo_pw程序哪个能算二维材料的弹性常数。如果都不能的话，对接QE的哪个程序可以算二维材料的弹性常数
作者
Author: biogon    时间: 2021-12-23 12:00
mxh 发表于 2021-12-23 10:51
请问大家，ElaStic和thermo_pw程序哪个能算二维材料的弹性常数。如果都不能的话，对接QE的哪个程序可以算二 ...

我用ElaStic算过BN，感觉结果还行
作者
Author: mxh    时间: 2021-12-23 20:42
biogon 发表于 2021-12-23 12:00
我用ElaStic算过BN，感觉结果还行

请问你算BN的时候，有没有遇到这样的问题。原始输入文件写的alat，但软件生成的那一些输入文件变成了cubic，这样对吗
作者
Author: mxh    时间: 2021-12-23 20:46
biogon 发表于 2021-12-23 12:00
我用ElaStic算过BN，感觉结果还行

我还有一个问题。用ElaStic算二维材料，它本身用来算三维的，我们现在用它算二维，那材料的真空层设多大比较合适
作者
Author: hdb    时间: 2023-5-10 11:32
本帖最后由 hdb 于 2023-5-10 12:09 编辑

Python提供了版本2to3的工具，一般叫做2to3.py，或者2to3-script.py，用以下命令转换一下2版本的脚本即可，上述的2020版确实有bug，会导致重写已有card，自己转换下更合适。下面的命令2选1：
  1. 2to3 -w filename  # 需要用which 2to3确定PATH中有这个命令
复制代码
  1. python 2to3.py -w filename  #需要找到py文件的绝对路径
复制代码


作者
Author: biogon    时间: 2023-5-10 21:29
hdb 发表于 2023-5-10 11:32
Python提供了版本2to3的工具，一般叫做2to3.py，或者2to3-script.py，用以下命令转换一下2版本的脚本即可， ...

用了这个以后原始版本就能用了？
作者
Author: hdb    时间: 2023-5-11 11:13
biogon 发表于 2023-5-10 21:29
用了这个以后原始版本就能用了？

是的，测试过了，没有出现2020版的bug
作者
Author: biogon    时间: 2023-5-11 14:15
hdb 发表于 2023-5-11 11:13
是的，测试过了，没有出现2020版的bug

可以传上来，置个顶
作者
Author: hdb    时间: 2023-5-11 14:45
biogon 发表于 2023-5-11 14:15
可以传上来，置个顶

ok，放git里面了，https://github.com/hdb106/Elastic-py3.git
作者
Author: 牧生    时间: 2023-6-8 18:48
不知道能不能结合CP2K进行计算和分析呢
作者
Author: biogon    时间: 2023-6-9 11:57
牧生 发表于 2023-6-8 18:48
不知道能不能结合CP2K进行计算和分析呢

你给这个脚本添加上用cp2k算的功能就行
作者
Author: chenmuyean    时间: 2024-10-19 22:44
本帖最后由 chenmuyean 于 2024-10-19 22:45 编辑

博主您好，请问QE+ElaStic能否计算材料在不同温度和不同压强下的弹性常数，有没有什么参数可以设置的？网上的资料太少了，所以来请教您，望回复，感激不尽！
作者
Author: ZZZSSSYYY    时间: 2025-10-9 15:11
ZZG 发表于 2021-10-21 15:22
嗯嗯，我又在重新看作者写的文章，看能不能用elastic1.1.请问下，在使用之前是不是要进行这样的配置环境 ...

我用conda创建了一个python2的虚拟环境，在其中运行可以



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3