|
|
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.";
|
|
发表于 Post on 2018-7-31 08:59:54
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author
发表于 Post on 2018-8-1 05:33:14
