CP2K输出的静电势问题

panger

       手册上的V_HARTREE_CUBE里面这句话是什么意思？“Note that by convention the potential has opposite sign than the expected physical one.”
计算得到的静电势的图是这样的。这和以往得到的静电势不太一样。为什么没有正负值不同的区域呢？对cp2k了解不深，不明白这张图能能得到什么信息。

丁越

卡开发发 发表于 2021-8-2 12:14
单位问题，V_h的积分应该dx和dy都是Bohr单位。但这样换算下来好像还是有些偏差，原因不明。

谢谢老师，我知道了

卡开发发

丁越 发表于 2021-8-2 11:43
太感谢老师了这是我用Multiwfn处理的步骤：载入静电势cube文件后，13- 18（Z,a）- 6，这里导出的loc ...

单位问题，V_h的积分应该dx和dy都是Bohr单位。但这样换算下来好像还是有些偏差，原因不明。

丁越

卡开发发 发表于 2021-8-2 11:24
这个python3代码(cube_average_z.py)我随便写写作为一个测试：
执行的方法是：
如果对平均静电势取最大 ...

太感谢老师了这是我用Multiwfn处理的步骤：载入静电势cube文件后，13- 18（Z,a）- 6，这里导出的locintcurve.txt数据中将第三列除以Axy面积（17.74*10.2425），然后以该txt中第二列数据Z为X轴，除了面积之后的第三列数据为Y轴作图，就得到上面的图像了。老师您看有什么问题么？

卡开发发

丁越 发表于 2021-8-2 08:29
好的，麻烦老师看一下

这个python3代码(cube_average_z.py)我随便写写作为一个测试：

1. #本程序只能用在正交晶格。
   
2. #非正交晶格理论上可以通过修改本程序实现。
   
3. import numpy as np
   
4. import sys
   
5. 
   
6. def average_cube_z(fname):
   
7.         f=open(fname,'r+')
   
8.         text=f.read()
   
9.         f.close()
   
10.         #cube文件按行进行切割。
    
11.         text=text.split('\n')
    
12.         #第3~6行为网格的轴向上的格点数和间距，两者对应相乘可以得到对应的晶格参数。
    
13.         lat_info=np.array([t.split() for t in text[3:6]])
    
14.         ngrids=np.array(lat_info[:,0],dtype=int)
    
15.         intv=np.array(lat_info[:,1:],dtype=float)
    
16.         lat=ngrids*intv
    
17.         #第2行为原子数和原点坐标信息，我们只取原子数的信息即可。
    
18.         natms=int(text[2].split()[0])
    
19.         #第6行起，后续有与原子数相等的行数为原子和原子坐标信息，这里我们不需要跳过去。
    
20.         grids=' '.join(text[6+natms:])
    
21.         #将数组的形状按照给定的网格尺寸重新排，注意cube格式变化最快的方向是z方向。
    
22.         grids=np.array(grids.split(),dtype=float)
    
23.         grids=grids.reshape(ngrids)
    
24.         #如果要更改为x为最快变化方向可以使用：
    
25.         #grids=np.transpose(grids,(2,1,0))
    
26.         #这里没必要这么处理。
    
27.         #---------------------------------------
    
28.         #求某个特定的z下静电势网格在xy面上的平均。
    
29.         #这么做可行的原因是，网格为矩形且均匀的情况下，利用矩形积分的公式：
    
30.         #∫∫ρ(x,y,z)dxdy/∫∫dxdy=(∑∑ρ_xyz)/n_x/n_y。
    
31.         #当然也可以用精度更高的积分方法或用FFT，这里只是定性说明。
    
32.         average=np.mean(grids,axis=(0,1))
    
33.         #产生Bohr单位下z的坐标网格。
    
34.         z=np.linspace(0,np.linalg.norm(lat[2]),ngrids[2],endpoint=False)
    
35.         #将数据按照z V_H的格式打印。
    
36.         head='%25s %25s\n'%('z','V_H')
    
37.         text='\n'.join(['% 16.8f % 16.8f'%(i,a) for i,a in zip(z,average)])
    
38.         text=head+text
    
39.         return text
    
40. 
    
41. text=average_cube_z(sys.argv[1])
    
42. print(text)

复制代码

执行的方法是：

1. python cube_average_z.py Cu_111-v_hartree-1_0.cube > potentials.dat

复制代码

如果对平均静电势取最大得到的结果是：Vh=0.10814 Hartree，
然后按照你前面给的数据：Ef=0.3029-0.3829=-0.08 Hartree，
然后按照eV单位计算功函数：(Vh-Ef)*27.211=5.11947754 eV
其他的说明在代码块部分说的应该也比较明确了，和官网给的程序结果的差异大概率来源于积分的方法，我这边用的是粗糙的矩形积分方法。可以用更高级别的基于插值的数值积分，也可以用FFT，原则上经过修改可以有更高的精度并可用于γ角非90°的单斜结构。至于之前为啥做不对，我也不知道你具体的处理步骤。

丁越

卡开发发 发表于 2021-8-1 20:17
我不确定你I(z)是否处理正确？要不然多给一些原始数据咱们可以一块讨论下问题在哪。

好的，麻烦老师看一下

卡开发发

丁越 发表于 2021-8-1 15:07
卡开发发老师，我使用您的方法计算了Cu(111)表面的功函数，真空层大小取了35Å，下面是做出的图像。 ...

我不确定你I(z)是否处理正确？要不然多给一些原始数据咱们可以一块讨论下问题在哪。

丁越

卡开发发 发表于 2017-4-3 09:17
1、Multiwfn计算local integral curve采用的公式是
I(z)=∫∫p(x,y,z)dxdy
平均静电势的量纲是能量的量 ...

卡开发发老师，我使用您的方法计算了Cu(111)表面的功函数，真空层大小取了35Å，下面是做出的图像。可是计算出的功函数大的离谱，不知道是不是对公示理解错了导致图画错了？（Axy就是Cu xy表面的面积这没错吧，建模时我已经把模型弄成正交晶系了，并且最后也除了这个表面积）。另外，我也用官网给的cubecruncher脚本处理计算了功函数，得到的值是5.25 eV（这是换算单位后的），与维基百科给出的Cu参考值4.53~5.10 eV接近。老师能指点下我哪里搞错了么？

丁越

panger 发表于 2017-4-2 18:27
另外还想请教您两个问题。文献中说功函计算是使用真空层中的静电势（选取真空层中部）减去材料的费米能级 ...用cp2k官网给的这个脚本按照如下方法也可以得到真空能级。

丁越

panger 发表于 2017-4-2 18:27
另外还想请教您两个问题。文献中说功函计算是使用真空层中的静电势（选取真空层中部）减去材料的费米能级 ...

panger

sobereva 发表于 2017-9-29 16:42
这看起来没问题

谢谢sob老师。

sobereva

panger 发表于 2017-9-29 16:16
谢谢您。目前cube文件还没算出来，等出来我马上看一下。可能我没表述清楚，我的意思是构建超胞的时候建成 ...

这看起来没问题

panger

sobereva 发表于 2017-9-28 23:58
你说的情况我不太清楚。你看一下.cub文件开头，只要三个平移矢量正分别对应x,y,z方向，就可以用Multiwf ...

谢谢您。目前cube文件还没算出来，等出来我马上看一下。可能我没表述清楚，我的意思是构建超胞的时候建成这样的，而不是软件默认的90，90，120的那种情况。

sobereva

panger 发表于 2017-9-28 16:51
sob老师，如果我在切面的时候直接修改一下，取成1*根号3的单胞，而不是1*1的单胞，这时是不是就可以直接 ...

你说的情况我不太清楚。你看一下.cub文件开头，只要三个平移矢量正分别对应x,y,z方向，就可以用Multiwfn分析

panger

sobereva 发表于 2017-4-3 01:46
可以

sob老师，如果我在切面的时候直接修改一下，取成1*根号3的单胞，而不是1*1的单胞，这时是不是就可以直接用multiwfn来分析结果了？