计算化学公社

标题: 谁有forcite计算玻璃化温度的脚本？ [打印本页]

作者
Author: zyj19831206 时间: 2016-3-8 18:38
标题: 谁有forcite计算玻璃化温度的脚本？
教程里计算玻璃化温度的过程是如图片所示，这个tg.pl脚本是怎么得到的？或者用synthia模块可以计算得到？

作者
Author: sobereva 时间: 2016-3-8 19:28
如果幻灯片上提到的那个文章及补充材料里没有这个文件，可能是开培训班的人从网上某处找来的

作者
Author: zyj19831206 时间: 2016-3-8 21:01
本帖最后由 zyj19831206 于 2016-3-8 21:04 编辑

sobereva 发表于 2016-3-8 19:28
如果幻灯片上提到的那个文章及补充材料里没有这个文件，可能是开培训班的人从网上某处找来的

google了半天都没有，我是想看看论坛里可有人有，我没有参加这个培训，

作者
Author: 卡开发发 时间: 2016-3-9 01:03
本帖最后由卡开发发于 2016-3-9 01:16 编辑

看起来应该就是优化、退火，然后不同温度下做分子动力学，最后在各自求平均密度，如果是这样大致可以使用这样的脚本：
(, 下载次数 Times of downloads: 345)
里面的选项可以按照自己的需求再进行修改。

作者
Author: zyj19831206 时间: 2016-3-9 10:24

卡开发发发表于 2016-3-9 01:03
看起来应该就是优化、退火，然后不同温度下做分子动力学，最后在各自求平均密度，如果是这样大致可以使用这 ...

请问卡兄，这个脚本哪里来的？你自己编写的？

作者
Author: 卡开发发 时间: 2016-3-9 21:22

zyj19831206 发表于 2016-3-9 10:24
请问卡兄，这个脚本哪里来的？你自己编写的？

自己写的，知道原理自己写一个简单的不难。

作者
Author: zyj19831206 时间: 2016-3-9 21:52

卡开发发发表于 2016-3-9 21:22
自己写的，知道原理自己写一个简单的不难。

原理是什么，是通过计算一系列的比表面，然后画图？

作者
Author: 卡开发发 时间: 2016-3-9 22:39

zyj19831206 发表于 2016-3-9 21:52
原理是什么，是通过计算一系列的比表面，然后画图？

比表面？不是密度的倒数？cm^3/g？看样子是不同温度下的

作者
Author: zyj19831206 时间: 2016-3-29 09:51

卡开发发发表于 2016-3-9 22:39
比表面？不是密度的倒数？cm^3/g？看样子是不同温度下的

搞到这个脚本了，问题是太长了，还不懂怎么用。
#!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;

}

作者
Author: sustal 时间: 2018-3-14 16:48
请问楼主，这个脚本怎么用，运行之前要对模型进行几何优化、退火、NVT\NPT然后再运行吗？我运行之后报错：
Try relaxing the structure with geometry optimization, reducing the time step,
reducing the temperature or a combination of these in Forcite.Dynamics (function/property "Run") at -e line 430.
求助

作者
Author: 易锦华 时间: 2021-3-18 09:15

sustal 发表于 2018-3-14 16:48
请问楼主，这个脚本怎么用，运行之前要对模型进行几何优化、退火、NVT\NPT然后再运行吗？我运行之后报错：
...

我也是，请问你最后解决了嘛

作者
Author: xiaohuidan 时间: 2022-11-11 11:38

易锦华发表于 2021-3-18 09:15
**** 作者被禁止或删除内容自动屏蔽 ****

请问您解决了吗，NVT+NPT的步骤应该怎么操作呢，万分感谢

作者
Author: xiaohuidan 时间: 2022-11-28 22:12
这个脚本咋用啊，求大神指点，跪谢

作者
Author: sxkl 时间: 2025-4-15 21:15
请问有没有这个脚本的使用教程呀

欢迎光临计算化学公社 (http://bbs.keinsci.com/)