Dear all,
I modified Stephen Todd's CoordinationNumberDistance script to calculate the evolution of coordination number for a trajectory. I didn't integrate the RDF curves because I'm calculating the CN for the surface species, the coordinating ions were not homogeneous distributed. I looped over the frames, created temporary structure files to calculate the coordination number for each frame, then tried to write the coordination number of each frame into a study table. However, this script didn't calculate the coordination number properly. It added the coordination number of all the frames, and put the sum into the cells where it was supposed to be the coordination number of each frame. Here is my script:
use strict;
use MaterialsScript qw(:all);
######################################################################################
# Begin User Editable Settings
my \\\$doc = \\\$Documents{"new-140000-200001.xtd"}; # Set this to the start document
my \\\$coordRadius =4.43; # The coordination distance radius defining the coordination sphere
my \\\$coordinatedAtom = "Ca1"; # The coordinated atom
my \\\$coordinatingAtom = "O*"; # The coordinating atoms
my \\\$maxAtoms = 1000; # Maximum number of atoms for the structure to be included
#
# End User Editable Settings
########################################################################################
# Sets the working document
# Sets the lattice 3D orientation
\\\$doc->Lattice3D->OrientationConvention = "A along X, B in XY plane";
# Check the number of atoms and sets variables based on them
my \\\$numAtoms = \\\$doc->SymmetrySystem->NumberOfAtoms;
print "The number of atoms in the system is \\\$numAtoms.\\n";
my \\\$addDoc = 0; # Variable used to determine whether the add structure to the study table
my \\\$arrayTerms = 3; # If you are adding the structure, the number of terms in the array changes
# from three to four (three coordinates + atom)
# Sets the values of addDoc and arrayTerms depending on the number of atoms in the system
if (\\\$numAtoms > \\\$maxAtoms) {
\\\$addDoc = 0;
\\\$arrayTerms = 3;
print "The number of atoms exceeds the maxAtoms, the document will not be added to the study table.\\n";
} else {
\\\$addDoc = 1;
\\\$arrayTerms = 4;
print "The number of atoms is less than the maxAtoms, the document will be added to the study table.\\n";
}
# Set up the study table and define the number of columns depending on whether a structure
# column is going to be added. This uses the createColumns subroutine.
my \\\$newStudyTable = Documents->New("CoordRadius\\\$coordRadius.std");
my \\\$calcSheet = \\\$newStudyTable->Sheets->Item(0);
if (\\\$addDoc == 1) {
createColumns(\\\$calcSheet, 1, 2, 3, 4);
\\\$calcSheet->ColumnHeading(0) = "Cell";
print "Created study table with column for cell.\\n";
} elsif (\\\$addDoc == 0) {
createColumns(\\\$calcSheet, 0, 1, 2, 3);
print "Created study table without column for cell.\\n";
} else {
die "There seems to be a problem with the value of addDoc";
}
# Gets the lengths of the axis of the system for dealing with periodic boundary conditions later
my \\\$lengthX = \\\$doc->Lattice3D->LengthA;
my \\\$lengthY = \\\$doc->Lattice3D->LengthB;
my \\\$lengthZ = \\\$doc->Lattice3D->LengthC;
# Gets all the atoms into a collection
my \\\$atoms = \\\$doc->UnitCell->Atoms;
# Finds the display style on the first atom
my \\\$displayStyle = \\\$doc->UnitCell->Atoms(0)->DisplayStyle;
############################################################################################################
# Go through all the atoms and assigns them to either coordinating or coordinated arrays depending on their
# element symbol. You can change this to forcefield type or whatever you wish.
my @coordinatingAtomsXYZ; # Array to hold coordinating atoms coordinates
my @coordinatedAtomsXYZ; # Array to hold coordinated atoms coordinates
my \\\$rowCounter = 0; # Simple row counter for writing to the study table - used later as well
for(my \\\$frame = 40001; \\\$frame <= 60002; \\\$frame = ++\\\$frame) {
#foreach my \\\$atom (@\\\$atoms){
\\\$doc->Trajectory->CurrentFrame = \\\$frame;
my \\\$tmpdoc = Documents->New("tmp.xsd");
\\\$tmpdoc->CopyFrom(\\\$doc);
my \\\$Aatoms = \\\$tmpdoc->UnitCell->Sets("\\\$coordinatedAtom")->Atoms;
my \\\$Batoms = \\\$tmpdoc->UnitCell->Sets("\\\$coordinatingAtom")->Atoms;
foreach my \\\$atom (@\\\$Aatoms) {
#if (\\\$atom->ElementSymbol eq \\\$coordinatedAtom) {
# Writes the coordinates into an array of the form (x1, Y1, Z1, Atom1, X2, Y2... Xn, Yn, Zn, Atomn)
push (@coordinatedAtomsXYZ, \\\$atom->X);
push (@coordinatedAtomsXYZ, \\\$atom->Y);
push (@coordinatedAtomsXYZ, \\\$atom->Z);
# Adds in the link to the atom if this is needed. If you are not adding the structure, this is ignored.
# Also adds the x, y, and z data to the study table, depending on whether the structure is needed.
if (\\\$addDoc == 1) {
push (@coordinatedAtomsXYZ, \\\$atom);
\\\$calcSheet->Cell(\\\$rowCounter, 1) = \\\$atom->XYZ->X;
\\\$calcSheet->Cell(\\\$rowCounter, 2) = \\\$atom->XYZ->Y;
\\\$calcSheet->Cell(\\\$rowCounter, 3) = \\\$atom->XYZ->Z;
} elsif (\\\$addDoc == 0) {
\\\$calcSheet->Cell(\\\$rowCounter, 0) = \\\$atom->XYZ->X;
\\\$calcSheet->Cell(\\\$rowCounter, 1) = \\\$atom->XYZ->Y;
\\\$calcSheet->Cell(\\\$rowCounter, 2) = \\\$atom->XYZ->Z;
}
++\\\$rowCounter; #Increment the row counter for the study table data
# Next part of the loop creates the coordinating atom array. This can be removed or replaced if you
# want to study something different like eg close contacts.
#} elsif (\\\$atom->ElmentSymbol eq \\\$coordinatingAtom) {
}foreach my \\\$atom (@\\\$Batoms){
# Writes the coordinates into an array of the form (x1, Y1, Z1, Atom1, X2, Y2... Xn, Yn, Zn, Atomn)
push (@coordinatingAtomsXYZ, \\\$atom->X);
push (@coordinatingAtomsXYZ, \\\$atom->Y);
push (@coordinatingAtomsXYZ, \\\$atom->Z);
# Adds in the link to the atom id if this is needed. If you are not adding the structure, this is ignored
if (\\\$addDoc == 1 ) {
push (@coordinatingAtomsXYZ, \\\$atom);
}
#}
}
# Calculate the number of atoms in each array. This is dependent on whether the atom id was added to the array
my \\\$numCoordinatedAtoms = @coordinatedAtomsXYZ / \\\$arrayTerms;
my \\\$numCoordinatingAtoms = @coordinatingAtomsXYZ / \\\$arrayTerms;
print "There are \\\$numCoordinatingAtoms coordinating atoms of element \\\$coordinatingAtom.\\n";
print "There are \\\$numCoordinatedAtoms coordinated atoms of element \\\$coordinatedAtom.\\n";
#####################################################################################################
# Go through each coordinated atom and check how many coordinating atoms are within the coordRadius.
my \\\$distanceX = 0; # Define some variables to store the distances in
my \\\$distanceY = 0;
my \\\$distanceZ = 0;
\\\$rowCounter = 0;
for (my \\\$coordinatedAtomCounter = 0; \\\$coordinatedAtomCounter < \\\$numCoordinatedAtoms; ++\\\$coordinatedAtomCounter) {
# Resets the coordination number counter to zero
my \\\$coordinationNumber = 0;
# Defines the working atom list array - this is used to store the atoms whose display style
# you will change. This is not really needed if the addAtoms variable is 1.
my @workingAtomList = ();
# Gets the atom coordinates from the array for the coordinated atoms
my \\\$coordinatedAtomX = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms)];
my \\\$coordinatedAtomY = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +1];
my \\\$coordinatedAtomZ = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +2];
if (\\\$addDoc == 1) {
push (@workingAtomList, @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +3]);
}
# Go through all the coordinating atoms to calculate the coordinatation number
for (my \\\$coordinatingAtomCounter = 0; \\\$coordinatingAtomCounter < \\\$numCoordinatingAtoms; ++\\\$coordinatingAtomCounter) {
# Reads the coordinates from the array which has the form (x1, Y1, Z1, X2, Y2... Xn, Yn, Zn)
my \\\$coordinatingAtomX = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms)];
my \\\$coordinatingAtomY = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms) + 1];
my \\\$coordinatingAtomZ = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms) + 2];
# Use the distanceCalc subroutine to calculate the distance between the atoms
\\\$distanceX = distanceCalc (\\\$coordinatingAtomX,\\\$coordinatedAtomX, \\\$lengthX);
\\\$distanceY = distanceCalc (\\\$coordinatingAtomY,\\\$coordinatedAtomY, \\\$lengthY);
\\\$distanceZ = distanceCalc (\\\$coordinatingAtomZ,\\\$coordinatedAtomZ, \\\$lengthZ);
# Using the different calculations x y, and z distances, calculate the actual vector
my \\\$vector = sqrt (\\\$distanceX*\\\$distanceX + \\\$distanceY*\\\$distanceY + \\\$distanceZ*\\\$distanceZ);
# Check if the vector is less than the coordRadius and update the coordinationNumber if successful
if (\\\$vector <= \\\$coordRadius) {
++\\\$coordinationNumber;
if (\\\$addDoc == 1) {
push (@workingAtomList, @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * 4) + 3]);
}
}
}
# Update the study table with the coordination number and/or document with atoms highlighted
if (\\\$addDoc == 1) {
# Use the working atom list array to change the display style of the coordinating and
# coordinated atoms.
foreach my \\\$workingAtom (@workingAtomList) {
\\\$workingAtom->DisplayStyle = "Ball and stick";
}
# Update the study table
\\\$calcSheet->Cell(\\\$rowCounter, 0) = \\\$tmpdoc;
\\\$calcSheet->Cell(\\\$rowCounter, 4) = \\\$coordinationNumber;
#Reset the display Style on the atoms that have just been changed
foreach my \\\$workingAtom (@workingAtomList){
\\\$workingAtom -> DisplayStyle = \\\$displayStyle;
}
} else {
# Just update the study table
\\\$calcSheet->Cell(\\\$rowCounter, 3) = \\\$coordinationNumber;
}
++\\\$rowCounter; # Updates the row counter
}
\\\$tmpdoc->Discard;
}
print "\\nCalculation Complete \\n";
####################################################################################################
# Subroutine section
# Subroutine to create the columns in the study table, depending on whether a structure is required
sub createColumns {
my (\\\$calcsheet, \\\$col1, \\\$col2, \\\$col3, \\\$col4) = @_;
\\\$calcSheet->ColumnHeading(\\\$col1) = "X coord";
\\\$calcSheet->ColumnHeading(\\\$col2) = "Y Coord";
\\\$calcSheet->ColumnHeading(\\\$col3) = "Z Coord";
\\\$calcSheet->ColumnHeading(\\\$col4) = "Coordination Number at \\\$coordRadius A";
}
# distanceCalc subroutine which calculates the distance between the atoms and allows for periodic boundary conditions
sub distanceCalc {
my ( \\\$hbAtom, \\\$donorAtom, \\\$length) = @_;
my \\\$d = \\\$hbAtom - \\\$donorAtom;
my \\\$x = int (\\\$d/\\\$length);
my \\\$dd = (\\\$d - (\\\$x * \\\$length));
if (\\\$dd > (\\\$length / 2)) {
\\\$dd = \\\$dd - \\\$length;
} elsif (\\\$dd < ((-1*\\\$length) / 2)) {
\\\$dd = \\\$dd + \\\$length;
}
return \\\$dd;
}
I also attached the the output study table, the last column is the coordination number. I calculated 5 frames. If I only calculated the first frame, the coordination number was 5, instead of 25.
X coord Y Coord Z Coord Coordination Number at 4.43 A
-20.50616956 28.78760147 48.69982107 25
-20.42534202 28.64093861 48.40114961 24
-20.34123772 28.35878403 48.32036261 25
-20.42949070 28.27443649 48.53806348 25
-20.66543016 28.18335445 48.68613123 24
Any suggestion is appreciated!
All the best,
Jingjing
I modified Stephen Todd's CoordinationNumberDistance script to calculate the evolution of coordination number for a trajectory. I didn't integrate the RDF curves because I'm calculating the CN for the surface species, the coordinating ions were not homogeneous distributed. I looped over the frames, created temporary structure files to calculate the coordination number for each frame, then tried to write the coordination number of each frame into a study table. However, this script didn't calculate the coordination number properly. It added the coordination number of all the frames, and put the sum into the cells where it was supposed to be the coordination number of each frame. Here is my script:
use strict;
use MaterialsScript qw(:all);
######################################################################################
# Begin User Editable Settings
my \\\$doc = \\\$Documents{"new-140000-200001.xtd"}; # Set this to the start document
my \\\$coordRadius =4.43; # The coordination distance radius defining the coordination sphere
my \\\$coordinatedAtom = "Ca1"; # The coordinated atom
my \\\$coordinatingAtom = "O*"; # The coordinating atoms
my \\\$maxAtoms = 1000; # Maximum number of atoms for the structure to be included
#
# End User Editable Settings
########################################################################################
# Sets the working document
# Sets the lattice 3D orientation
\\\$doc->Lattice3D->OrientationConvention = "A along X, B in XY plane";
# Check the number of atoms and sets variables based on them
my \\\$numAtoms = \\\$doc->SymmetrySystem->NumberOfAtoms;
print "The number of atoms in the system is \\\$numAtoms.\\n";
my \\\$addDoc = 0; # Variable used to determine whether the add structure to the study table
my \\\$arrayTerms = 3; # If you are adding the structure, the number of terms in the array changes
# from three to four (three coordinates + atom)
# Sets the values of addDoc and arrayTerms depending on the number of atoms in the system
if (\\\$numAtoms > \\\$maxAtoms) {
\\\$addDoc = 0;
\\\$arrayTerms = 3;
print "The number of atoms exceeds the maxAtoms, the document will not be added to the study table.\\n";
} else {
\\\$addDoc = 1;
\\\$arrayTerms = 4;
print "The number of atoms is less than the maxAtoms, the document will be added to the study table.\\n";
}
# Set up the study table and define the number of columns depending on whether a structure
# column is going to be added. This uses the createColumns subroutine.
my \\\$newStudyTable = Documents->New("CoordRadius\\\$coordRadius.std");
my \\\$calcSheet = \\\$newStudyTable->Sheets->Item(0);
if (\\\$addDoc == 1) {
createColumns(\\\$calcSheet, 1, 2, 3, 4);
\\\$calcSheet->ColumnHeading(0) = "Cell";
print "Created study table with column for cell.\\n";
} elsif (\\\$addDoc == 0) {
createColumns(\\\$calcSheet, 0, 1, 2, 3);
print "Created study table without column for cell.\\n";
} else {
die "There seems to be a problem with the value of addDoc";
}
# Gets the lengths of the axis of the system for dealing with periodic boundary conditions later
my \\\$lengthX = \\\$doc->Lattice3D->LengthA;
my \\\$lengthY = \\\$doc->Lattice3D->LengthB;
my \\\$lengthZ = \\\$doc->Lattice3D->LengthC;
# Gets all the atoms into a collection
my \\\$atoms = \\\$doc->UnitCell->Atoms;
# Finds the display style on the first atom
my \\\$displayStyle = \\\$doc->UnitCell->Atoms(0)->DisplayStyle;
############################################################################################################
# Go through all the atoms and assigns them to either coordinating or coordinated arrays depending on their
# element symbol. You can change this to forcefield type or whatever you wish.
my @coordinatingAtomsXYZ; # Array to hold coordinating atoms coordinates
my @coordinatedAtomsXYZ; # Array to hold coordinated atoms coordinates
my \\\$rowCounter = 0; # Simple row counter for writing to the study table - used later as well
for(my \\\$frame = 40001; \\\$frame <= 60002; \\\$frame = ++\\\$frame) {
#foreach my \\\$atom (@\\\$atoms){
\\\$doc->Trajectory->CurrentFrame = \\\$frame;
my \\\$tmpdoc = Documents->New("tmp.xsd");
\\\$tmpdoc->CopyFrom(\\\$doc);
my \\\$Aatoms = \\\$tmpdoc->UnitCell->Sets("\\\$coordinatedAtom")->Atoms;
my \\\$Batoms = \\\$tmpdoc->UnitCell->Sets("\\\$coordinatingAtom")->Atoms;
foreach my \\\$atom (@\\\$Aatoms) {
#if (\\\$atom->ElementSymbol eq \\\$coordinatedAtom) {
# Writes the coordinates into an array of the form (x1, Y1, Z1, Atom1, X2, Y2... Xn, Yn, Zn, Atomn)
push (@coordinatedAtomsXYZ, \\\$atom->X);
push (@coordinatedAtomsXYZ, \\\$atom->Y);
push (@coordinatedAtomsXYZ, \\\$atom->Z);
# Adds in the link to the atom if this is needed. If you are not adding the structure, this is ignored.
# Also adds the x, y, and z data to the study table, depending on whether the structure is needed.
if (\\\$addDoc == 1) {
push (@coordinatedAtomsXYZ, \\\$atom);
\\\$calcSheet->Cell(\\\$rowCounter, 1) = \\\$atom->XYZ->X;
\\\$calcSheet->Cell(\\\$rowCounter, 2) = \\\$atom->XYZ->Y;
\\\$calcSheet->Cell(\\\$rowCounter, 3) = \\\$atom->XYZ->Z;
} elsif (\\\$addDoc == 0) {
\\\$calcSheet->Cell(\\\$rowCounter, 0) = \\\$atom->XYZ->X;
\\\$calcSheet->Cell(\\\$rowCounter, 1) = \\\$atom->XYZ->Y;
\\\$calcSheet->Cell(\\\$rowCounter, 2) = \\\$atom->XYZ->Z;
}
++\\\$rowCounter; #Increment the row counter for the study table data
# Next part of the loop creates the coordinating atom array. This can be removed or replaced if you
# want to study something different like eg close contacts.
#} elsif (\\\$atom->ElmentSymbol eq \\\$coordinatingAtom) {
}foreach my \\\$atom (@\\\$Batoms){
# Writes the coordinates into an array of the form (x1, Y1, Z1, Atom1, X2, Y2... Xn, Yn, Zn, Atomn)
push (@coordinatingAtomsXYZ, \\\$atom->X);
push (@coordinatingAtomsXYZ, \\\$atom->Y);
push (@coordinatingAtomsXYZ, \\\$atom->Z);
# Adds in the link to the atom id if this is needed. If you are not adding the structure, this is ignored
if (\\\$addDoc == 1 ) {
push (@coordinatingAtomsXYZ, \\\$atom);
}
#}
}
# Calculate the number of atoms in each array. This is dependent on whether the atom id was added to the array
my \\\$numCoordinatedAtoms = @coordinatedAtomsXYZ / \\\$arrayTerms;
my \\\$numCoordinatingAtoms = @coordinatingAtomsXYZ / \\\$arrayTerms;
print "There are \\\$numCoordinatingAtoms coordinating atoms of element \\\$coordinatingAtom.\\n";
print "There are \\\$numCoordinatedAtoms coordinated atoms of element \\\$coordinatedAtom.\\n";
#####################################################################################################
# Go through each coordinated atom and check how many coordinating atoms are within the coordRadius.
my \\\$distanceX = 0; # Define some variables to store the distances in
my \\\$distanceY = 0;
my \\\$distanceZ = 0;
\\\$rowCounter = 0;
for (my \\\$coordinatedAtomCounter = 0; \\\$coordinatedAtomCounter < \\\$numCoordinatedAtoms; ++\\\$coordinatedAtomCounter) {
# Resets the coordination number counter to zero
my \\\$coordinationNumber = 0;
# Defines the working atom list array - this is used to store the atoms whose display style
# you will change. This is not really needed if the addAtoms variable is 1.
my @workingAtomList = ();
# Gets the atom coordinates from the array for the coordinated atoms
my \\\$coordinatedAtomX = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms)];
my \\\$coordinatedAtomY = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +1];
my \\\$coordinatedAtomZ = @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +2];
if (\\\$addDoc == 1) {
push (@workingAtomList, @coordinatedAtomsXYZ[(\\\$coordinatedAtomCounter * \\\$arrayTerms) +3]);
}
# Go through all the coordinating atoms to calculate the coordinatation number
for (my \\\$coordinatingAtomCounter = 0; \\\$coordinatingAtomCounter < \\\$numCoordinatingAtoms; ++\\\$coordinatingAtomCounter) {
# Reads the coordinates from the array which has the form (x1, Y1, Z1, X2, Y2... Xn, Yn, Zn)
my \\\$coordinatingAtomX = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms)];
my \\\$coordinatingAtomY = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms) + 1];
my \\\$coordinatingAtomZ = @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * \\\$arrayTerms) + 2];
# Use the distanceCalc subroutine to calculate the distance between the atoms
\\\$distanceX = distanceCalc (\\\$coordinatingAtomX,\\\$coordinatedAtomX, \\\$lengthX);
\\\$distanceY = distanceCalc (\\\$coordinatingAtomY,\\\$coordinatedAtomY, \\\$lengthY);
\\\$distanceZ = distanceCalc (\\\$coordinatingAtomZ,\\\$coordinatedAtomZ, \\\$lengthZ);
# Using the different calculations x y, and z distances, calculate the actual vector
my \\\$vector = sqrt (\\\$distanceX*\\\$distanceX + \\\$distanceY*\\\$distanceY + \\\$distanceZ*\\\$distanceZ);
# Check if the vector is less than the coordRadius and update the coordinationNumber if successful
if (\\\$vector <= \\\$coordRadius) {
++\\\$coordinationNumber;
if (\\\$addDoc == 1) {
push (@workingAtomList, @coordinatingAtomsXYZ[(\\\$coordinatingAtomCounter * 4) + 3]);
}
}
}
# Update the study table with the coordination number and/or document with atoms highlighted
if (\\\$addDoc == 1) {
# Use the working atom list array to change the display style of the coordinating and
# coordinated atoms.
foreach my \\\$workingAtom (@workingAtomList) {
\\\$workingAtom->DisplayStyle = "Ball and stick";
}
# Update the study table
\\\$calcSheet->Cell(\\\$rowCounter, 0) = \\\$tmpdoc;
\\\$calcSheet->Cell(\\\$rowCounter, 4) = \\\$coordinationNumber;
#Reset the display Style on the atoms that have just been changed
foreach my \\\$workingAtom (@workingAtomList){
\\\$workingAtom -> DisplayStyle = \\\$displayStyle;
}
} else {
# Just update the study table
\\\$calcSheet->Cell(\\\$rowCounter, 3) = \\\$coordinationNumber;
}
++\\\$rowCounter; # Updates the row counter
}
\\\$tmpdoc->Discard;
}
print "\\nCalculation Complete \\n";
####################################################################################################
# Subroutine section
# Subroutine to create the columns in the study table, depending on whether a structure is required
sub createColumns {
my (\\\$calcsheet, \\\$col1, \\\$col2, \\\$col3, \\\$col4) = @_;
\\\$calcSheet->ColumnHeading(\\\$col1) = "X coord";
\\\$calcSheet->ColumnHeading(\\\$col2) = "Y Coord";
\\\$calcSheet->ColumnHeading(\\\$col3) = "Z Coord";
\\\$calcSheet->ColumnHeading(\\\$col4) = "Coordination Number at \\\$coordRadius A";
}
# distanceCalc subroutine which calculates the distance between the atoms and allows for periodic boundary conditions
sub distanceCalc {
my ( \\\$hbAtom, \\\$donorAtom, \\\$length) = @_;
my \\\$d = \\\$hbAtom - \\\$donorAtom;
my \\\$x = int (\\\$d/\\\$length);
my \\\$dd = (\\\$d - (\\\$x * \\\$length));
if (\\\$dd > (\\\$length / 2)) {
\\\$dd = \\\$dd - \\\$length;
} elsif (\\\$dd < ((-1*\\\$length) / 2)) {
\\\$dd = \\\$dd + \\\$length;
}
return \\\$dd;
}
I also attached the the output study table, the last column is the coordination number. I calculated 5 frames. If I only calculated the first frame, the coordination number was 5, instead of 25.
X coord Y Coord Z Coord Coordination Number at 4.43 A
-20.50616956 28.78760147 48.69982107 25
-20.42534202 28.64093861 48.40114961 24
-20.34123772 28.35878403 48.32036261 25
-20.42949070 28.27443649 48.53806348 25
-20.66543016 28.18335445 48.68613123 24
Any suggestion is appreciated!
All the best,
Jingjing
