Welcome To Support Community

Materials Studio

Advanced Search
Ask Search:

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.


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
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.

                 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"); }

                 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"); }

                 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

            # Calculate hydrogen bonds
            if ($ObjectType eq "Atom")

            # 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;
            $calcSheet->Cell($row, 20) = $bothDoc->UnitCell->HydrogenBonds->Count
                           - $set1Doc->UnitCell->HydrogenBonds->Count
                           - $set2Doc->UnitCell->HydrogenBonds->Count if ($ObjectType eq "Atom");
                 # Discard the temporary documents

             } #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 "");
        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];
    if ($total < $tot_in_sets)
        print "More $objects in sets than in document\n";
        return 0;
    return 1;

Thank you 
Jason DeJoannisJason DeJoannis
Hi Jingjing,

The thread for that script (which I wrote) is linked below and includes what might be an answer to your question:
Interaction energies between more than two groups

Basically it is looking for a Materials Studio settings file with the name "Layer".
Hi Jason
Thanks for your reply. The setting problem was solved. However, I had another problem like another person asked earlier. If I chose different amount of solute atoms as one set, I'll get different interaction energies. I did use the latest version you posted on April 24, 2013, and the date on the script was April 22, 2013 (new the third one you have posted). My system contains a layer of liquid solution, and a layer of solid phase. I have two sets: solid layer and 1 solute atom. I also noticed that if I chose part of the solid layer, I also got different interaction energies. Could you please help me out of this? Thank you very much!

All the best,
Jason DeJoannisJason DeJoannis
Hi Jingjing,
The script is designed to return the interaction energies between the user-defined sets of atoms. Thus if the user has specified two sets A and B, the script ignores all other atoms in the system. So in your case, it sounds like it is working like expected. I do not see a problem. When you change the members of either A or B, you should indeed expect different energies. What were you expecting? You should probably normalize the interaction energy by surface area.