请选择 进入手机版 | 继续访问电脑版
第9届北京科音分子动力学与GROMACS培训班将于4月17~20日于北京举办,请点击此链接查看培训详情,欢迎参加和相互转告!

计算化学公社

 找回密码
 现在注册!
查看: 3942|回复: 11

[其它量化程序] 二次差分电荷密度分析问题

[复制链接]

547

帖子

2

威望

3508

eV
积分
4095

Level 6 (一方通行)

发表于 2018-7-31 08:59:54 | 显示全部楼层 |阅读模式
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.";


2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

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

评分

参与人数 1eV +3 收起 理由
wangyj + 3 谢谢

查看全部评分

满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

547

帖子

2

威望

3508

eV
积分
4095

Level 6 (一方通行)

 楼主| 发表于 2018-7-31 10:04:43 | 显示全部楼层
卡开发发 发表于 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.

DMol3All.xsd

103.87 KB, 下载次数: 7

DMol3COOnly.xsd

77 KB, 下载次数: 3

DMol3PdOnly.xsd

101.28 KB, 下载次数: 4

2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

发表于 2018-7-31 12:18:39 | 显示全部楼层
zyj19831206 发表于 2018-7-31 10:04
分别是这三个文件,但是我把程序与这些文件放在一起,运行出来啥也没有啊,就这个结果:
DMol3All field ...

检查一下三个文件的display sytle-field是否被导入,如果没有的话到analysis里面把电子密度import一下,按道理应该会生成DensityDifferenceDmol3.xsd并且包含电荷密度差的field。
满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

547

帖子

2

威望

3508

eV
积分
4095

Level 6 (一方通行)

 楼主| 发表于 2018-7-31 14:04:51 | 显示全部楼层
卡开发发 发表于 2018-7-31 12:18
检查一下三个文件的display sytle-field是否被导入,如果没有的话到analysis里面把电子密度import一下, ...

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

2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

发表于 2018-7-31 15:38:15 | 显示全部楼层
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.

看下有没有注意到这个问题。
满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

3万

帖子

25

威望

3万

eV
积分
65563

管理员

公社社长+计算化学玩家

发表于 2018-8-1 05:33:14 | 显示全部楼层
注意求助帖应当在帖子标题明确体现出提问、求助(见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等绘图。
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395(已满),2号:466017436(已满)。3号:764390338(可加),合计8000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

547

帖子

2

威望

3508

eV
积分
4095

Level 6 (一方通行)

 楼主| 发表于 2018-8-1 21:17:34 | 显示全部楼层
sobereva 发表于 2018-8-1 05:33
注意求助帖应当在帖子标题明确体现出提问、求助(见http://bbs.keinsci.com/thread-9348-1-1.html),给你 ...

好的,多谢社长。

13

帖子

0

威望

42

eV
积分
55

Level 2 能力者

发表于 2021-3-4 16:26:10 | 显示全部楼层
同问题呀楼主,生成了densitydifference,但是没有密度图示,请问楼主最后如何解决的?

2265

帖子

3

威望

8626

eV
积分
10951

Level 6 (一方通行)

Ab Initio Amateur

发表于 2021-3-4 18:49:18 | 显示全部楼层
loading 发表于 2021-3-4 16:26
同问题呀楼主,生成了densitydifference,但是没有密度图示,请问楼主最后如何解决的?

你用的是CASTEP还是DMol3?
满招损,谦受益。热衷于理论和方法研究水平不高但欢迎讨论。

13

帖子

0

威望

42

eV
积分
55

Level 2 能力者

发表于 2021-3-5 15:02:28 | 显示全部楼层
卡开发发 发表于 2021-3-4 18:49
你用的是CASTEP还是DMol3?

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

13

帖子

0

威望

42

eV
积分
55

Level 2 能力者

发表于 2021-3-5 15:03:20 | 显示全部楼层
zyj19831206 发表于 2018-8-1 21:17
好的,多谢社长。

楼主,我解决了,右键,display style,isosurface,右下角,type选±,isovalue值改小一些就可以显示了
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )

GMT+8, 2021-4-15 07:06 , Processed in 0.225766 second(s), 29 queries .

快速回复 返回顶部 返回列表