求助分波态密度计算输出文件E缺失

Clifford

请教各位老师，我最近在计算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，请教各位大佬是否遇到过相似的问题，望赐教！

zjxitcc

如果你确定你的计算无误，那就是程序输出格式没考虑周全的问题，你手动把它改成0.000E+00即可

Clifford

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

您好！感谢解答，我尝试这样修改过，但这样的问题还有很多，在多个输出文件里都存在。请教老师有没有办法调节输出数据的小数点显示位数，或者在输入文件做某种定义使负多少次方以下的数为零？

卡开发发

Clifford 发表于 2022-2-5 14:27
您好！感谢解答，我尝试这样修改过，但这样的问题还有很多，在多个输出文件里都存在。请教老师有没有办法 ...

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

1. import re
   
2. import sys
   
3. 
   
4. """
   
5. 使用方法: pdos.py [name].pdos_atm[atm]_wfc[orb]
   
6. """
   
7. 
   
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截断的影响。

Clifford

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

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

卡开发发

Clifford 发表于 2022-2-5 23:24
感谢卡开发发大佬指教，太厉害了！我用您的方法替换了相应字符串，就能顺利作出了pdos图。
但 ...

这个问题说不好，搞不好是什么特定的数值误差所引起的。要是实在不需要展宽稳妥还是直接用四面体积分方法来搞。

Clifford

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

嗯嗯，是的，学习了

sobereva

Fortran里可以明确指定指数部分输出多少位而且保留E符号
比如5E100用E16.5来输出，就会导致输出是0.50000+101，即E没了。如果用E16.5E3，则输出为0.50000E+101。因此只需要改一行源代码即可解决

Clifford

sobereva 发表于 2022-2-6 02:41
Fortran里可以明确指定指数部分输出多少位而且保留E符号
比如5E100用E16.5来输出，就会导致输出是0.50000+ ...

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