计算化学公社

标题: CP2K输出的静电势问题 [打印本页]

作者
Author: panger 时间: 2017-4-1 20:56
标题: CP2K输出的静电势问题
手册上的V_HARTREE_CUBE里面这句话是什么意思？“Note that by convention the potential has opposite sign than the expected physical one.”
计算得到的静电势的图是这样的。这和以往得到的静电势不太一样。为什么没有正负值不同的区域呢？对cp2k了解不深，不明白这张图能能得到什么信息。

作者
Author: sobereva 时间: 2017-4-2 01:31
静电势正值区域的数量级范围往往显著大于负值区域的数量级范围。在显示静电势等值面图的时候，如果isovalue不恰当设置，而且没有让程序同时输出两种符号的话，会看不到负值区域。
在此之前应该先看一下cube文件的最负的数值是多少，载入Multiwfn的时候屏幕上会立刻提示负值区域的加和（不过是否兼容CP2K的cube文件没测试过）

作者
Author: panger 时间: 2017-4-2 08:16
本帖最后由 panger 于 2017-4-2 08:19 编辑

sobereva 发表于 2017-4-2 01:31
静电势正值区域的数量级范围往往显著大于负值区域的数量级范围。在显示静电势等值面图的时候，如果isovalue ...

谢谢sob老师，Multiwfn能读cp2k的cube。您说的是这个吧？我该怎么设置，负值区域才能显示呢？还是本身结果就有问题？实在是迷糊了。

作者
Author: sobereva 时间: 2017-4-2 08:28

panger 发表于 2017-4-2 08:16
谢谢sob老师，Multiwfn能读cp2k的cube。您说的是这个吧？我该怎么设置，负值区域才能显示呢？还是本身结 ...

你可以压缩一下传到网盘上，我下载看看。如果压缩后小于3MB传论坛里也行

作者
Author: panger 时间: 2017-4-2 08:53

sobereva 发表于 2017-4-2 08:28
你可以压缩一下传到网盘上，我下载看看。如果压缩后小于3MB传论坛里也行

我给您sina的邮箱发了邮件可以吗？您受累看看。

作者
Author: sobereva 时间: 2017-4-2 08:55

panger 发表于 2017-4-2 08:53
我给您sina的邮箱发了邮件可以吗？您受累看看。

可

作者
Author: sobereva 时间: 2017-4-2 16:11
你的格点文件没问题，是你显示设置问题。这里用VMD读入，分别建立正值和负值的显示方式，等值面数值分别为0.2（红色）和-0.2（青色）
(, 下载次数 Times of downloads: 63)

作者
Author: panger 时间: 2017-4-2 18:21

sobereva 发表于 2017-4-2 16:11
你的格点文件没问题，是你显示设置问题。这里用VMD读入，分别建立正值和负值的显示方式，等值面数值分别为0 ...

太感谢您了，不过我做出和您一样的图等值面的数值选的是0.3和-0.1才能出现您这个效果？是我的软件版本有问题吗？我用VMD和VESTA做出来的等值面都是这个数值才能出您这样的效果。

作者
Author: panger 时间: 2017-4-2 18:27
本帖最后由 panger 于 2017-4-3 00:56 编辑

sobereva 发表于 2017-4-2 16:11
你的格点文件没问题，是你显示设置问题。这里用VMD读入，分别建立正值和负值的显示方式，等值面数值分别为0 ...

另外还想请教您两个问题。文献中说功函计算是使用真空层中的静电势（选取真空层中部）减去材料的费米能级，这个真空层中部的静电势能拿到吗？从cube文件。我用Multiwfn手册4.13.6里面的方法绘制的local integral curve纵坐标好大，这种图有没有意义呢？

作者
Author: panger 时间: 2017-4-3 01:01
跟这张图是不是有相似的意思呢？

作者
Author: sobereva 时间: 2017-4-3 01:18

panger 发表于 2017-4-2 18:21
太感谢您了，不过我做出和您一样的图等值面的数值选的是0.3和-0.1才能出现您这个效果？是我的软件版本有 ...

能显示出来就行，也许是我记错了

作者
Author: sobereva 时间: 2017-4-3 01:19

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

当前的cube文件的三个平移矢量不是恰好是在x,y,z方向，即格点不是立方的，这种情况cube文件虽然能载入Multiwfn，但是数值不正常。

功函数这些固体物理相关的概念还是问开卡发发吧。

作者
Author: panger 时间: 2017-4-3 01:22
本帖最后由 panger 于 2017-4-3 01:25 编辑

sobereva 发表于 2017-4-3 01:19
当前的cube文件的三个平移矢量不是恰好是在x,y,z方向，即格点不是立方的，这种情况cube文件虽然能载入M ...

谢谢您。太感谢了。您的意思是说，这种非立方的格点，还不能用Multiwfn分析？

作者
Author: sobereva 时间: 2017-4-3 01:37

panger 发表于 2017-4-3 01:22
谢谢您。太感谢了。您的意思是说，这种非立方的格点，还不能用Multiwfn分析？

对

作者
Author: panger 时间: 2017-4-3 01:41
本帖最后由 panger 于 2017-4-3 01:47 编辑

sobereva 发表于 2017-4-3 01:19
当前的cube文件的三个平移矢量不是恰好是在x,y,z方向，即格点不是立方的，这种情况cube文件虽然能载入M ...

另外，如果对于我这种情况，如果我想用Multiwfn来做，我可不可以在排除周围活性中心(我指金属团簇)影响的情况下，比如我把载体建的比较大。人为的建立一个立方的晶格来算呢？类似于加了个框的团簇模型。这样是否能得到相对合理的结果？

作者
Author: sobereva 时间: 2017-4-3 01:46

panger 发表于 2017-4-3 01:41
另外，如果对于我这种情况，如果我想用Multiwfn来做，我可不可以在排除周围活性中心(我指金属团簇)影响影 ...

可以

作者
Author: panger 时间: 2017-4-3 02:05
本帖最后由 panger 于 2017-4-3 02:09 编辑

谢谢sob老师，这么晚，打扰您了。

作者
Author: 卡开发发 时间: 2017-4-3 09:17

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

1、Multiwfn计算local integral curve采用的公式是
I(z)=∫∫p(x,y,z)dxdy
平均静电势的量纲是能量的量纲，和上述定义差一个面积的量纲，当Z方向与XY面垂直：
VH=∫∫p(x,y,z)dxdy/∫∫dxdy=I(z)/Axy，其中Axy是XY的面积。
最后注意单位换算。

2、Multiwfn应该只支持立方的格子，就目前而言，方便的做法就是把晶格基矢取为正交的就可以了，对于二维系统只要将A和B的两个方向进行线性组合就行。

作者
Author: panger 时间: 2017-9-28 16:51

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

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

作者
Author: sobereva 时间: 2017-9-28 23:58

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

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

作者
Author: panger 时间: 2017-9-29 16:16
本帖最后由 panger 于 2017-9-29 16:17 编辑

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

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

作者
Author: sobereva 时间: 2017-9-29 16:42

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

这看起来没问题

作者
Author: panger 时间: 2017-9-29 16:45

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

谢谢sob老师。

作者
Author: 丁越 时间: 2021-7-15 19:17
本帖最后由丁越于 2021-7-29 08:53 编辑

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

作者
Author: 丁越 时间: 2021-7-28 23:04
本帖最后由丁越于 2021-8-1 15:12 编辑

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

作者
Author: 丁越 时间: 2021-8-1 15:07

卡开发发发表于 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接近。老师能指点下我哪里搞错了么？

作者
Author: 卡开发发 时间: 2021-8-1 20:17
本帖最后由卡开发发于 2021-8-1 20:39 编辑

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

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

作者
Author: 丁越 时间: 2021-8-2 08:29

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

好的，麻烦老师看一下

作者
Author: 卡开发发 时间: 2021-8-2 11:24

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

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

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

复制代码

执行的方法是：

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°的单斜结构。至于之前为啥做不对，我也不知道你具体的处理步骤。

作者
Author: 丁越 时间: 2021-8-2 11:43

卡开发发发表于 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轴作图，就得到上面的图像了。老师您看有什么问题么？

作者
Author: 卡开发发 时间: 2021-8-2 12:14

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

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

作者
Author: 丁越 时间: 2021-8-2 14:04

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

谢谢老师，我知道了

欢迎光临计算化学公社 (http://bbs.keinsci.com/)

Powered by Discuz! X3.3