Welcome To Support Community

Materials Studio

Advanced Search
Ask Search:
jjbjjb 

Error for interaction energy calculation

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