计算化学公社

标题: 二次差分电荷密度分析问题 [打印本页]

作者
Author: zyj19831206 时间: 2018-7-31 08:59
标题: 二次差分电荷密度分析问题
CASTEP里面貌似能够直接给出差分电荷密度，但是dmol里面貌似只能给出total density和spin density以及deform density,怎么通过这三个密度给出差分电荷密度呢，搞到一个脚本，但是具体不知道怎么使用它，脚本是perl语言的，希望高手指教下。
#!perl

use strict;
use MaterialsScript qw(:all);

# Authors: Stephen Todd and John Lear

# This script reads in the values of three fields, does some maths, and then
# writes the new field values to a new field. For this example, it calculates the
# electron density difference between a molecule and a surface.

# new field value = xsdDoc1 - ($xsdDoc2 + $xsdDoc3)

# The initial documents require a little setting up as you need to create field
# probes in each field you intend to probe. The example documents already contain field
# probes.

###############################################################################
# Begin user editable settings
#

my $xsdDoc1 = "DMol3All";                                              # Set the name of the first XSD doc

my $xsdDoc2 = "DMol3PdOnly";                                     # Set the name of the second XSD doc

my $xsdDoc3 = "DMol3COOnly";                                     # Set the name of the third XSD doc

my $fieldName = "DMol3 total electron density";             # Set the name of the original field

my $newXSD = "DensityDifferenceDmol3";                      # The name of the xsd for the new field

#
# End user editable settings
################################################################################

# Set the documents and fields up

my $doc1 = $Documents{"$xsdDoc1.xsd"};
my $field1 = $doc1->UnitCell->Fields("$fieldName");

my $doc2 = $Documents{"$xsdDoc2.xsd"};
my $field2 = $doc2->UnitCell->Fields("$fieldName");

my $doc3 = $Documents{"$xsdDoc3.xsd"};
my $field3 = $doc3->UnitCell->Fields("$fieldName");

# Copy the doc1 into a new xsd to work on - this copies the structure and field data

my $docNew = Documents->New("$newXSD.xsd");
$docNew->CopyFrom($doc1);
my $fieldNew = $docNew->UnitCell->Fields("$fieldName");

# Find extent of two of the original fields and check they match

$field1->GridSize =~ m/(\d+)\s*x\s*(\d+)\s*x\s*(\d+)/;
my $xExtent1 = $1;
my $yExtent1 = $2;
my $zExtent1 = $3;

print "$xsdDoc1 field size is:\n";
print "$xExtent1 $yExtent1 $zExtent1\n";

$field2->GridSize =~ m/(\d+)\s*x\s*(\d+)\s*x\s*(\d+)/;
my $xExtent2 = $1;
my $yExtent2 = $2;
my $zExtent2 = $3;

print "$xsdDoc2 field  size is:\n";
print "$xExtent2 $yExtent2 $zExtent2\n";

if (($xExtent1 != $xExtent2) || ($yExtent1 != $yExtent2) || ($zExtent1 != $zExtent2)) {
      die "The fields are not the same size and hence cannot be worked on";
}

# Define the field probes with the probe mode of nearest grid point. This means that the field
# probe jumps to the nearest grid point in the field. You can only edit field values if you
# are on grid points.

my $fieldProbe1;
eval q($fieldProbe1 = $field1->FieldProbe); die "Please add a field probe to $xsdDoc1" if $@;

$fieldProbe1->ProbeMode = "NearestGridPoint";

my $fieldProbe2;
eval q($fieldProbe2 = $field2->FieldProbe); die "Please add a field probe to $xsdDoc2" if $@;

$fieldProbe2->ProbeMode = "NearestGridPoint";

my $fieldProbe3;
eval q($fieldProbe3 = $field3->FieldProbe); die "Please add a field probe to $xsdDoc3" if $@;

$fieldProbe3->ProbeMode = "NearestGridPoint";

my $fieldProbeNew = $fieldNew->FieldProbe;
$fieldProbeNew->ProbeMode = "NearestGridPoint";

# Loop over the different fields, reading in the values and setting the new field value.

for (my $x=0; $x<=$xExtent1; ++$x) {

      for (my $y=0; $y<=$yExtent1; ++$y) {

            for (my $z=0; $z<=$zExtent1; ++$z) {

            # Set the field probe to the ProbeVoxel Position and get the field value

            $fieldProbe1->ProbeVoxelPosition = Point(X => $x, Y=> $y, Z => $z);
            my $fieldValue1 = $fieldProbe1->FieldValue;

            $fieldProbe2->ProbeVoxelPosition = Point(X => $x, Y=> $y, Z => $z);
            my $fieldValue2 = $fieldProbe2->FieldValue;

            $fieldProbe3->ProbeVoxelPosition = Point(X => $x, Y=> $y, Z => $z);
            my $fieldValue3 = $fieldProbe3->FieldValue;

            # From the field values, calculate the new field value

            my $newFieldValue = $fieldValue1 - ( $fieldValue2 + $fieldValue3 );

            # Move the field probe in the new field and set the value

            $fieldProbeNew->ProbeVoxelPosition = Point(X => $x, Y=> $y, Z => $z);
            $fieldProbeNew->FieldValue = $newFieldValue;

         }
      }
}

print "Calculation is complete.";

作者
Author: 卡开发发 时间: 2018-7-31 09:40
三个结构（相互作用结构+两个片段）分别用相同的晶格、截断半径、k点等计算密度，然后分别给xsdDoc1、2、3赋值你三个任务的名字运行脚本应该就行了。实际也不是非得采用perl做不可，计算密度所生成的grd文件里面的数据对应做差导出新的grd文件导入到ms也可以实现。

作者
Author: zyj19831206 时间: 2018-7-31 10:04

卡开发发发表于 2018-7-31 09:40
三个结构（相互作用结构+两个片段）分别用相同的晶格、截断半径、k点等计算密度，然后分别给xsdDoc1、2、3 ...

分别是这三个文件，但是我把程序与这些文件放在一起，运行出来啥也没有啊，就这个结果：
DMol3All field size is:
16 12 44
DMol3PdOnly field size is:
16 12 44
Calculation is complete.

作者
Author: 卡开发发 时间: 2018-7-31 12:18

zyj19831206 发表于 2018-7-31 10:04
分别是这三个文件，但是我把程序与这些文件放在一起，运行出来啥也没有啊，就这个结果：
DMol3All field ...

检查一下三个文件的display sytle-field是否被导入，如果没有的话到analysis里面把电子密度import一下，按道理应该会生成DensityDifferenceDmol3.xsd并且包含电荷密度差的field。

作者
Author: zyj19831206 时间: 2018-7-31 14:04

卡开发发发表于 2018-7-31 12:18
检查一下三个文件的display sytle-field是否被导入，如果没有的话到analysis里面把电子密度import一下， ...

只有densitydifference这个set，没有具体的密度图示。

作者
Author: 卡开发发 时间: 2018-7-31 15:38

zyj19831206 发表于 2018-7-31 14:04
只有densitydifference这个set，没有具体的密度图示。

# The initial documents require a little setting up as you need to create field
# probes in each field you intend to probe.

看下有没有注意到这个问题。

作者
Author: sobereva 时间: 2018-8-1 05:33
注意求助帖应当在帖子标题明确体现出提问、求助（见http://bbs.keinsci.com/thread-9348-1-1.html），给你改了

你可以产生出不同体系/电子态的.grd文件（一种记录格点数据的文件，详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》http://sobereva.com/379），载入到Multiwfn里，利用Multiwfn的主功能13的子功能11对格点数据进行相互求差，Multiwfn手册4.13节有例子。算完的格点数据可以在Multiwfn里直接观看，也可以导出成cube文件，用VMD、gview等绘图。

作者
Author: zyj19831206 时间: 2018-8-1 21:17

sobereva 发表于 2018-8-1 05:33
注意求助帖应当在帖子标题明确体现出提问、求助（见http://bbs.keinsci.com/thread-9348-1-1.html），给你 ...

好的，多谢社长。

作者
Author: loading 时间: 2021-3-4 16:26
同问题呀楼主，生成了densitydifference，但是没有密度图示，请问楼主最后如何解决的？

作者
Author: 卡开发发 时间: 2021-3-4 18:49

loading 发表于 2021-3-4 16:26
同问题呀楼主，生成了densitydifference，但是没有密度图示，请问楼主最后如何解决的？

你用的是CASTEP还是DMol3？

作者
Author: loading 时间: 2021-3-5 15:02

卡开发发发表于 2021-3-4 18:49
你用的是CASTEP还是DMol3？

dmol3，问题已经解决啦，右键，display style,isosurface,右下角，type选±，isovalue值改小一些就可以显示了

作者
Author: loading 时间: 2021-3-5 15:03

zyj19831206 发表于 2018-8-1 21:17
好的，多谢社长。

楼主，我解决了，右键，display style,isosurface,右下角，type选±，isovalue值改小一些就可以显示了

作者
Author: xiapin 时间: 2021-7-28 16:29
本帖最后由 xiapin 于 2021-7-28 16:31 编辑

sobereva 发表于 2018-8-1 05:33
注意求助帖应当在帖子标题明确体现出提问、求助（见http://bbs.keinsci.com/thread-9348-1-1.html），给你 ...

社长，请问借助Multiwfn软件做差得到的.cub文件能做出差分电荷图吗？
将.cub导进VMD，Gview都没有做出差分电荷密度图。请问能用VESTA作图吗？导进VESTA后发现得到一些黄色的点点，如图，也是没有得到想要的图。请社长帮忙解答一下作图的疑惑。

作者
Author: sobereva 时间: 2021-7-29 02:09

xiapin 发表于 2021-7-28 16:29
社长，请问借助Multiwfn软件做差得到的.cub文件能做出差分电荷图吗？
将.cub导进VMD，Gview都没有做出差 ...

Multiwfn产生的记录密度差的cube文件不可能用VMD作不出来等值面图，没出来肯定是你操作不对。诸如cube文件产生流程有问题，isovalue设得不当，等等

你当前帖的如果是等值面图，肯定是cube文件里记录的根本就不是密度差数据

作者
Author: 南晓 时间: 2022-4-21 16:34
楼主请问在运行中出现该问题是什么意思该咋解决请指教
1 field size is:
99 79 44
2 field size is:
48 30 43
The fields are not the same size and hence cannot be worked on

作者
Author: 南晓 时间: 2022-4-22 09:18

loading 发表于 2021-3-4 16:26
同问题呀楼主，生成了densitydifference，但是没有密度图示，请问楼主最后如何解决的？

您好想请问一下我在运行该脚本时，在第82行停止了，显示
please add a field probe to 1 at -e line 82
请问您知道这个该怎么解决吗，多谢

作者
Author: 南晓 时间: 2022-4-22 09:30

卡开发发发表于 2018-7-31 15:38
看下有没有注意到这个问题。

您好想请教一个问题我在运行该脚本时，出现问题显示
1 field size is:
282 219 50
2 field size is:
282 219 50
please add a fied probe to 1 at -e line 82.
请问一下这个是什么意思

作者
Author: 卡开发发 时间: 2022-4-22 10:44

南晓发表于 2022-4-22 09:30
您好想请教一个问题我在运行该脚本时，出现问题显示
1 field size is:
282 219 50

在ms的visualizer里面加个probe就行了。

作者
Author: 南晓 时间: 2022-4-22 11:15

卡开发发发表于 2022-4-22 10:44
在ms的visualizer里面加个probe就行了。

嗯嗯谢谢我去试试

作者
Author: gaoming 时间: 2024-3-21 20:50
您好，您这脚本能发我一份吗？我复制这里面的代码到MS里面运行报错，万分感谢！

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