Dear All
I'm using the following script to calculate the interaction energy between an ion and a solid phase. This script works between a monovalent metallic ion and a solid phase, however, if I changed the ion to a sulfate, the solid phase is the same, I got the following error:
Unable to find settings file Layer (function/property "LoadSettings") at -e line 66.
I attached the script and highlighted line 66. I appreciated any suggestion.
#!perl
use strict;
use Getopt::Long;
use MaterialsScript qw(:all);
# Begin user settings
# These are for the user menu
#my \\\$doc = Documents->ActiveDocument;
#my %Arg;
#GetOptions(\\%Arg, "Forcefield_file=s", "First_frame=i", "Last_frame=i",
# "Write_structure_frequency=i");
#my \\\$forcefield = \\\$Arg{"Forcefield_file"}; # Leave blank for standard MS ff specified via settings file
#my \\\$firstframe = \\\$Arg{"First_frame"}; # Frame to start from
#my \\\$lastframe = \\\$Arg{"Last_frame"}; # Frame to end on, 0=last
#my \\\$structFreq = \\\$Arg{"Write_structure_frequency"}; # Frequency for writing structure to table, 0=none
# Comment the above lines and uncomment these to run from script editor
my \\\$doc = \\\$Documents{"Layer.xtd"};
my \\\$forcefield = "";
my \\\$firstframe = 80000;
my \\\$lastframe = 81000;
my \\\$structFreq = 0;
#
# End user settings
#################################################################################
# The trajectory document
my \\\$filename = \\\$doc->Name;
my \\\$doc = \\\$Documents{"\\\$filename.xtd"};
# Automatically detect Forcite or Mesocite
my \\\$module = Modules->Forcite;
my \\\$ObjectType = "Atom";
if (\\\$doc->UnitCell->Beads->Count > \\\$doc->UnitCell->Atoms->Count)
{
\\\$module = Modules->Mesocite;
\\\$ObjectType = "Bead";
}
my \\\$ObjectTypes = \\\$ObjectType . "s";
# Check sets are ok
die "The sets must be mutually exclusive and non-empty\\n" unless SetsAreOK(\\\$doc, \\\$ObjectType);
# The energy settings
\\\$module->LoadSettings(\\\$filename);
if (\\\$forcefield ne "")
{
\\\$forcefield =~ s/\\.off\\\$//;
\\\$module->ChangeSettings([CurrentForcefield => "/\\\$forcefield"]);
print "Using OFF forcefield \\\$forcefield\\n";
}
# Create a new study table to hold the structures and energies and set the headings
my \\\$newStudyTable = Documents->New("\\\$filename"."_IntEnergy.std");
my \\\$avgSheet = \\\$newStudyTable->Sheets->Item(0);
\\\$avgSheet->Title = "Statistics";
# Create an array with all the set names
my \\\$maxsets = 4; # adjustable
my \\\$nsets = \\\$doc->UnitCell->Sets->Count;
if (\\\$nsets > \\\$maxsets) { \\\$nsets = \\\$maxsets; }
my @setname;
for (my \\\$i=0; \\\$i<\\\$nsets; \\\$i++)
{
push @setname, \\\$doc->UnitCell->Sets(\\\$i)->Name;
}
# Loop over every pair of sets
for (my \\\$i=0; \\\$i<@setname; \\\$i++)
{
my \\\$setname1 = \\\$setname[\\\$i];
for (my \\\$j=\\\$i+1; \\\$j<@setname; \\\$j++)
{
my \\\$setname2 = \\\$setname[\\\$j];
printf "Interactions between sets %s and %s\\n", \\\$setname1, \\\$setname2;
my \\\$calcSheet = \\\$newStudyTable->InsertSheet;
\\\$calcSheet->Title = "\\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(0) = "Frame";
\\\$calcSheet->ColumnHeading(1) = "Time (ps)";
\\\$calcSheet->ColumnHeading(2) = "\\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(3) = "Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(4) = "vdW Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(5) = "Coulomb Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(6) = "Hbonds \\\$setname1 + \\\$setname2" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(7) = "\\\$setname1";
\\\$calcSheet->ColumnHeading(8) = "Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(9) = "vdW Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(10) = "Coulomb Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(11) = "Hbonds \\\$setname1" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(12) = "\\\$setname2";
\\\$calcSheet->ColumnHeading(13) = "Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(14) = "vdW Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(15) = "Coulomb Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(16) = "Hbonds \\\$setname2" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(17) = "Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(18) = "vdW Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(19) = "Coulomb Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(20) = "Hbonds between \\\$setname1 \\\$setname2" if (\\\$ObjectType eq "Atom");
# Get the total number of frames in the trajectory
my \\\$numFrames = \\\$doc->Trajectory->NumFrames;
if (\\\$lastframe == 0) { \\\$lastframe = \\\$numFrames; }
print "Calculating interaction energy from frame \\\$firstframe to \\\$lastframe on \\\$filename trajectory\\n";
# Starts to loop over the frames in the trajectory
for ( my \\\$count=\\\$firstframe; \\\$count<=\\\$lastframe; \\\$count=\\\$count+10)
{
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 \\\$bothDoc = Documents->New("both.xsd");
my \\\$set1Doc = Documents->New("set1.xsd");
my \\\$set2Doc = Documents->New("set2.xsd");
# Create a copy of the structure in temporary documents and delete
# the atoms or beads we don't need. This way bonds are preserved.
\\\$bothDoc->CopyFrom(\\\$doc);
my \\\$cell = \\\$bothDoc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$i)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$cell->Sets(\\\$j)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
\\\$set1Doc->CopyFrom(\\\$doc);
my \\\$cell = \\\$set1Doc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$i)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
\\\$set2Doc->CopyFrom(\\\$doc);
my \\\$cell = \\\$set2Doc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$j)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
# Calculate the energies for each grouping
\\\$module->Energy->Run(\\\$bothDoc);
\\\$module->Energy->Run(\\\$set1Doc);
\\\$module->Energy->Run(\\\$set2Doc);
# Calculate hydrogen bonds
if (\\\$ObjectType eq "Atom")
{
\\\$bothDoc->CalculateHBonds;
\\\$set1Doc->CalculateHBonds;
\\\$set2Doc->CalculateHBonds;
}
# Frame number
my \\\$row = \\\$count - \\\$firstframe;
\\\$calcSheet->Cell(\\\$row, 0) = \\\$doc->Trajectorys->Item(0)->CurrentFrame;
\\\$calcSheet->Cell(\\\$row, 1) = \\\$doc->Trajectorys->Item(0)->FrameTime;
# Put the structures in a study table for safekeeping
if (\\\$structFreq > 0 and \\\$row % \\\$structFreq == 0)
{
\\\$calcSheet->Cell(\\\$row, 2) = \\\$bothDoc;
\\\$calcSheet->Cell(\\\$row, 7) = \\\$set1Doc;
\\\$calcSheet->Cell(\\\$row, 12) = \\\$set2Doc;
}
# Individual energies
\\\$calcSheet->Cell(\\\$row, 3) = \\\$bothDoc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 4) = \\\$bothDoc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 5) = \\\$bothDoc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 6) = \\\$bothDoc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
\\\$calcSheet->Cell(\\\$row, 8) = \\\$set1Doc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 9) = \\\$set1Doc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 10) = \\\$set1Doc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 11) = \\\$set1Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
\\\$calcSheet->Cell(\\\$row, 13) = \\\$set2Doc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 14) = \\\$set2Doc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 15) = \\\$set2Doc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 16) = \\\$set2Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
# Interaction energies
my \\\$col=17;
foreach my \\\$component ("PotentialEnergy", "vanderWaalsEnergy", "ElectrostaticEnergy")
{
my \\\$interaction = \\\$bothDoc->\\\$component - \\\$set1Doc->\\\$component - \\\$set2Doc->\\\$component;
\\\$calcSheet->Cell(\\\$row, \\\$col) = \\\$interaction;
\\\$col++;
}
\\\$calcSheet->Cell(\\\$row, 20) = \\\$bothDoc->UnitCell->HydrogenBonds->Count
- \\\$set1Doc->UnitCell->HydrogenBonds->Count
- \\\$set2Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
# Discard the temporary documents
\\\$bothDoc->Discard;
\\\$set1Doc->Discard;
\\\$set2Doc->Discard;
} #next frame
}# next j
}# next i
# Calculate the averages
\\\$avgSheet->ColumnHeading(0) = "Interaction";
\\\$avgSheet->ColumnHeading(1) = "Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(2) = "n";
\\\$avgSheet->ColumnHeading(3) = "StdDev";
\\\$avgSheet->ColumnHeading(4) = "StdErr";
\\\$avgSheet->ColumnHeading(5) = "vdW Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(6) = "n";
\\\$avgSheet->ColumnHeading(7) = "StdDev";
\\\$avgSheet->ColumnHeading(8) = "StdErr";
\\\$avgSheet->ColumnHeading(9) = "Coulomb Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(10) = "n";
\\\$avgSheet->ColumnHeading(11) = "StdDev";
\\\$avgSheet->ColumnHeading(12) = "StdErr";
if (\\\$ObjectType eq "Atom")
{
\\\$avgSheet->ColumnHeading(13) = "Hbonds";
\\\$avgSheet->ColumnHeading(14) = "n";
\\\$avgSheet->ColumnHeading(15) = "StdDev";
\\\$avgSheet->ColumnHeading(16) = "StdErr";
}
my \\\$col_a = 17;
my \\\$col_b = 20;
\\\$col_b-- unless (\\\$ObjectType eq "Atom");
for (my \\\$pair=1; \\\$pair<\\\$newStudyTable->Sheets->Count; \\\$pair++)
{
my \\\$sheet = \\\$newStudyTable->Sheets(\\\$pair);
\\\$avgSheet->Cell(\\\$pair-1, 0) = \\\$sheet->Title;
for (my \\\$col=\\\$col_a; \\\$col<=\\\$col_b; \\\$col++)
{
my \\\$n=0; my \\\$avg=0; my \\\$msq=0;
for (my \\\$i=0; \\\$i<\\\$sheet->RowCount; \\\$i++)
{
my \\\$value;
eval { \\\$value = \\\$sheet->Cell(\\\$i, \\\$col); };
last if (\\\$@);
last if (\\\$value eq "");
\\\$avg+=\\\$value;
\\\$msq+=(\\\$value*\\\$value);
\\\$n++;
}
\\\$avg/=\\\$n;
\\\$msq/=\\\$n;
my \\\$var = \\\$msq - \\\$avg*\\\$avg;
my \\\$colout = 4*(\\\$col-\\\$col_a)+1;
\\\$avgSheet->Cell(\\\$pair-1, \\\$colout) = \\\$avg;
\\\$avgSheet->Cell(\\\$pair-1, 1+\\\$colout) = \\\$n;
\\\$avgSheet->Cell(\\\$pair-1, 2+\\\$colout) = sqrt(\\\$var);
\\\$avgSheet->Cell(\\\$pair-1, 3+\\\$colout) = sqrt(\\\$var/\\\$n);
}
}
print "Calculation complete!\\n";
sub SetsAreOK
{
my (\\\$doc, \\\$ObjectType) = @_;
my \\\$objects = \\\$ObjectType ."s";
my \\\$items = \\\$doc->UnitCell->\\\$objects;
my \\\$sets = \\\$doc->UnitCell->Sets;
# Get the total number of atoms or beads
my \\\$total = \\\$items->Count;
# Cache their names
my @name;
foreach my \\\$item (@{\\\$items})
{
push @name, \\\$item->Name;
}
# Go through each set and each item and rename them
my \\\$tot_in_sets=0;
foreach my \\\$set (@{\\\$sets})
{
my \\\$items = \\\$set->\\\$objects;
my \\\$n = \\\$items->Count;
if (\\\$n == 0)
{
print "Set ".\\\$set->Name." is empty (UnitCell)\\n";
return 0;
}
\\\$tot_in_sets += \\\$n;
foreach my \\\$item (@{\\\$items})
{
if (\\\$item->Name =~ /^Member of set /)
{
my \\\$otherset = \\\$item->Name;
\\\$otherset =~ s/^Member of set //;
print "Overlap between sets ".\\\$set->Name." and \\\$otherset\\n";
return 0;
};
\\\$item->Name = "Member of set ".\\\$set->Name;
}
}
# Put their names back
my \\\$i=0;
foreach my \\\$item (@{\\\$items})
{
\\\$item->Name = \\\$name[\\\$i];
\\\$i++;
}
if (\\\$total < \\\$tot_in_sets)
{
print "More \\\$objects in sets than in document\\n";
return 0;
}
return 1;
}
Thank you
Jingjing
I'm using the following script to calculate the interaction energy between an ion and a solid phase. This script works between a monovalent metallic ion and a solid phase, however, if I changed the ion to a sulfate, the solid phase is the same, I got the following error:
Unable to find settings file Layer (function/property "LoadSettings") at -e line 66.
I attached the script and highlighted line 66. I appreciated any suggestion.
#!perl
use strict;
use Getopt::Long;
use MaterialsScript qw(:all);
# Begin user settings
# These are for the user menu
#my \\\$doc = Documents->ActiveDocument;
#my %Arg;
#GetOptions(\\%Arg, "Forcefield_file=s", "First_frame=i", "Last_frame=i",
# "Write_structure_frequency=i");
#my \\\$forcefield = \\\$Arg{"Forcefield_file"}; # Leave blank for standard MS ff specified via settings file
#my \\\$firstframe = \\\$Arg{"First_frame"}; # Frame to start from
#my \\\$lastframe = \\\$Arg{"Last_frame"}; # Frame to end on, 0=last
#my \\\$structFreq = \\\$Arg{"Write_structure_frequency"}; # Frequency for writing structure to table, 0=none
# Comment the above lines and uncomment these to run from script editor
my \\\$doc = \\\$Documents{"Layer.xtd"};
my \\\$forcefield = "";
my \\\$firstframe = 80000;
my \\\$lastframe = 81000;
my \\\$structFreq = 0;
#
# End user settings
#################################################################################
# The trajectory document
my \\\$filename = \\\$doc->Name;
my \\\$doc = \\\$Documents{"\\\$filename.xtd"};
# Automatically detect Forcite or Mesocite
my \\\$module = Modules->Forcite;
my \\\$ObjectType = "Atom";
if (\\\$doc->UnitCell->Beads->Count > \\\$doc->UnitCell->Atoms->Count)
{
\\\$module = Modules->Mesocite;
\\\$ObjectType = "Bead";
}
my \\\$ObjectTypes = \\\$ObjectType . "s";
# Check sets are ok
die "The sets must be mutually exclusive and non-empty\\n" unless SetsAreOK(\\\$doc, \\\$ObjectType);
# The energy settings
\\\$module->LoadSettings(\\\$filename);
if (\\\$forcefield ne "")
{
\\\$forcefield =~ s/\\.off\\\$//;
\\\$module->ChangeSettings([CurrentForcefield => "/\\\$forcefield"]);
print "Using OFF forcefield \\\$forcefield\\n";
}
# Create a new study table to hold the structures and energies and set the headings
my \\\$newStudyTable = Documents->New("\\\$filename"."_IntEnergy.std");
my \\\$avgSheet = \\\$newStudyTable->Sheets->Item(0);
\\\$avgSheet->Title = "Statistics";
# Create an array with all the set names
my \\\$maxsets = 4; # adjustable
my \\\$nsets = \\\$doc->UnitCell->Sets->Count;
if (\\\$nsets > \\\$maxsets) { \\\$nsets = \\\$maxsets; }
my @setname;
for (my \\\$i=0; \\\$i<\\\$nsets; \\\$i++)
{
push @setname, \\\$doc->UnitCell->Sets(\\\$i)->Name;
}
# Loop over every pair of sets
for (my \\\$i=0; \\\$i<@setname; \\\$i++)
{
my \\\$setname1 = \\\$setname[\\\$i];
for (my \\\$j=\\\$i+1; \\\$j<@setname; \\\$j++)
{
my \\\$setname2 = \\\$setname[\\\$j];
printf "Interactions between sets %s and %s\\n", \\\$setname1, \\\$setname2;
my \\\$calcSheet = \\\$newStudyTable->InsertSheet;
\\\$calcSheet->Title = "\\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(0) = "Frame";
\\\$calcSheet->ColumnHeading(1) = "Time (ps)";
\\\$calcSheet->ColumnHeading(2) = "\\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(3) = "Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(4) = "vdW Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(5) = "Coulomb Energy of \\\$setname1 + \\\$setname2";
\\\$calcSheet->ColumnHeading(6) = "Hbonds \\\$setname1 + \\\$setname2" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(7) = "\\\$setname1";
\\\$calcSheet->ColumnHeading(8) = "Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(9) = "vdW Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(10) = "Coulomb Energy of \\\$setname1";
\\\$calcSheet->ColumnHeading(11) = "Hbonds \\\$setname1" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(12) = "\\\$setname2";
\\\$calcSheet->ColumnHeading(13) = "Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(14) = "vdW Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(15) = "Coulomb Energy of \\\$setname2";
\\\$calcSheet->ColumnHeading(16) = "Hbonds \\\$setname2" if (\\\$ObjectType eq "Atom");
\\\$calcSheet->ColumnHeading(17) = "Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(18) = "vdW Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(19) = "Coulomb Interaction Energy \\\$setname1 \\\$setname2";
\\\$calcSheet->ColumnHeading(20) = "Hbonds between \\\$setname1 \\\$setname2" if (\\\$ObjectType eq "Atom");
# Get the total number of frames in the trajectory
my \\\$numFrames = \\\$doc->Trajectory->NumFrames;
if (\\\$lastframe == 0) { \\\$lastframe = \\\$numFrames; }
print "Calculating interaction energy from frame \\\$firstframe to \\\$lastframe on \\\$filename trajectory\\n";
# Starts to loop over the frames in the trajectory
for ( my \\\$count=\\\$firstframe; \\\$count<=\\\$lastframe; \\\$count=\\\$count+10)
{
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 \\\$bothDoc = Documents->New("both.xsd");
my \\\$set1Doc = Documents->New("set1.xsd");
my \\\$set2Doc = Documents->New("set2.xsd");
# Create a copy of the structure in temporary documents and delete
# the atoms or beads we don't need. This way bonds are preserved.
\\\$bothDoc->CopyFrom(\\\$doc);
my \\\$cell = \\\$bothDoc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$i)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$cell->Sets(\\\$j)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
\\\$set1Doc->CopyFrom(\\\$doc);
my \\\$cell = \\\$set1Doc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$i)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
\\\$set2Doc->CopyFrom(\\\$doc);
my \\\$cell = \\\$set2Doc->UnitCell;
my \\\$objects = \\\$cell->\\\$ObjectTypes;
foreach (@{\\\$objects}) { \\\$_->Name = "Delete"; }
foreach (@{\\\$cell->Sets(\\\$j)->\\\$ObjectTypes}) { \\\$_->Name = "Keep"; }
foreach (@{\\\$objects}) { \\\$_->Delete if (\\\$_->Name eq "Delete"); }
# Calculate the energies for each grouping
\\\$module->Energy->Run(\\\$bothDoc);
\\\$module->Energy->Run(\\\$set1Doc);
\\\$module->Energy->Run(\\\$set2Doc);
# Calculate hydrogen bonds
if (\\\$ObjectType eq "Atom")
{
\\\$bothDoc->CalculateHBonds;
\\\$set1Doc->CalculateHBonds;
\\\$set2Doc->CalculateHBonds;
}
# Frame number
my \\\$row = \\\$count - \\\$firstframe;
\\\$calcSheet->Cell(\\\$row, 0) = \\\$doc->Trajectorys->Item(0)->CurrentFrame;
\\\$calcSheet->Cell(\\\$row, 1) = \\\$doc->Trajectorys->Item(0)->FrameTime;
# Put the structures in a study table for safekeeping
if (\\\$structFreq > 0 and \\\$row % \\\$structFreq == 0)
{
\\\$calcSheet->Cell(\\\$row, 2) = \\\$bothDoc;
\\\$calcSheet->Cell(\\\$row, 7) = \\\$set1Doc;
\\\$calcSheet->Cell(\\\$row, 12) = \\\$set2Doc;
}
# Individual energies
\\\$calcSheet->Cell(\\\$row, 3) = \\\$bothDoc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 4) = \\\$bothDoc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 5) = \\\$bothDoc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 6) = \\\$bothDoc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
\\\$calcSheet->Cell(\\\$row, 8) = \\\$set1Doc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 9) = \\\$set1Doc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 10) = \\\$set1Doc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 11) = \\\$set1Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
\\\$calcSheet->Cell(\\\$row, 13) = \\\$set2Doc->PotentialEnergy;
\\\$calcSheet->Cell(\\\$row, 14) = \\\$set2Doc->vanderWaalsEnergy;
\\\$calcSheet->Cell(\\\$row, 15) = \\\$set2Doc->ElectrostaticEnergy;
\\\$calcSheet->Cell(\\\$row, 16) = \\\$set2Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
# Interaction energies
my \\\$col=17;
foreach my \\\$component ("PotentialEnergy", "vanderWaalsEnergy", "ElectrostaticEnergy")
{
my \\\$interaction = \\\$bothDoc->\\\$component - \\\$set1Doc->\\\$component - \\\$set2Doc->\\\$component;
\\\$calcSheet->Cell(\\\$row, \\\$col) = \\\$interaction;
\\\$col++;
}
\\\$calcSheet->Cell(\\\$row, 20) = \\\$bothDoc->UnitCell->HydrogenBonds->Count
- \\\$set1Doc->UnitCell->HydrogenBonds->Count
- \\\$set2Doc->UnitCell->HydrogenBonds->Count if (\\\$ObjectType eq "Atom");
# Discard the temporary documents
\\\$bothDoc->Discard;
\\\$set1Doc->Discard;
\\\$set2Doc->Discard;
} #next frame
}# next j
}# next i
# Calculate the averages
\\\$avgSheet->ColumnHeading(0) = "Interaction";
\\\$avgSheet->ColumnHeading(1) = "Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(2) = "n";
\\\$avgSheet->ColumnHeading(3) = "StdDev";
\\\$avgSheet->ColumnHeading(4) = "StdErr";
\\\$avgSheet->ColumnHeading(5) = "vdW Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(6) = "n";
\\\$avgSheet->ColumnHeading(7) = "StdDev";
\\\$avgSheet->ColumnHeading(8) = "StdErr";
\\\$avgSheet->ColumnHeading(9) = "Coulomb Interaction Energy (kcal/mol)";
\\\$avgSheet->ColumnHeading(10) = "n";
\\\$avgSheet->ColumnHeading(11) = "StdDev";
\\\$avgSheet->ColumnHeading(12) = "StdErr";
if (\\\$ObjectType eq "Atom")
{
\\\$avgSheet->ColumnHeading(13) = "Hbonds";
\\\$avgSheet->ColumnHeading(14) = "n";
\\\$avgSheet->ColumnHeading(15) = "StdDev";
\\\$avgSheet->ColumnHeading(16) = "StdErr";
}
my \\\$col_a = 17;
my \\\$col_b = 20;
\\\$col_b-- unless (\\\$ObjectType eq "Atom");
for (my \\\$pair=1; \\\$pair<\\\$newStudyTable->Sheets->Count; \\\$pair++)
{
my \\\$sheet = \\\$newStudyTable->Sheets(\\\$pair);
\\\$avgSheet->Cell(\\\$pair-1, 0) = \\\$sheet->Title;
for (my \\\$col=\\\$col_a; \\\$col<=\\\$col_b; \\\$col++)
{
my \\\$n=0; my \\\$avg=0; my \\\$msq=0;
for (my \\\$i=0; \\\$i<\\\$sheet->RowCount; \\\$i++)
{
my \\\$value;
eval { \\\$value = \\\$sheet->Cell(\\\$i, \\\$col); };
last if (\\\$@);
last if (\\\$value eq "");
\\\$avg+=\\\$value;
\\\$msq+=(\\\$value*\\\$value);
\\\$n++;
}
\\\$avg/=\\\$n;
\\\$msq/=\\\$n;
my \\\$var = \\\$msq - \\\$avg*\\\$avg;
my \\\$colout = 4*(\\\$col-\\\$col_a)+1;
\\\$avgSheet->Cell(\\\$pair-1, \\\$colout) = \\\$avg;
\\\$avgSheet->Cell(\\\$pair-1, 1+\\\$colout) = \\\$n;
\\\$avgSheet->Cell(\\\$pair-1, 2+\\\$colout) = sqrt(\\\$var);
\\\$avgSheet->Cell(\\\$pair-1, 3+\\\$colout) = sqrt(\\\$var/\\\$n);
}
}
print "Calculation complete!\\n";
sub SetsAreOK
{
my (\\\$doc, \\\$ObjectType) = @_;
my \\\$objects = \\\$ObjectType ."s";
my \\\$items = \\\$doc->UnitCell->\\\$objects;
my \\\$sets = \\\$doc->UnitCell->Sets;
# Get the total number of atoms or beads
my \\\$total = \\\$items->Count;
# Cache their names
my @name;
foreach my \\\$item (@{\\\$items})
{
push @name, \\\$item->Name;
}
# Go through each set and each item and rename them
my \\\$tot_in_sets=0;
foreach my \\\$set (@{\\\$sets})
{
my \\\$items = \\\$set->\\\$objects;
my \\\$n = \\\$items->Count;
if (\\\$n == 0)
{
print "Set ".\\\$set->Name." is empty (UnitCell)\\n";
return 0;
}
\\\$tot_in_sets += \\\$n;
foreach my \\\$item (@{\\\$items})
{
if (\\\$item->Name =~ /^Member of set /)
{
my \\\$otherset = \\\$item->Name;
\\\$otherset =~ s/^Member of set //;
print "Overlap between sets ".\\\$set->Name." and \\\$otherset\\n";
return 0;
};
\\\$item->Name = "Member of set ".\\\$set->Name;
}
}
# Put their names back
my \\\$i=0;
foreach my \\\$item (@{\\\$items})
{
\\\$item->Name = \\\$name[\\\$i];
\\\$i++;
}
if (\\\$total < \\\$tot_in_sets)
{
print "More \\\$objects in sets than in document\\n";
return 0;
}
return 1;
}
Thank you
Jingjing
