Dear alll,

I used the following code to analysis water dipole angle from the surface. This code works for some system, but not work for another system. I don't know why. It seems the "foreach my \\\$setMolecule(@\\\$setMolecules)" doesn't work for specific system. If it's available for you, would please help me out?

I insert "print "I am here\\n";" in the foreach cycle, it doesn't show any result in output.
-----------------------------------------
 #!perl
use strict;
use MaterialsScript qw(:all);
use Math::Trig;

my \\\$doc = \\\$Documents{"12354687892.xtd"};
my \\\$setA = "ai";
my \\\$setB = "water";

# get the trajectory data
my \\\$trj = \\\$doc->Trajectory;
my \\\$numFrames = \\\$trj->NumFrames;

# create a study table for the results
my \\\$std = Documents->New(\\\$doc->Name.".std");
\\\$std->ColumnHeading(0) = "Angle (deg)";
\\\$std->ColumnHeading(1) = "Distance (angstrom)";
my \\\$setMolecules = \\\$doc->UnitCell->Sets(\\\$setB)->Molecules;

# run over all frames in the trajectory
my \\\$irow = 0;
my \\\$setAtoms = \\\$doc->UnitCell->Sets(\\\$setA)->Atoms;
my \\\$plane = \\\$doc->UnitCell->CreateBestFitPlane(\\\$setAtoms);
my \\\$bestFitPlane = \\\$plane->BestFitPlane;
for(my \\\$frame = 201; \\\$frame <= 300; \\\$frame++)
{
   \\\$trj->CurrentFrame = \\\$frame;

   # get the normal vector of the plane
   my \\\$normal = \\\$plane->BestFitPlane->UnitNormalVector;
   
   foreach my \\\$setMolecule(@\\\$setMolecules)
   { 
     print "I am here\\n"; 
     # z position of the molecule
     my \\\$distance = \\\$setMolecule->Center->Z;

     # get dipole moment
     my \\\$dipoleMoment = \\\$setMolecule->DipoleMoment;
     my \\\$direction=Normalize(\\\$dipoleMoment);

     # calculate the angle between the dipole and the surface normal  
     my \\\$cosAngle = InProduct(\\\$normal, \\\$direction);
     my \\\$angle =180/pi*acos(\\\$cosAngle);
     
     print "\\\$distance\\n";
     print "\\\$angle\\n";
       
     # store the data in the study table
     \\\$std->Cell(\\\$irow,1) = \\\$distance;       
     \\\$std->Cell(\\\$irow++,0) = \\\$angle;
   }
}
\\\$doc->Discard; # discard any changes made

sub InProduct(){
my(\\\$a, \\\$b) = @_;
return \\\$a->X*\\\$b->X+\\\$a->Y*\\\$b->Y+\\\$a->Z*\\\$b->Z
}

sub Normalize
{
    my \\\$p = shift;
    my \\\$length = sqrt(\\\$p->X*\\\$p->X +    \\\$p->Y*\\\$p->Y + \\\$p->Z*\\\$p->Z);
    \\\$p->X /= \\\$length;
    \\\$p->Y /= \\\$length;
    \\\$p->Z /= \\\$length;
    return \\\$p;
}
-----------------------------------------------------------------------