首页 >> 专栏博客 >>其他未分类 >> 相互作用能脚本
详细内容

相互作用能脚本

时间:2023-05-25     作者:qcl【转载】

#!perl


# Author: Stephen Todd


# This script calculates the interaction energy between two layers in a structure.

# The layers have to be defined as sets called Layer 1 and Layer 2 - these are the 

# default names assigned by the Layer builder. 


# This also uses the default settings for Forcite - ie Universal with current charges. 

# You should change the settings if you want to use another forcefield or summation method.


# This script also assumes that the surface area is in the AB-plane.


# A study table is produced containing the layers, energies, interaction energy, and then

# the interaction energy per angstrom^2.


# Note. The input trajectory should have a large vaccuum, greater than the

# cut-off you are using in the non-bond calculation.


use strict;

use MaterialsScript qw(:all);


#################################################################################

# Begin user editable settings

#


my $filename = "200-321-6-1"; # Name of the trajectory without extension

my $energySettings = [

CurrentForcefield => "COMPASSIII",

Quality => "Ultra-fine",

AssignForcefieldTypes => "No", 

ChargeAssignment => "Use current"];

Modules->Forcite->ChangeSettings( $energySettings );

#

# End user editable settings

#################################################################################



# Defines the trajectory document


my $doc = $Documents{"$filename.xtd"};


# Create a new study table to hold the structures and energies and set the headings


my $newStudyTable = Documents->New("$filename"."_IntEnergy.std");

my $calcSheet = $newStudyTable->Sheets->Item(0);


$calcSheet->ColumnHeading(0) = "$filename Cell";

$calcSheet->ColumnHeading(1) = "Energy of  $filename Cell";

$calcSheet->ColumnHeading(2) = "Layer 1";

$calcSheet->ColumnHeading(3) = "Energy of Layer 1";

$calcSheet->ColumnHeading(4) = "Layer 2";

$calcSheet->ColumnHeading(5) = "Energy of layer 2";

$calcSheet->ColumnHeading(6) = "Layer 3";

$calcSheet->ColumnHeading(7) = "Energy of layer 3";

$calcSheet->ColumnHeading(8) = "Interaction Energy";

$calcSheet->ColumnHeading(9) = "Interaction Energy per angstrom^2";


# Get the total number of frames in the trajectory


my $numFrames = $doc->Trajectory->NumFrames;


print "Calculating interaction energy for $numFrames frames on $filename trajectory\n";


# Calculates the area of the surface based on the surface being in the A-B plane


my $length1 = $doc->Lattice3D->LengthA;

my $length2 = $doc->Lattice3D->LengthB;

my $surfaceArea = $length1 * $length2;

print "The surface area is $surfaceArea angstroms^2\n";


# Initialise Forcite. If you want to change the settings, you can load them or

# use a change settings command here.


my $forcite = Modules->Forcite;


# Starts to loop over the frames in the trajectory


for ( my $count=1000; $count<=$numFrames; $count++) {

 print "Calculating energy for frame $count\n";

     # Define the frame of the trajectory

     $doc->Trajectorys->Item(0)->CurrentFrame = $count;


     # Create three documents to hold the temporary layers

     my $allDoc = Documents->New("all.xsd");     

     my $layer1Doc = Documents->New("Layer1.xsd");

     my $layer2Doc = Documents->New("Layer2.xsd"); 

     my $layer3Doc = Documents->New("Layer3.xsd");

     # Create a copy of the structure in temporary documents and delete

     # the layers that you don't need.


     $allDoc->CopyFrom($doc);


     $layer1Doc->CopyFrom($doc);

     $layer1Doc->DisplayRange->Sets("surface")->Atoms->Delete;   


     $layer2Doc->CopyFrom($doc);

     $layer2Doc->DisplayRange->Sets("EDTA")->Atoms->Delete;

          

     $layer3Doc->CopyFrom($doc);

     $layer3Doc->DisplayRange->Sets("EDTA")->Atoms->Delete;

     $layer3Doc->DisplayRange->Sets("surface")->Atoms->Delete;


     # Put the structures in a study table for safekeeping


     $calcSheet->Cell($count-1, 0) = $allDoc;

     $calcSheet->Cell($count-1, 2) = $layer1Doc;

     $calcSheet->Cell($count-1, 4) = $layer2Doc;

     $calcSheet->Cell($count-1, 6) = $layer3Doc;


     #Calculate the energies for all the layers and put in the study table

     

     $forcite->Calculation->Run($allDoc, [Task => 'Energy']);

     $calcSheet->Cell($count-1, 1) = $allDoc->PotentialEnergy;

     

     $forcite->Calculation->Run($layer1Doc, [Task => 'Energy']);

     $calcSheet->Cell($count-1, 3) = $layer1Doc->PotentialEnergy;

     

     $forcite->Calculation->Run($layer2Doc, [Task => 'Energy']);

     $calcSheet->Cell($count-1, 5) = $layer2Doc->PotentialEnergy;

     

     $forcite->Calculation->Run($layer3Doc, [Task => 'Energy']);

     $calcSheet->Cell($count-1, 7) = $layer3Doc->PotentialEnergy;


     # Calculate the Interaction Energy from the cells in the Study Table

     

     my $totalEnergy = $calcSheet->Cell($count-1,1);    

     my $layer1Energy = $calcSheet->Cell($count-1,3);

     my $layer2Energy = $calcSheet->Cell($count-1,5);

     my $layer3Energy = $calcSheet->Cell($count-1,7);    

     

     my $interactionEnergy = $totalEnergy - ($layer1Energy + $layer2Energy) + $layer3Energy;

     $calcSheet->Cell($count-1, 8) = $interactionEnergy;

     

     my $interactionEnergyArea = $interactionEnergy / $surfaceArea;

     $calcSheet->Cell($count-1, 9) = $interactionEnergyArea;

          

     # Discard the temporary documents

     

     $allDoc->Discard;

     $layer1Doc->Discard;

     $layer2Doc->Discard;

     $layer3Doc->Discard;

}

print "Calculation complete.\n";


最新评论
请先登录才能进行回复登录
技术支持: CLOUD | 管理登录
seo seo