pDynamo3/pSimulation/GrowingStringPath.py at main · pdynamo/pDynamo3 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
"""Path generation using the growing string method."""

from pBabel                                 import ExportTrajectory
from pCore                                  import Align         , \
                                                   Clone         , \
                                                   logFile       , \
                                                   LogFileActive
from pMolecule                              import SystemGeometryObjectiveFunction
from pScientific.Arrays                     import Array
from pScientific.ObjectiveFunctionIterators import LBFGSMinimizer

#===================================================================================================================================
# . Notes:
#
# . The number of images is nimages and are numbered from 0 to nimages-1. The first (0) and last images (nimages-1) are fixed so
#   only those in the range [1,nimages-2] move.
#
# . The overall characteristics of the optimization seem to be similar to that of the original although some of the rms gradient
#   tolerance and stepSize formulae may need playing with. The module, however, should provide a foundation for future work.
#
#===================================================================================================================================

# . Errors.
class ArgumentError ( Exception ): pass

#===================================================================================================================================
# . "Growing string" method for generating an initial path.
#===================================================================================================================================
# . Basic minimization options.
_MinimizerOptions = { "logFrequency"         :   1 ,
                      "maximumIterations"    :  50 ,
                      "rmsGradientTolerance" : 1.5 }

# . Change path to trajectory object (as for all other scripts)?

def GrowingStringInitialPath ( system, npoints, point0, pointn, path, **options ):
    """Use a \"growing string\" method for generating an initial path."""

    # . Get some options.
    log      = options.pop ( "log", logFile )
    sequence = options.pop ( "sequence", "Forward" ).lower ( )

    # . Get the number of points.
    npoints = max ( npoints, 2 )

    # . Save the current system coordinates.
    temporary = Clone ( system.coordinates3 )

    # . Create the trajectory and write the first frame.
    trajectory = ExportTrajectory ( path, system )
    trajectory.WriteHeader ( )
    trajectory.WriteFrame ( point0, frame = 0 )

    # . Get the first and last energies.
    energies = {}
    point0.CopyTo ( system.coordinates3 )
    energies[0] = system.Energy ( log = log )
    pointn.CopyTo ( system.coordinates3 )
    if system.freeAtoms is None: system.coordinates3.Superimpose ( point0 )
    energies[npoints-1] = system.Energy ( log = log )

    # . Save the last point.
    trajectory.WriteFrame ( system.coordinates3, frame = npoints - 1 )

    # . Generate intermediate points.
    if npoints > 2:

        # . Set up the optimizer.
        algorithm = options.pop ( "optimizer", LBFGSMinimizer )
        options.update ( _MinimizerOptions )
        optimizer = algorithm.WithOptions ( **options )
        optimizer.Summary ( log = log )

        # . Set up the objective function.
        of = SystemGeometryObjectiveFunction.FromSystem ( system )

        # . Get variables for the first and last points and the tangent constraint.
        first   = Array.WithExtent ( of.nVariables )
        last    = Array.WithExtent ( of.nVariables )
        tangent = Array.WithExtent ( of.nVariables )

        # . Save the last point (given it is already in system.coordinates3).
        of.VariablesGet ( last )

        # . Save the first point.
        point0.CopyTo ( system.coordinates3 )
        of.VariablesGet ( first )

        # . Generate the sequence over which iterations are to occur.
        sequenceData = []
        if   sequence == "alternate" :
            for i in range ( ( npoints - 2 ) // 2 ):
                sequenceData.append ( ( first, last, i + 1, npoints - 2 - ( 2 * i ) ) )
                sequenceData.append ( ( last, first, npoints - 2 - i, npoints - 2 - ( 2 * i + 1 ) ) )
            if ( npoints - 2 ) % 2 != 0: sequenceData.append ( ( first, last, ( npoints - 2 ) // 2 + 1, 1 ) )
        elif sequence == "backward"  :
            for i in range ( npoints - 2 ): sequenceData.append ( ( last, first, npoints - 2 - i, npoints - 2 - i ) )
        elif sequence == "forward"   :
            for i in range ( npoints - 2 ): sequenceData.append ( ( first, last, i + 1, npoints - 2 - i ) )
        else: raise ArgumentError ( "Unknown sequence option - " + sequence + "." )

        # . Loop over points.
        for ( start, stop, index, remainder ) in sequenceData:

            # . Get the new point to optimize and the tangent constraint.
            stop.CopyTo ( tangent )
            tangent.Add ( start, scale = -1.0 )
            start.Add ( tangent, scale = 1.0 / float ( remainder + 1 ) )
            tangent.Normalize ( )

            # . Set up the object function.
            of.RemoveRotationTranslation ( reference = point0 )
            of.AddLinearConstraint ( tangent )
            of.VariablesPut ( start )

            # . Minimize and save the point.
            optimizer.Iterate ( of, log = log )
            of.VariablesGet ( start )
            energies[index] = system.Energy ( log = log )
            trajectory.WriteFrame ( system.coordinates3, frame = index )

    # . Finish up.
    system.coordinates3 = temporary
    trajectory.WriteFooter ( )
    trajectory.Close ( )

    # . Output.
    if LogFileActive ( log ):
        items = [ ( "{:d}".format ( key ), "{:.3f}".format ( energies[key] ) ) for key in sorted ( energies.keys ( ) ) ]
        c1    = max ( [ len ( k ) for ( k, _ ) in items ] )
        c2    = max ( [ len ( v ) for ( _, v ) in items ] )
        table = log.GetTable ( columns = [ max ( 12, c1 ), max ( 14, c2 ) ] )
        table.Start   ( )
        table.Title   ( "Growing String Energies" )
        table.Heading ( "Structure" )
        table.Heading ( "Energy"    )
        for ( k, v ) in items:
            table.Entry ( k, align = Align.Left )
            table.Entry ( v )
        table.Stop ( )

#===================================================================================================================================
# . Testing.
#===================================================================================================================================
if __name__ == "__main__" :
    pass