|
|
搞到这个脚本了，问题是太长了，还不懂怎么用。
#!perl
use strict;
use MaterialsScript qw(:all);
# Author: Stephen Todd
# Version: 5
# Materials Studio Version: 5.0
# Modules: Materials Visualizer, Forcite Plus, COMPASS
# This script runs numerous cycles of NPT dynamics at different temperatures and plots the density,
# temperature, energy and optionally, the percentage of torsions that have change, in a study table.
# At the end of the simulations, the values are averaged and a new tab is created in the study table.
# This data can be used to estimate the Tg of a structure but this should be run multiple times and
# values should be averaged to get better statistics.
# A trajectory is required as input and the script iterates over each frame in the trajectory, calculating
# a separate Tg for each structure.
# After the density is equilibrated, the production NPT dynamics is performed and the density and energy
# are averaged over the production run. An optional torsional analysis can be performed to estimate how
# many torsions have changed over the production run. This is based on the idea that below the glass
# transition temperature, the backbone torsions are frozen and do not change as often.
# Note that the length of equilibrationDuration is scaled by the temperature so that it increases as
# temperature decreases to longer equilibration times at lower temperatures.
###################################################################################################
# Begin editable settings
my $filename = "Polyethylene"; # Trajectory filename
# Tg settings
my $startTemp = 300; # Start temperature for the Tg scan
my $endTemp = 100; # End temperature for the Tg scan
my $tempStepSize = 25; # Step size for each run
# Perform equilibration - set this to 0 if you don't want to run equilibration
my $numberCycles = 1; # this is the number of temperature cycles in the initial equilibration
# Density equilibration settings
my $equilibrationPressure = 0.0001; # Pressure for the density equilibration
my $timeStep = 1; # Dynamics time step in fs
my $nonBondMethod = "Group based"; # Set the non-bond method - applies to both VdW and Coulomb
my $densityTolerance = 0.02; # Sets the tolerance for the density check
my $equilibrationDuration = 20.0; # Total simulation time (in ps) for the dynamics equilibration phase.
my $equilibrationFrameEvery = 1000; # Writes out a frame every X steps
my $blockSize = 4; # Average the density over this number of frames
my $maxRestarts = 20; # The maximum number of restarts for the NPT run
# Settings for the torsion angle analysis and NPT production run
my $calculateTorsions = 0; # Choose whether to perform torsion analysis
my $productionDuration = 50; # Length of production calculation (ps)
my $productionFrameEvery = 5000; # Writes out a frame every X steps
my $angleWindow = 30; # Angle window for defining cis/trans etc
# End editable settings
####################################################################################################
#my $doc = $Documents{"$xsdDoc.xtd"};
my $studyTable = Documents->New("$filename.std");
my $textDoc = Documents->New("$filename"."_summary.txt");
# Converts the dynamics time to number of steps
use constant FEMTO_TO_PICO => 1000; # multiply by this to convert from femto-unit to pico-unit.
# Define a hash list to hold the torsion data
my %allTorsions;
#################################################################################
# Initialises Forcite with the basic energy settings
my $forcite = Modules->Forcite;
$forcite->ChangeSettings([ Quality => "Medium",
'3DPeriodicElectrostaticSummationMethod' => $nonBondMethod,
'3DPeriodicvdWSummationMethod' => $nonBondMethod,
CurrentForcefield => "COMPASS",
TimeStep => $timeStep,
TrajectoryRestart => 'false',
AppendTrajectory => 'false',
WriteLevel => "Silent",
Ensemble3D => "NPT",
Thermostat => "Nose",
Barostat => "Andersen",
]);
###################################################################################
# Split the trajectory and set up the calculations
my @frameNames = splitTheTrajectory($filename);
# Now go over each frame in frameNames and run the calculations.....
foreach my $frameName (@frameNames) {
# Create the new tab in the study table
my $newSheet = $studyTable->InsertSheet;
$newSheet->Title = "$frameName";
$newSheet->ColumnHeading(0) = "Temperature (K)";
$newSheet->ColumnHeading(1) = "Density g/cm3";
$newSheet->ColumnHeading(2) = "Total Energy kcal/mol";
$newSheet->ColumnHeading(3) = "%change in backbone torsions over $productionDuration ps.";
# Specify the document to run the calculation on
my $doc = $Documents{"$frameName.xsd"};
# Run an initial equilibration run at the start temperature
$textDoc->Append( "\n#######################################################\n");
$textDoc->Append( "Starting temperature equilibration run on $frameName\n");
$textDoc->Append( "#######################################################\n\n");
$textDoc->Save;
my $initialTrajectory = initialEquilibration($doc, $numberCycles);
################################################################################
# Now go through the run the different temperature calculation on the individual
# frame.
my $studyTableRowCounter = 0;
for (my $temperature = $startTemp; $temperature >= $endTemp; $temperature = $temperature - $tempStepSize) {
# Initialise some variables
my $averageDensity = undef;
my $densityCondition = undef;
my $totalEnergy = undef;
my $firstCycle = undef;
my $equilibrationTrajectory = undef;
# Calculate the new equilibration number of steps, weighted by the temperature
my $adjustedEquilibrationDuration = $equilibrationDuration *$startTemp/$temperature;
my $equilibrationSteps = ($adjustedEquilibrationDuration * FEMTO_TO_PICO / $timeStep);
my $productionSteps = $productionDuration * FEMTO_TO_PICO / $timeStep;
# Create the first frame from the trajectory
$doc = Documents->New("$frameName"."_$temperature.xsd");
$initialTrajectory->Trajectory->CurrentFrame = $initialTrajectory->Trajectory->EndFrame;
$doc->CopyFrom($initialTrajectory);
$initialTrajectory->Close;
$textDoc->Append( "\n############################################################\n");
$textDoc->Append( "Starting NPT equilibration run on $frameName at $temperature K\n");
$textDoc->Append( "#############################################################\n\n");
$textDoc->Save;
# Set the extra settings for the equilibration NPT dynamics run
$forcite->ChangeSettings([ Pressure => $equilibrationPressure,
Temperature => $temperature,
InitialVelocities => "Random",
NumberOfSteps => $equilibrationSteps,
TrajectoryFrequency => $equilibrationFrameEvery,
]);
eval {$firstCycle = $forcite->Dynamics->Run($doc)};
if ($@) {
# Something bad has happened, move onto the next one
$textDoc->Append( "\n!!!!!!!!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!!!!\n");
$textDoc->Append( "Failed in first stage of dynamics. Error is\n");
$textDoc->Append( $@);
$textDoc->Save;
my $densityCondition = "failed";
} else {
$doc->Close;
$equilibrationTrajectory = $firstCycle->Trajectory;
$totalEnergy = $firstCycle->TotalEnergy;
$textDoc->Append( "Initial equilibration run of $adjustedEquilibrationDuration"."ps complete, checking density\n");
$textDoc->Save;
################################################################################
# Calculates the density for the trajectory and see if it has converged
# Defines the numframes and startframes variables. These are used to work through sections of the
# trajectory later on.
my $numFrames = $equilibrationTrajectory->Trajectory->NumFrames;
my $startFrame = 1;
# Uses the subroutine density_check to check for convergence of the density
($densityCondition, $averageDensity) = density_check($equilibrationTrajectory, $startFrame, $numFrames);
$textDoc->Append( "Density has $densityCondition convergence criteria\n");
$textDoc->Append( "Average density is $averageDensity\n");
$textDoc->Save;
##############################################################################
# While loop that keeps running NPT dynamics until the density has converged.
my $restartCycles = 1; # A counter for the number of NPT cycles
while ($densityCondition eq "failed" ) {
$textDoc->Append( "\nRestarting another $equilibrationDuration"."ps NPT cycle number $restartCycles \n");
$textDoc->Save;
my $results = undef;
eval {$results = $forcite->Dynamics->Run($equilibrationTrajectory, [ TrajectoryRestart => 'true',
AppendTrajectory => 'true',
])};
if ($@) {
# Something bad has happened, move onto the next one
$textDoc->Append( "\n!!!!!!!!!!!!!!!!!!! ERROR !!!!!!!!!!!!!!!!!!!\n");
$textDoc->Append( "Failed in dynamics cycle stage $restartCycles\n");
$textDoc->Append( $@);
$textDoc->Save;
last;
}
# Resets the values of startframe and numframe so that the density check is only performed on the
# new frames added in the last trajectory cycle. Uses density_check subroutine to check for convergence.
$totalEnergy = $results->TotalEnergy;
$startFrame = $numFrames;
$numFrames = $equilibrationTrajectory->Trajectory->NumFrames;
$textDoc->Append( "Starting frame is $startFrame and number of frames is $numFrames.\n");
$textDoc->Append( "Total energy for this section is $totalEnergy\n");
($densityCondition, $averageDensity) = density_check($equilibrationTrajectory, $startFrame, $numFrames);
$textDoc->Append( "Density has $densityCondition convergence criteria\n");
$textDoc->Append( "Average density is $averageDensity\n");
$textDoc->Save;
# Increments restartCycles to keep track of how many restarts are occuring
++$restartCycles;
# Extra check so it doesn't run forever if convergence is tricky....
if (($restartCycles == $maxRestarts) && ($densityCondition eq "failed")) {
$textDoc->Append( "Warning: Exceeded maximum number of NPT cycles without density convergence\n");
$textDoc->Save;
# The last command breaks out of the current loop, ie the while loop
# and returns to the foreach loop so that it starts on the next frame
last;
}
}
}
################################################################################
# If the density has converged, put the values in the study table
$textDoc->Append( "\nDensity convergence check has $densityCondition for trajectory.\n");
$textDoc->Save;
####################################################################################
# Now do the Torsions check
my $percentageChange = -1;
#Run a production NPT dynamics calculation
my $startFrame = $equilibrationTrajectory->Trajectory->NumFrames;
$forcite->Dynamics->Run($equilibrationTrajectory, [ TrajectoryRestart => 'true',
AppendTrajectory => 'true',
NumberOfSteps => $productionSteps,
TrajectoryFrequency => $productionFrameEvery,
]);
# find last frame
my $endFrame = $equilibrationTrajectory->Trajectory->NumFrames;
# If the document has torsions, evaluate the torsion changes
my $hasTorsions = $equilibrationTrajectory->UnitCell->Torsions->Count;
if ($hasTorsions == 0) {
$textDoc->Append( "There are no torsions defined in the document.\nNot performing NVT run\n");
$percentageChange = "No torsions";
$textDoc->Save;
} else {
$percentageChange = evaluateTorsionChanges($equilibrationTrajectory, $startFrame, $endFrame);
}
my ($averageProductionDensity, $averageProductionEnergy) = averageProperties($equilibrationTrajectory, $startFrame, $endFrame);
# Update the study table with the values. If it fails to converge, write a warning
$newSheet->Cell($studyTableRowCounter, 0) = $temperature;
if ($densityCondition eq "failed") {
$newSheet->Cell($studyTableRowCounter, 1) = $averageProductionDensity;
$newSheet->Cell($studyTableRowCounter, 2) = $averageProductionEnergy;
$newSheet->Cell($studyTableRowCounter, 3) = $percentageChange;
$newSheet->Cell($studyTableRowCounter, 4) = "Warning: Failed to equilibrate";
} else {
$newSheet->Cell($studyTableRowCounter, 1) = $averageProductionDensity;
$newSheet->Cell($studyTableRowCounter, 2) = $averageProductionEnergy;
$newSheet->Cell($studyTableRowCounter, 3) = $percentageChange;
}
++$studyTableRowCounter;
# Reset initial trajectory to copy the last frame over.
$initialTrajectory = $equilibrationTrajectory;
}
}
# Calculates averages for all columns in the study table apart from the one passed in the
# subroutine - in this case, column 0 which is the energies.
calculateAveragesForStudyTable($studyTable, 0);
$textDoc->Append( "\n Calculation complete.\n");
##################################################################################
# Splits the trajectory into individual XSD's in appropriate folders, optionally
# calculate the backbone torsions
sub splitTheTrajectory {
my $trjName = shift;
my $trj = $Documents{"$trjName.xtd"};
# Calculate the torsions in the backbone atom (optionally)
if ($calculateTorsions == 1) {
mainFindBackboneAtoms($trj);
}
$textDoc->Save;
for (my $frameCounter = 1; $frameCounter <= $trj->Trajectory->NumFrames; ++$frameCounter) {
#Sets the current frame and creates a new document
$trj->Trajectory->CurrentFrame = $frameCounter;
my $frameDoc = Documents->New("$filename"."$frameCounter.xsd");
# Copies from current frame into new document
$frameDoc->CopyFrom($trj);
$frameDoc->Close;
# Stores the frame name
push (@frameNames, "$filename"."$frameCounter");
}
return @frameNames;
}
#######################################################################################
# The initial equilibration run, a short dynamics run with a small timestep followed
# by a temperature cycle.
sub initialEquilibration {
my $doc = shift;
my $cycles = shift;
# run a quick minimization
my $optResults;
eval {$optResults = $forcite->GeometryOptimization->Run($doc);};
if ($@) {
$textDoc->Append( "There is a problem with the geometry optimization\n");
$textDoc->Append( $@);
$textDoc->Save;
} else {
# Now run a dynamics with a short timestep, NVE ensemble
my $shortDynamics = $forcite->Dynamics->Run($doc, ([ TimeStep => 0.5,
NumberOfSteps => 1000,
Ensemble3D => "NVE",
Quality => "Coarse"
]));
$textDoc->Append( "Short timestep dynamics is complete\n");
$textDoc->Save;
# now take the results and put them through a simulated annealing run
for (my $cycleCounter = 0; $cycleCounter < $cycles; ++$cycleCounter) {
# take the temperature up
for (my $temperature = $endTemp; $temperature <= $startTemp; $temperature += $tempStepSize){
$textDoc->Append( "Running at $temperature\n");
$textDoc->Save;
# The actual dynamics run
my $equilibrationTempCycle = $forcite->Dynamics->Run($shortDynamics->Trajectory, ([ Ensemble3D => "NPT",
Pressure => 0.001,
Temperature => $temperature,
NumberOfSteps => 5000,
TrajectoryFrequency => 2500,
TrajectoryRestart => 1,
AppendTrajectory => 1,
Quality => "Coarse"
]));
}
# take the temperature back down
for (my $temperature = $startTemp; $temperature >= $endTemp; $temperature -= $tempStepSize){
$textDoc->Append( "Running at $temperature\n");
$textDoc->Save;
# The actual dynamics run
my $equilibrationTempCycle = $forcite->Dynamics->Run($shortDynamics->Trajectory, ([ Ensemble3D => "NPT",
Pressure => 0.001,
Temperature => $temperature,
NumberOfSteps => 5000,
TrajectoryFrequency => 2500,
TrajectoryRestart => 1,
AppendTrajectory => 1,
Quality => "Coarse"
]));
}
}
return ($shortDynamics->Trajectory);
}
}
##############################################################################
# Density check subroutine
sub density_check {
my ($newTrj, $startFrame, $numFrames) = @_;
# Define the array to hold the density values
my @density;
# Pulls the density values from each frame and sticks them in @density
for (my $frame = $startFrame ; $frame <= $numFrames; ++$frame){
$newTrj->Trajectory->CurrentFrame = int($frame);
my $density = $newTrj->SymmetrySystem->Density;
push (@density, $density);
}
my $numValuesInArray = @density;
my @averagedBlocks;
# Go through each value in the array and add them up for each block
# If the number of frames divided by the block size is not an integer,
# the remainder is averaged, unless it is a single value.
my $blockCounter = 0;
my $densitySum = 0;
foreach (my $valueCounter = 0; $valueCounter < $numValuesInArray; ++$valueCounter) {
if (($blockCounter < $blockSize) && ($valueCounter != ($numValuesInArray -1))) {
$densitySum = $densitySum + @density[$valueCounter];
++$blockCounter;
# print " $valueCounter $blockCounter @density[$valueCounter] \n";
}
elsif ($blockCounter == $blockSize){
# reset the block counter and push the value into the array
my $blockDensity = $densitySum/$blockSize;
push (@averagedBlocks, $blockDensity);
#print "blockdensity is $blockDensity\n";
$blockCounter = 0;
$densitySum = 0;
$densitySum = $densitySum + @density[$valueCounter];
++$blockCounter;
# print "$valueCounter $blockCounter @density[$valueCounter] \n";
} else {
# This gets the last density and averages it.
$densitySum = $densitySum + @density[$valueCounter];
my $blockDensity = $densitySum/($blockCounter+1);
push (@averagedBlocks, $blockDensity);
# print "$valueCounter $blockCounter @density[$valueCounter] \n";
#print "blockdensity is is $blockDensity\n";
}
}
# Now work out the difference between the blocks to check for convergence
my @differenceBlocks;
my $numberBlocks = @averagedBlocks;
for (my $blockCounter = 1; $blockCounter < $numberBlocks; ++$blockCounter) {
my $difference = sqrt((@averagedBlocks[$blockCounter] - @averagedBlocks[$blockCounter-1])**2);
push (@differenceBlocks, $difference);
}
# Finally, go through each difference and see if they are smaller than the
my $densityCondition = "passed";
foreach my $value (@differenceBlocks) {
if ($value > $densityTolerance) {
$densityCondition = "failed";
}
}
# And calculate the average density over the cycle
my $totalDensity = 0;
foreach my $value (@averagedBlocks) {
$totalDensity = $totalDensity + $value;
}
my $averageDensity = $totalDensity/$numberBlocks;
# Returns the density_check variables
return $densityCondition, $averageDensity;
}
###############################################################################
# Find the main backbone atoms in the system
sub mainFindBackboneAtoms {
my $doc = shift;
my $atoms = $doc->UnitCell->Atoms;
# Define a global array to store the branch atoms - used later
my @branches;
# Go through all the atoms and store their names in an arrray called elements
# If a backbone atom is found, change it's name to backbone.
my @elements;
my $molecules = $doc->UnitCell->Molecules;
my $moleculeCounter = 0;
foreach my $molecule (@$molecules) {
# rename the molecules
$molecule->Name = "$moleculeCounter";
++$moleculeCounter;
# cycle through the atoms renaming the backbone atoms
my $atoms = $molecule->Descendants->Atoms;
foreach my $atom (@$atoms) {
push (@elements, $atom->Name);
if ($atom->IsBackboneAtom) {
$atom->Name = "backbone";
}
}
}
# now find the molecules that have backbone atoms and run on those only.
# if no molecules exist like that.... don't run anything.
my @backboneMolecules;
foreach my $molecule (@$molecules) {
my $atoms = $molecule->Descendants->Atoms;
foreach my $atom (@$atoms) {
push (@elements, $atom->Name);
if ($atom->IsBackboneAtom) {
push (@backboneMolecules, $molecule->Name);
last;
}
}
}
# check to see how many molecules have backbone atoms
my $anyBackboneMolecules = 0;
my $numberBackboneMolecules = @backboneMolecules;
if ($numberBackboneMolecules < 1) {
$anyBackboneMolecules = 0;
$textDoc->Append( "Can't find any any backbone molecules\n");
} else {
$anyBackboneMolecules = 1;
$textDoc->Append( "There are $numberBackboneMolecules molecules containing backbone atoms \n");
}
$textDoc->Save;
foreach my $backboneMolecule (@backboneMolecules) {
# Create an array to hold the torsions for each molecule and get the name of the molecule
my $torsionCounter = 0;
my @moleculeTorsions;
my $molecule = $doc->UnitCell->Molecules("$backboneMolecule");
$textDoc->Append( "Processing molecule $backboneMolecule\n");
my $continue = 1; # $continue defines whether all the atoms have been examined
my $counter = 0; # counter for the naming of the atoms
# Put all the backbone atoms in an array
my @backboneAtoms = getAllBackboneAtoms($molecule);
my $numberBackboneAtoms = @backboneAtoms;
$textDoc->Append( "There are $numberBackboneAtoms backbone atoms in $backboneMolecule molecule\n");
# Find the terminal atom and assign it to atom1
my $atom1 = findTerminalAtom(@backboneAtoms);
$atom1->Name = "$counter";
# Create some variables to hold the other 3 atoms in the torsion
my $atom2;
my $atom3;
my $atom4;
my $endOfChain; # Variable to tell whether you have reached the end of the chain
# For the first set of atoms, define the atoms based on the naming of the previous
# atoms. The findNextBackboneAtom looks for an attached atom called backbone. If
# finds one, it is defined as the next atom in the torsion.
++$counter;
($atom2, $endOfChain) = findNextBackboneAtom($atom1,$counter,@branches);
++$counter;
($atom3, $endOfChain) = findNextBackboneAtom($atom2,$counter, @branches);
++$counter;
($atom4, $endOfChain) = findNextBackboneAtom($atom3,$counter,@branches);
# Create a torsion from the atoms and push to store it in the molecule array and the all torsions array
my $torsion = $doc->CreateTorsion([$atom1, $atom2, $atom3, $atom4]);
$torsion->Name = "$backboneMolecule"."_$torsionCounter";
++$torsionCounter;
# push the torsion name into the %allTorsions to store
$allTorsions{$torsion->Name};
# Now move the torsions along the backbone by one atom and create a new torsion
# while there are backbone atoms present.
while ($continue == 1) {
++$counter;
$atom1 = $atom2;
$atom2 = $atom3;
$atom3 = $atom4;
($atom4, $endOfChain) = findNextBackboneAtom($atom3,$counter,@branches);
# If the last atom found was valid, create a new torsion
if ($endOfChain == 0) {
$torsion = $doc->CreateTorsion([$atom1, $atom2, $atom3, $atom4]);
$torsion->Name = "$backboneMolecule"."_$torsionCounter";
++$torsionCounter;
# push the torsion name into the allTorsions store
$allTorsions{$torsion->Name};
# If the last atom wasn't valid, you have found the end of the chain
# Now you have to start looking for branches.
} elsif ($endOfChain == 1) {
# Check that the @braches array contains some information.
my $numberBranches = @branches;
# print "Run out of atoms, looking for branch points.\n";
if ($numberBranches > 0) {
# If a branch point if found, grab the atom name and set this to atom2
# Set the name of atom1name to the atom attached to atom2 that is not
# a named backbone atom;
$textDoc->Append( "Found a branch point, now working along the branch point\n");
my $atom2Name = @branches[0];
my $atom1Name = $atom2Name - 1;
# Remove the atom from the array;
shift(@branches);
# Go through all the backbone atoms to find the one with the branch name
foreach my $branchBackboneAtom (@backboneAtoms) {
if ($branchBackboneAtom->Name eq "$atom2Name") {
$atom2 = $branchBackboneAtom;
} elsif ($branchBackboneAtom->Name eq "$atom1Name") {
$atom1 = $branchBackboneAtom;
}
}
# Create the first torsion in the branch
($atom3, $endOfChain) = findNextBackboneAtom($atom2,$counter,@branches);
($atom4, $endOfChain) = findNextBackboneAtom($atom3,$counter, @branches);
$torsion = $doc->CreateTorsion([$atom1, $atom2, $atom3, $atom4]);
$torsion->Name = "$backboneMolecule"."_$torsionCounter";
++$torsionCounter;
# push the torsion name into the allTorsions store
$allTorsions{$torsion->Name};
} else {
# If there are no branches in the array, breakout of the while
# loop and stopping looking for atoms
$textDoc->Append( "There are no more branches.\n\n");
$continue = 0;
}
}
}
}
# And then go through the atoms renaming them again
my $atomCounter = 0;
foreach my $atom (@$atoms) {
$atom->Name = "@elements[$atomCounter]";
++$atomCounter;
}
return $anyBackboneMolecules;
}
#################################################################################
# Subroutine to get all the backbone atoms based on the IsBackboneAtom property
sub getAllBackboneAtoms {
my ($molecule) = @_;
# Get all the atoms for the particular molecule
my $atoms = $molecule->Descendants->Atoms;
my @backboneAtoms;
# If the IsBackboneAtom is true, store it in the array
foreach my $atom (@$atoms) {
if ($atom->IsBackboneAtom) {
push (@backboneAtoms, $atom);
}
}
return @backboneAtoms;
}
###############################################################################
# A subroutine to find a starting atom, or terminal backbone atom.
sub findTerminalAtom {
my (@backboneAtoms) = @_;
my $terminalAtom;
# Simple loop to look for the attached atoms and see how many have the
# IsBackboneAtom property
foreach my $atom (@backboneAtoms) {
my $attachedAtoms = $atom->AttachedAtoms;
my $backboneCounter = 0;
foreach my $attAtoms (@$attachedAtoms) {
if ($attAtoms->IsBackboneAtom) {
++$backboneCounter;
}
}
# If the backbone atom counter is 1, then stop looking
# as you have found one terminal backbone atom
if ($backboneCounter == 1) {
$terminalAtom = $atom;
# print "Found a terminal backbone atom\n";
last;
}
}
return $terminalAtom;
}
#################################################################################
# Subroutine to find the next backbone atom attached to a given atom
sub findNextBackboneAtom {
my ($currentBackboneAtom, $counter, @branches) = @_;
# Get the attached atoms to the current backbone atom
my $attachedAtoms = $currentBackboneAtom->AttachedAtoms;
my $newAtom;
my $endOfChain;
my $backboneAtomCounterCheck = 0;
# Look to see if there are any backbone atoms and count the number
# in the $backboneAtomCounterCheck variable
foreach my $tempAtom (@$attachedAtoms) {
if ($tempAtom->Name eq "backbone") {
$newAtom = $tempAtom;
++$backboneAtomCounterCheck;
}
}
# Now decide what to do. If there is more than one backbone atom, this is
# a branch point and you need to add it to the branch point array for checking
# later on.
if ($backboneAtomCounterCheck > 1) {
push(@branches, $currentBackboneAtom->Name);
$newAtom->Name = "$counter";
# printf ("Found a branch atom at %0.0f\n",$currentBackboneAtom->Name);
} else {
# If there is one or less backbone atoms, use the eval statement to see
# whether you can change the name of the atom. If you can't, you haven't found
# an atom and are at the end of the chain. In this case, you have to return a
# new value for endOfChain so it can start looking for branches.
# If there is an atom, just increment it's name property with a new value of counter.
eval{$newAtom->Name = "$counter"};
if ($@) {
$endOfChain = 1; # Condition for failing to find a new atom
} else {
$endOfChain = 0; # Condition for finding a new atom
}
}
# Return the new atom and the value for endOfChain
return $newAtom, $endOfChain;
}
#########################################################################################
# Evaluate how many torsions have changed as a percentage of the total number of torsions
sub evaluateTorsionChanges {
my ($equilibrationTrajectory,$startFrame, $lastFrame) = @_;
my $percentageChange = -1;
my $numFramesToAnalyse = $lastFrame - $startFrame;
$textDoc->Append( "\nAnalysing $numFramesToAnalyse frames\n");
$textDoc->Save;
# go through the frames and grab the torsion name to put in a hash of arrays
my $torsions = $equilibrationTrajectory->UnitCell->Torsions;
my $numberTorsions = $torsions->Count;
for (my $frameCounter = $startFrame; $frameCounter <= $lastFrame; ++$frameCounter) {
# Set the trajectory frame to the first or last frame
$equilibrationTrajectory->Trajectory->CurrentFrame = int($frameCounter);
$textDoc->Append( "Searching for torsions in frame $frameCounter\n");
# Go through all the torsions and put the cis/trans/gauche plus into the hash of arrays
foreach my $torsion (@$torsions) {
#my $angleDegrees = $torsion->Angle;
my $angle = findCisTrans($torsion->Angle);
#print "$angle : $angleDegrees\n";
push (@{$allTorsions{$torsion->Name}},$angle);
}
}
$textDoc->Save;
# Get the arrays from the hash list and compare the values in the arrays
my $changedCounter = 0;
foreach my $value (keys %allTorsions) {
my @values = @{ $allTorsions{$value}};
my $hasChanged = 0;
for (my $torsionFrameCounter = 1; $torsionFrameCounter < $numFramesToAnalyse; ++$torsionFrameCounter) {
# Compares the two values of the torsion angles. Note this could
# be extended to look at numerous values.
if (@values[$torsionFrameCounter] ne @values[$torsionFrameCounter-1]) {
# print "$value torsion has changed from frame $torsionFrameCounter from @values[$torsionFrameCounter-1] to @values[$torsionFrameCounter]\n";
$hasChanged = 1;
}
}
if ($hasChanged == 1) { ++$changedCounter;}
}
$percentageChange = ($changedCounter/$numberTorsions)*100;
$textDoc->Append( "$percentageChange torsions have changed\n");
$textDoc->Save;
# Empty the hash list for the next round of calculations
%allTorsions = ();
return $percentageChange;
}
##############################################################################
# A subroutine to read in the value of the angle and convert
# it into cis/trans/gauche plus or gauche minus
sub findCisTrans {
my $angle = shift;
my $newAngle;
use constant CIS => 0;
use constant TRANS => 180;
use constant GM => -60;
use constant GP => 60;
if (($angle > CIS-$angleWindow) && ($angle < CIS+$angleWindow)) {
$newAngle = "C";
} elsif (($angle > GM-$angleWindow) && ($angle < GM+$angleWindow)) {
$newAngle = "GM";
} elsif (($angle > GP-$angleWindow) && ($angle < GP+$angleWindow)) {
$newAngle = "GP";
} elsif (($angle < ((TRANS*-1) + $angleWindow)) || ( $angle > TRANS-$angleWindow)) {
$newAngle = "T";
} else {
# if the angle is none of the above, specify it as Other
$newAngle = "O";
}
return $newAngle;
}
#################################################################
# Put the properties to average into arrays and calculate the stats
sub averageProperties {
my ($trj, $start, $end) = @_;
# Use analysis to view the data in the study table
$forcite->ChangeSettings([ ActiveDocumentFrameRange => "$start - $end"]);
my $results = $forcite->Analysis->ViewInStudyTable($trj);
my $resultsST = $results->StudyTable;
my @energies = ();
my @densities = ();
for (my $row = 0; $row < $resultsST->RowCount; ++$row) {
push (@energies, $resultsST->Cell($row, "Hamiltonian"));
push (@densities, $resultsST->Cell($row, "Density"));
}
$resultsST->Delete;
# Use the calculateStatistics to calculate stats for the properties
my ($averageEnergy, $sdEnergy) = calculateStatistics(@energies);
my ($averageDensity, $sdDensity) = calculateStatistics(@densities);
$textDoc->Append(sprintf ("Average density is %0.3F, standard deviation is %0.3F\n",$averageDensity, $sdDensity));
$textDoc->Append(sprintf ("Average energy is %0.3F, standard deviation is %0.3F\n",$averageEnergy, $sdEnergy));
return ($averageDensity, $averageEnergy);
}
##########################################################
# Calculate the mean and standard deviation
sub calculateStatistics {
my @property = @_;
my $numberValues = @property;
my $meanSum = 0;
# calculate Mean
foreach my $value (@property) {
$meanSum += $value;
}
my $mean = $meanSum/$numberValues;
# Calculate standard deviation
my $stdDevSum = 0;
foreach my $value (@property) {
my $diff2 = ($value - $mean)**2;
$stdDevSum += $diff2;
}
my $stdDev = ($stdDevSum/($numberValues - 1))**0.5;
return $mean, $stdDev;
}
############################################################################
# Calculates the averages for a number of sheets in a study table
sub calculateAveragesForStudyTable {
my ($table, @ignoreColumns) = @_;
# Find the number of rows and columns in the one of the new tabs
my $rowCount = $table->Sheets("@frameNames[0]")->RowCount;
my $columnCount = $table->Sheets("@frameNames[0]")->ColumnCount;
# Insert a sheet and call it averages
my $averageSheet = $table->InsertSheet;
$averageSheet->Title = "Averages";
# Go through the columns first and, if the index matches one in @ignoreColumns, just make a copy.
# Otherwise, calculate the average value from the sheets.
# The averageColumnCounter is used to track the number of columns in the
# average tab. As errors are calculated, these will add an extra column.
my $averageColumnCounter =0;
for (my $columnCounter = 0; $columnCounter < $columnCount; ++$columnCounter) {
# Check to see whether the averages should be calculated for this column, set to 0 if they are required
my $isIgnored = 0;
foreach my $ignoredColumn (@ignoreColumns) {
if ($columnCounter == $ignoredColumn) {
$textDoc->Append( "Just copying the values for column ".$averageSheet->ColumnHeading($columnCounter)."\n");
$isIgnored = 1;
}
}
# Set the column name, create a new column heading for error it needs to be created
$averageSheet->ColumnHeading($averageColumnCounter) = $table->Sheets("@frameNames[0]")->ColumnHeading($columnCounter);
if ($isIgnored == 0) {$averageSheet->ColumnHeading($averageColumnCounter+1) = "Absolute Error"; }
# Now go through the rows for each tab and average them
for (my $rowCounter = 0 ; $rowCounter < $rowCount; ++$rowCounter) {
# If the column is ignored, make a copy of the values from another tab, else do the averaging thing
if ($isIgnored == 1) {
$averageSheet->Cell($rowCounter, $averageColumnCounter) = sprintf ("%0.1f", $table->Sheets("@frameNames[0]")->Cell($rowCounter,$columnCounter));
} else {
# get the values and stick them in an array
my @valuesToAverage = ();
foreach my $sheetName (@frameNames) {
push (@valuesToAverage , $table->Sheets("$sheetName")->Cell($rowCounter,$columnCounter));
}
my $numValues = @valuesToAverage;
my $sum = 0;
foreach my $value (@valuesToAverage) {
$sum+=$value;
}
my $mean = $sum/$numValues;
# Calculate the differences for each value in the array from the mean
# If the difference between the mean is negative, multiply it by -1
my @differences = ();
foreach my $value (@valuesToAverage) {
my $difference = $value - $mean;
if ($difference <= 0) { $difference = $difference *-1;}
push (@differences, $difference);
}
my @sortedDifferences = sort {$a <=> $b} @differences;
# Put the average value in the cell
$averageSheet->Cell($rowCounter, $averageColumnCounter) = sprintf ("%0.2f", $mean);
# Calculate the error and put it in the study table
$averageSheet->Cell($rowCounter, $averageColumnCounter+ 1) = sprintf ("%0.4f", @sortedDifferences[-1]);
}
}
if ($isIgnored == 0) { ++$averageColumnCounter};
++$averageColumnCounter;
}
$textDoc->Save;
} |
|
发表于 Post on 2016-3-8 18:38:28
收藏 Add to favorites
发表于 Post on 2016-3-8 19:28:13
楼主 Author
