Dear Reinier,
I read your post about "How to measure dihedral angles between two vectors" (Updated Link)
And I modified the following script according to your recommendation in above post. However, the angle is still between 60-120 degree as Jingjing mentioned before. If it's available for you, would you please update the correct script? Following image is the perfect result I am looking for.
----------------------------------
#!perl
use strict;
use MaterialsScript qw(:all);
use Math::Trig;
my \\\$doc = \\\$Documents{"Layer.xtd"};
my \\\$setA = "Si";
my \\\$setB = "H2O";
# 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 = 1; \\\$frame <= 1000; \\\$frame++)
{
\\\$trj->CurrentFrame = \\\$frame;
# get the normal vector of the plane
my \\\$normal = \\\$plane->BestFitPlane->UnitNormalVector;
foreach my \\\$setMolecule(@\\\$setMolecules)
{
# z position of the molecule
my \\\$distance = \\\$setMolecule->Center->Z;
# get dipole moment
my \\\$dipoleMoment = \\\$setMolecule->DipoleMoment;
# calculate the angle between the dipole and the surface normal
my \\\$cosAngle = InProduct(\\\$normal, \\\$dipoleMoment);
my \\\$angle =180/pi*acos(\\\$cosAngle);
# store the data in the study table
\\\$std->Cell(\\\$irow,0) = \\\$angle;
\\\$std->Cell(\\\$irow++,1) = \\\$distance;
}
}
\\\$doc->Discard; # discard any changes made
sub InProduct(){
my(\\\$a, \\\$b) = @_;
return \\\$a->X*\\\$b->X+\\\$a->Y*\\\$b->Y+\\\$a->Z*\\\$b->Z
}
Regards,
Jingjing
reinier
The problem is probably in your InProduct function;it has a prototype () indicating the function takes no parameters, whereas it takes two,
sub InProduct(){
my(\\\$a, \\\$b) = @_;
return \\\$a->X*\\\$b->X+\\\$a->Y*\\\$b->Y+\\\$a->Z*\\\$b->Z
}
Deleting () should resolve the problem.
reinier
You need to normalize the dipole moment vector before taking the inproduct, as above,
# get the direction of the dipole
my \\\$dipoleMoment = \\\$molecule->DipoleMoment;
my \\\$direction = Normalize(\\\$dipoleMoment);
# calculate the angle between the dipole and the z-axis
my \\\$angle = 180/pi*acos(InProduct(\\\$direction, \\\$axis));
------------------------------------
