|
相互作用能脚本时间:2023-05-25 #!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"; |