计算化学公社

标题: 求助分波态密度计算输出文件E缺失 [打印本页]
作者
Author: Clifford    时间: 2022-2-5 14:13
标题: 求助分波态密度计算输出文件E缺失
请教各位老师，我最近在计算pdos时遇到一个奇怪的问题，正常进行过scf，projwfc计算后得到了一系列pwscf.pdos_atm#2(xx)_wfc#1(s)文件，做图时提示could not convert string to float: '0.198-100'，经查找发现pwscf.pdos_atm#1(Pb)_wfc#1(d)中存在一些缺失E的string，导致数据被判断为string，截取部分如下：（在其它分波文件中还有很多类似的情况）
   9.383  0.798E-84  0.213E-85  0.965E-85  0.965E-85  0.581E-84  0.238E-86
   9.393  0.224E-83  0.335E-85  0.215E-84  0.215E-84  0.178E-83  0.255E-86
   9.403  0.226E-83  0.431E-85  0.217E-84  0.217E-84  0.178E-83  0.497E-86
   9.413  0.156E-83  0.247E-85  0.131E-84  0.131E-84  0.127E-83  0.287E-86
   9.423  0.167E-85  0.953E-86  0.239E-86  0.239E-86  0.198-100  0.241E-86
   9.433  0.195E-84  0.434E-86  0.329E-85  0.329E-85  0.112E-84  0.131E-85
   9.443  0.264E-84  0.316E-85  0.535E-85  0.535E-85  0.112E-84  0.135E-85
   9.453  0.264E-84  0.316E-85  0.535E-85  0.535E-85  0.112E-84  0.135E-85
   9.463  0.950E-85  0.338E-85  0.247E-85  0.247E-85  0.113E-85  0.575E-87
   9.473  0.469E-84  0.509E-85  0.692E-85  0.692E-85  0.275E-84  0.463E-86

之前很长时间做这样的计算没遇到过这个问题，qe版本为7.0，请教各位大佬是否遇到过相似的问题，望赐教！

作者
Author: zjxitcc    时间: 2022-2-5 14:19
如果你确定你的计算无误，那就是程序输出格式没考虑周全的问题，你手动把它改成0.000E+00即可
作者
Author: Clifford    时间: 2022-2-5 14:27
zjxitcc 发表于 2022-2-5 14:19
如果你确定你的计算无误，那就是程序输出格式没考虑周全的问题，你手动把它改成0.000E+00即可

您好！感谢解答，我尝试这样修改过，但这样的问题还有很多，在多个输出文件里都存在。请教老师有没有办法调节输出数据的小数点显示位数，或者在输入文件做某种定义使负多少次方以下的数为零？
作者
Author: 卡开发发    时间: 2022-2-5 20:50
本帖最后由 卡开发发 于 2022-2-5 21:28 编辑
Clifford 发表于 2022-2-5 14:27
您好！感谢解答，我尝试这样修改过，但这样的问题还有很多，在多个输出文件里都存在。请教老师有没有办法 ...

这个缺陷是Fortran的'E'输出格式对指数通常只能输出2位，有两种途径可以处理这个问题。
一种途径是，可以写个python程序，例如：
  1. import re
  2. import sys

  4. """
  5. 使用方法: pdos.py [name].pdos_atm[atm]_wfc[orb]
  6. """

  8. #读取文件
  9. f=open(sys.argv[1],'r+')
  10. inp=f.read()
  11. f.close()
  12. """
  13. 将inp中的0.XXX-XXX更换为0,000E+00。可以这么做的原因如下：
  14. 第一列为F8.3格式，负号不可能出现在两个数字间。
  15. 第二列之后的浮点数，负号不应该出现在两个数字间，只可能：
  16. （1）要么是在一个空格和数字之间，要么是在E和数字之间。
  17. 对+E100的情况可能需要另外讨论。
  18. """
  19. out=re.sub('0.\d{3}-\d{3}','0.000E+00',inp)
  20. #写入文件
  21. f=open(sys.argv[1],'w+')
  22. f.write(out)
  23. f.close()
复制代码
然后对异常数字进行替换。当然上述也不是唯一的做法，比如利用python的try-except结构对异常数据进行处理。
另一种方式则是通过修改pp.x自身输出规则来实现，在PP/src/partialdos.f90下，将pdos累加为ldos这一段，
  1.         !
  2.         ! ldos = PDOS summed over m (m=-l:+l)
  3.         !
  4.         ldos  (:,:,:) = 0.d0
  5.         DO ik=1,nkseff
  6.            DO ie= 0, ne
  7.               DO is=1, nspin
  8.                  DO m=1,2 * nlmchi(nwfc)%l + 1
  9.                     ldos  (ie, is, ik) = ldos  (ie, is, ik) + pdos(ie,nwfc+m-1,is,ik)
  10.                  ENDDO
  11.               ENDDO
  12.            ENDDO
  13.         ENDDO
复制代码

这里过后可以简单写个pdos和ldos的循环，将pdos和ldos<1e-100的全部替换为0.0d0即可。在这之后去做处理的原因是我们不希望LDOS在累加的时候收到PDOS截断的影响。



作者
Author: Clifford    时间: 2022-2-5 23:24
卡开发发 发表于 2022-2-5 20:50
这个缺陷是Fortran的'E'输出格式对指数通常只能输出2位，有两种途径可以处理这个问题。
一种途径是，可 ...

感谢卡开发发大佬指教，太厉害了！我用您的方法替换了相应字符串，就能顺利作出了pdos图。
但新的问题出现了，所有pdos曲线的态密度值基本为零，多番尝试后发现这和scf中degauss太小(1.0d-9)有关，增大degauss(1.0d-2)，或者最好采用Tetrahedron method就根本不会有指数项低于-100这个问题。
作者
Author: 卡开发发    时间: 2022-2-5 23:59
Clifford 发表于 2022-2-5 23:24
感谢卡开发发大佬指教，太厉害了！我用您的方法替换了相应字符串，就能顺利作出了pdos图。
但 ...

这个问题说不好，搞不好是什么特定的数值误差所引起的。要是实在不需要展宽稳妥还是直接用四面体积分方法来搞。
作者
Author: Clifford    时间: 2022-2-6 01:01
卡开发发 发表于 2022-2-5 23:59
这个问题说不好，搞不好是什么特定的数值误差所引起的。要是实在不需要展宽稳妥还是直接用四面体积分方法 ...

嗯嗯，是的，学习了
作者
Author: sobereva    时间: 2022-2-6 02:41
Fortran里可以明确指定指数部分输出多少位而且保留E符号
比如5E100用E16.5来输出，就会导致输出是0.50000+101，即E没了。如果用E16.5E3，则输出为0.50000E+101。因此只需要改一行源代码即可解决
作者
Author: Clifford    时间: 2022-2-6 08:38
sobereva 发表于 2022-2-6 02:41
Fortran里可以明确指定指数部分输出多少位而且保留E符号
比如5E100用E16.5来输出，就会导致输出是0.50000+ ...

好办法，谢谢sob老师指点！



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