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