#!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 = "Layer";		# Name of the trajectory without extension

#
# 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) = "Interaction Energy";
$calcSheet->ColumnHeading(7) = "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=1; $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"); 

     # 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("Layer 2")->Atoms->Delete;   

     $layer2Doc->CopyFrom($doc);
     $layer2Doc->DisplayRange->Sets("Layer 1")->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;

     #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;

     # 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 $interactionEnergy = $totalEnergy - ($layer1Energy + $layer2Energy);
     $calcSheet->Cell($count-1, 6) = $interactionEnergy;
     
     my $interactionEnergyArea = $interactionEnergy / $surfaceArea;
     $calcSheet->Cell($count-1, 7) = $interactionEnergyArea;
          
     # Discard the temporary documents
     
     $allDoc->Discard;
     $layer1Doc->Discard;
     $layer2Doc->Discard;
}
print "Calculation complete.\n";