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