Welcome To Support Community

Materials Studio

Advanced Search
Ask Search:
jjbjjb 

Coordination number evolution script

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
Jason DeJoannisJason DeJoannis
Hi Jingjing,
The problem is that you are calculating the coordination number for each A atom with every B atom from every frame! I suggest rewriting the script to carry out the calculation within the main loop so that you can wipe the A and B arrays with each iteration.