sbco.tcl

::pwtk::sbco::run

::pwtk::sbco::run

SYNOPSIS

proc ::pwtk::sbco::run {head {init_ion_dynamics "bfgs"}} {

PURPOSE

Run the SBCO calculation, i.e., a series of pw.x contrained optimizations with a specified bond step-wise elongated, optionally followed by a neb.x calculation to verify the transition-state and activation energy.

BEWARE

This method is aplicable only for reactions whose reaction coordinate is predominantly along a given bond. A serious caveat of this method is that it may completely miss the transition state, yielding a considerably underestimated activation energy. For this reason, a post-NEB calculations can be performed.

ARGUMENTS

head -- rootname for the I/O files
init_ion_dynamics -- which ion_dynamics to use for optimizing the 1st SBCO image

RETURN VALUE

if post_neb == false, the "E(max) - E(0)" difference in eV (aka the estimation of the activation energy)
if post_neb == true, the NEB corrected activation energy, E(neb.TS) - E(0)

SOURCE

    variable sbco 

    set head sbco.[::pwtk::trim_prefix sbco $head]

    # initialize & parse SBCO ...
    
    parseSBCO_

    # calculation must be 'relax'
    
    set calculation [::pwtk::input::namelistGetVarValue CONTROL calculation trim]
    if { $calculation ne "relax" } {
        ::pwtk::error "for SBCO, the calculation must be 'relax', but the current value is '$calculation'" 1
    }

    # printout
    
    printTitle_ "Scan Bond Constrained Optimization"    
    print_ "bond (atom1 atom2)         ==  $sbco(atom1) $sbco(atom2)
mass (atom1 atom2)         ==  $sbco(mass1) $sbco(mass2)
number of steps            ==  $sbco(nsteps)
bond increment             ==  $sbco(increment) $sbco(unit)
bond tolerance             ==  $sbco(tolerance) bohr
full opt. of the 1st image ==  $sbco(first_opt)
auto mode                  ==  $sbco(auto)
post-relax calculation     ==  $sbco(post_relax)
post-NEB calculation       ==  $sbco(post_neb)"
    if { $sbco(post_neb) } {
        print_ "post-NEB CI_scheme         ==  $sbco(CI_scheme)"
    }
    ::pwtk::print ""
    
    set nsteps        $sbco(nsteps)
    set ion_dynamics  [::pwtk::input::namelistGetVarValue IONS ion_dynamics trim]
    set pwiAtmPosUnit [::pwtk::pwi::getAtmPosUnit enforce]
    switch $pwiAtmPosUnit {
        alat - bohr - angstrom - crystal {}
        default {
            ::pwtk::error "the '$pwiAtmPosUnit' units not imlpemented" 1
        }
    }
    
    # prepare for *.dat and *.axsf files
    
    set dat_fid [open "$head.dat" w]
    chan configure $dat_fid -buffering line
    puts $dat_fid "#  increment ($sbco(unit))   total energy (Ry)   total energy (eV)   dE (eV)   image-No."
    set xsf [::pwtk::pwi::getXSFPrimVec]\n\n

    #
    # run the SBCO calculations
    #
    set de_ev_max 0.0
    set iMax 0
    
    for {set istep 0} {$istep <= $nsteps} {incr istep} {

        # set I/O filenames
        set input  $head.$istep.in
        set output $head.$istep.out

        ::pwtk::print "step \#.$istep"
        
        if { $istep == 0 && $sbco(first_opt) } {

            # the initial step is special
            
            ::pwtk::print "   initial step will be fully optimized"
            if { $init_ion_dynamics ne "" } {
                ::pwtk::print "   ion_dynamics   = $init_ion_dynamics"          
                IONS "ion_dynamics = [::pwtk::squote $init_ion_dynamics]"
            } else {
                ::pwtk::print "   ion_dynamics   = $ion_dynamics"               
            }            
        } else {
            # the rest of the steps are the same
            
            if { $ion_dynamics eq "bfgs" || $ion_dynamics eq {} } {
                set ion_dynamics damp
                ::pwtk::print "Warning: ion_dynamics == 'bfgs' is incompatible with CONSTRAINTS, using ion_dynamics == '$ion_dynamics' instead"
            }
            ::pwtk::print "   ion_dynamics   = ${ion_dynamics}"
            IONS "ion_dynamics = [::pwtk::squote ${ion_dynamics}]"

            if { $istep > 0 } {
                # load the coordinates for the new SBCO run
                ::pwtk::pwi::replaceCoor [getNewCoor_ $pwiAtmPosUnit $old_output] $pwiAtmPosUnit
            }

            # make a CONSTRAINTS card
            
            CONSTRAINTS "1 $sbco(tolerance)\n'distance' $sbco(atom1) $sbco(atom2)"
        }
        ::pwtk::print "   bond increment = [format %.4f [expr {$istep * $sbco(increment)}]] $sbco(unit)\n"

        # run the pw.x calculation

        ::pwtk::outdir_postfix $head.$istep
        ::pwtk::runPW $input $output
        
        # write total energy to stdout and update data files

        writeData_
        
        if { [info exist old_totene] && $sbco(auto) && $totene < $old_totene } {
            ::pwtk::print "Auto SBCO condition is met (energy new < energy old).
   Stopping the SBCO calculation.
      energy new = $totene Ry
      energy old = $old_totene Ry
"
            set nsteps $istep
            incr istep
            break
        }
        
        # save the name of input/output file

        set old_output $output
        set old_totene $totene
    }

    # if post_relax & iMax < nsteps, run the full relaxation
    
    if { $sbco(post_relax) } {
        if { $iMax < $nsteps } {            
            ::pwtk::input::pushpop {
                print_ "Performing post-relax SBCO calculation\n"
                if { $init_ion_dynamics ne "" } {
                    IONS "ion_dynamics = [::pwtk::squote $init_ion_dynamics]"
                }
                ::pwtk::input::clear CONSTRAINTS
                ::pwtk::pwi::ATOMIC_POSITIONS_fromPWO $output
                
                incr nsteps
                set name   relax.[::pwtk::trim_prefix sbco $head].$istep
                set output $name.out
                
                ::pwtk::outdir_postfix $name
                ::pwtk::runPW $name

                # calculate the bond length
                set dist_ [postRelaxDist_ $pwiAtmPosUnit $output]
                writeData_
            }
        } else {
            ::pwtk::infoMsg "Post-relax SBCO calculation skipped because the last image\nhas the highest energy."
        }
    }

    incr nsteps; # account for 0th image
    
    close $dat_fid
    writeFile $head.axsf "ANIMSTEPS $nsteps\n\nCRYSTAL\n\n$xsf"

    # write some final info

    ::pwtk::print "SBCO calculation completed.

SBCO energies written to:    $head.dat
SBCO structures written to:  $head.axsf

Summary:
   Number of the SBCO images: $nsteps
   SBCO data (bond increments & energies):\n"
    ::pwtk::print [::pwtk::readFile $head.dat]

    # if requested, run post-NEB
    
    if { $sbco(post_neb) } {
        return [postNeb $head]
    } else {
        ::pwtk::print "   SBCO image with the highest energy:  $iMax
   * dE(max) = E($iMax) - E(0) = [format %.4f $de_ev_max] eV\n"
        return $de_ev_max
    }
}

TABLE OF CONTENTS

::pwtk::sbco::run