#!perl

 

use strict;
use Getopt::Long;
use MaterialsScript qw(:all);

 

#################################################################################
#
# Authors: Stephen Todd, Jason DeJoannis
# Date:  Apr. 22, 2013
# Purpose: Interaction energies between each pair of sets
# Requirements: Trajectory with sets defined
#  The sets must be non-overlapping (no two can contain the same atom or bead)
#  A Forcite/Mesocite settings file with the same name as the trajectory
#
# Works for Forcite or Mesocite and with standard MS or OFF forcefields
#
# Outputs: A study table with a summary page and one page for each interaction
#   van der Waals and Coulomb energy components
#   Hydrogen bond counts
#
#################################################################################
# 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{"test.xtd"};
#my \$forcefield = "";
#my \$firstframe = 1;
#my \$lastframe  = 0;
#my \$structFreq = 10;

 

#
# 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++)
  {
 
   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;
}