二次差分电荷密度分析问题

zyj19831206

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.";

gaoming

您好，您这脚本能发我一份吗？我复制这里面的代码到MS里面运行报错，万分感谢！

南晓

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

嗯嗯谢谢我去试试

卡开发发

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

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

南晓

卡开发发 发表于 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.
请问一下这个是什么意思

南晓

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

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

南晓

楼主请问在运行中出现该问题是什么意思该咋解决请指教
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

sobereva

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

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

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

xiapin

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

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

loading

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

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

loading

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

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

卡开发发

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

你用的是CASTEP还是DMol3？

loading

同问题呀楼主，生成了densitydifference，但是没有密度图示，请问楼主最后如何解决的？

zyj19831206

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

好的，多谢社长。

sobereva

注意求助帖应当在帖子标题明确体现出提问、求助（见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等绘图。