TABLE OF CONTENTS

::pwtk::sbco::postNeb

SYNOPSIS

proc ::pwtk::sbco::postNeb {head {n_images -3}} {

USAGE

   postNeb  head  ?n_images?

REMARK

This command is called by ::pwtk::sbco::run if post_neb == true.

PURPOSE

Perform the NEB calculation on top of the SBCO calculated data to locate the transition state (TS).

The NEB calculation is performed if and only if the SBCO was able to bracket the TS. There should be an intermediate image with the maximum energy in file $head.dat. This "maximum" image is taken as a TS guess, the image before the "maximum" SBCO image is uses as IS, and the image after the SBCO "maximum" image is used as FS. The NEB calculation is performed with $n_images between the so-located first and last image.

ARGUMENTS

BEWARE

The PATH's num_of_images is set only if the specified n_images >= 3.

If n_images <= -3 and PATH's num_of_images is not defined, then PATH's num_of_images is set to abs(n_images).

RETURN VALUE

NEB corrected activation energy (measured with respect to the 1st SBCO image).

SOURCE

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

    if { ! [file exists $head.dat] } {
        ::pwtk::error "cannot perform SBCO post-NEB calculation, data file $head.dat does not exist" 1
    }
    if { ! [string is integer -strict $n_images] } {
        ::pwtk::error "expected an integer number for n_images, but got $n_images" 1
    }
    if { $n_images >= 3 } {        
        PATH " num_of_images = $n_images"        
    } elseif { [::pwtk::input::namelistGetVarValue PATH num_of_images] eq {} && $n_images <= -3 } {
        PATH " num_of_images = [expr abs($n_images)] "
    }

    printTitle_ "Post-NEB SBCO calculation"
        
    # find the image with the maximum energy

    set ind    0
    set uphill 0
    set maxInd 0
    set old_ene_ev 1e10

    foreach line [split [::pwtk::readFile $head.dat] \n] {
        
        if { ! [regexp -- {(^[ \t]*#|^[ \t]*$)} $line] } {
            # not a comment and not an empty line
            
            set ene_ev($ind) [lindex $line 2]
                    
            if { $ind > 0 } {
                if { $uphill && $ene_ev($ind) < $old_ene_ev } {
                    # maximum located
                    set maxInd [expr {$ind - 1}]
                    break
                } elseif { $ene_ev($ind) > $old_ene_ev } {
                    set uphill 1
                }
            }       
            set old_ene_ev $ene_ev($ind)
            incr ind
        }
    }
    if { $maxInd == 0 } {
        ::pwtk::error "cannot perform post-NEB SBCO calculation, maximum not found in '$head.dat'" 1
    }
    set im [expr {$maxInd - 1}]
    set iM [expr {$maxInd + 1}]

    ::pwtk::fileMustExist $head.$im.out
    ::pwtk::fileMustExist $head.$maxInd.out
    ::pwtk::fileMustExist $head.$iM.out
    
    # print some useful informations

    print_ "Image with the maximum energy taken from file: $head.$maxInd.out

The following images will be used as starting points for the NEB calculation:
   * first-image from        : $head.$im.out
   * intermediate-image from : $head.$maxInd.out
   * last-image from         : $head.$iM.out\n"
    
    # extract the images 

    for {set i $im} {$i <= $iM} {incr i} {
        set image($i) [::pwtk::pwo::getAtmPos -k $head.$i.out]
        set unit($i)  [::pwtk::trimCardFlag [lrange [lindex [split [string trim $image($i)] \n] 0] 1 end]]
        set coor($i)  [::pwtk::atmPosToCoor $image($i)]
    }

    # set the NEB input

    variable sbco

    set neb_head post-neb.[::pwtk::trim_prefix sbco $head]

    CONTROL " prefix  = '$neb_head' , calculation = "
    input_clear IONS CONSTRAINTS

    POSITIONS "
FIRST_IMAGE
$image($im)

INTERMEDIATE_IMAGE
ATOMIC_POSITIONS $unit($maxInd)
$coor($maxInd)

LAST_IMAGE
ATOMIC_POSITIONS $unit($iM)
$coor($iM)"

    ::pwtk::ifexist sbco(CI_scheme) {
        PATH " CI_scheme = [::pwtk::squote $sbco(CI_scheme)] "
    }
    set CI_scheme [::pwtk::input::namelistGetVarValue PATH CI_scheme trim]
    ::pwtk::ifset CI_scheme no-CI
    set input  $neb_head.$CI_scheme.in
    set output $neb_head.$CI_scheme.out
    
    # run the NEB calculation

    ::pwtk::outdir_postfix $neb_head
    ::pwtk::runNEB $input $output

    ::pwtk::print "Post-NEB SBCO calculation finished. See the \"$output\" file for details."

    # get the cumulative activation energy
    # (neb.x calculates activation energy with respect to the SBCO image No.$im)

    foreach {var value} [::pwtk::path::energies $neb_head.path] { set $var $value }
    # BEWARE: energies in the path file are in Hartrees
    set Eact_neb [expr $::pwtk::ha2ev*($Emax - $E(1))]

    set Eact_sbco [expr $ene_ev($maxInd) - $ene_ev(0)]
    set Eact_corr [expr $ene_ev($im) + $Eact_neb - $ene_ev(0)]

    ::pwtk::print "
Post-NEB SBCO summary:
   * SBCO image with the highest energy: $maxInd
   * SBCO only activation energy = [format %.4f $Eact_sbco] eV (measured wrt SBCO image No.0)

   * NEB image with the highest energy: $iTS
   * NEB only activation energy  = [format %.4f $Eact_neb] eV (measured wrt SBCO image No.$im)

   * NEB corrected SBCO activation energy = [format %.4f $Eact_corr] eV (measured wrt SBCO image No.0)
"    
    return $Eact_corr
}