geometric/ic_tools.py
def check_internal_hess(coords, molecule, IC, engine, dirname, verbose=0):
    """ Calculate the Cartesian Hessian using finite difference,
    transform to internal coordinates, then check the internal coordinate
    Hessian using finite difference. """
    # Initial energy and gradient
    spcalc = engine.calc(coords, dirname)
    E = spcalc['energy']
    gradx = spcalc['gradient']
    # Finite difference step
    h = 1.0e-3

    # Calculate Hessian using finite difference
    nc = len(coords)
    Hx = np.zeros((nc, nc), dtype=float)
    logger.info("Calculating Cartesian Hessian using finite difference on Cartesian gradients\n")
    for i in range(nc):
        logger.info(" coordinate %i/%i\n" % (i+1, nc))
        coords[i] += h
        gplus = engine.calc(coords, dirname)['gradient']
        coords[i] -= 2*h
        gminus = engine.calc(coords, dirname)['gradient']
        coords[i] += h
        Hx[i] = (gplus-gminus)/(2*h)

    # Internal coordinate Hessian using analytic transformation
    Hq = IC.calcHess(coords, gradx, Hx)

    # Initial internal coordinates and gradient
    q0 = IC.calculate(coords)
    Gq = IC.calcGrad(coords, gradx)

    Hq_f = np.zeros((len(q0), len(q0)), dtype=float)
    logger.info("-=# Now checking Hessian of the energy in internal coordinates using finite difference on gradient #=-\n")
    logger.info("%20s %20s : %14s %14s %14s\n" % ('IC1 Name', 'IC2 Name', 'Analytic', 'Numerical', 'Abs-Diff'))
    h = 1.0e-2
    for i in range(len(q0)):
        dq = np.zeros_like(q0)
        dq[i] += h
        x1 = IC.newCartesian(coords, dq, verbose)
        qplus = IC.calculate(x1)
        gplus = engine.calc(x1, dirname)['gradient']
        gqplus = IC.calcGrad(x1, gplus)
        dq[i] -= 2*h
        x1 = IC.newCartesian(coords, dq, verbose)
        qminus = IC.calculate(x1)
        gminus = engine.calc(x1, dirname)['gradient']
        gqminus = IC.calcGrad(x1, gminus)
        fdiffg = (gqplus-gqminus)/(2*h)
        for j in range(len(q0)):
            fdiff = fdiffg[j]
            Hq_f[i, j] = fdiff
            logger.info("%20s %20s : % 14.6f % 14.6f % 14.6f\n" % (IC.Internals[i], IC.Internals[j], Hq[i, j], fdiff, Hq[i, j]-fdiff))

    Eigx = sorted(np.linalg.eigh(Hx)[0])
    logger.info("Hessian Eigenvalues (Cartesian):\n")
    for i in range(len(Eigx)):
        logger.info("% 10.5f %s" % (Eigx[i], "\n" if i%9 == 8 else ""))
    Eigq = sorted(np.linalg.eigh(Hq)[0])
    logger.info("Hessian Eigenvalues (Internal):\n")
    for i in range(len(Eigq)):
        logger.info("% 10.5f %s" % (Eigq[i], "\n" if i%9 == 8 else ""))
    return Hq, Hq_f
 Calculate the Cartesian Hessian using finite difference,
    transform to internal coordinates, then check the internal coordinate
    Hessian using finite difference. 
Calculate the Cartesian Hessian using finite difference,
transform to internal coordinates, then check the internal coordinate
Hessian using finite difference.
Calculate the Cartesian Hessian using finite difference transform to internal coordinates then check the internal coordinate Hessian using finite difference .
def check_internal_hess(coords, molecule, IC, engine, dirname, verbose=0):
    spcalc = engine.calc(coords, dirname)
    E = spcalc['energy']
    gradx = spcalc['gradient']
    h = 1.0e-3
    nc = len(coords)
    Hx = np.zeros((nc, nc), dtype=float)
    logger.info("Calculating Cartesian Hessian using finite difference on Cartesian gradients\n")
    for i in range(nc):
        logger.info(" coordinate %i/%i\n" % (i+1, nc))
        coords[i] += h
        gplus = engine.calc(coords, dirname)['gradient']
        coords[i] -= 2*h
        gminus = engine.calc(coords, dirname)['gradient']
        coords[i] += h
        Hx[i] = (gplus-gminus)/(2*h)
    Hq = IC.calcHess(coords, gradx, Hx)
    q0 = IC.calculate(coords)
    Gq = IC.calcGrad(coords, gradx)
    Hq_f = np.zeros((len(q0), len(q0)), dtype=float)
    logger.info("-=# Now checking Hessian of the energy in internal coordinates using finite difference on gradient #=-\n")
    logger.info("%20s %20s : %14s %14s %14s\n" % ('IC1 Name', 'IC2 Name', 'Analytic', 'Numerical', 'Abs-Diff'))
    h = 1.0e-2
    for i in range(len(q0)):
        dq = np.zeros_like(q0)
        dq[i] += h
        x1 = IC.newCartesian(coords, dq, verbose)
        qplus = IC.calculate(x1)
        gplus = engine.calc(x1, dirname)['gradient']
        gqplus = IC.calcGrad(x1, gplus)
        dq[i] -= 2*h
        x1 = IC.newCartesian(coords, dq, verbose)
        qminus = IC.calculate(x1)
        gminus = engine.calc(x1, dirname)['gradient']
        gqminus = IC.calcGrad(x1, gminus)
        fdiffg = (gqplus-gqminus)/(2*h)
        for j in range(len(q0)):
            fdiff = fdiffg[j]
            Hq_f[i, j] = fdiff
            logger.info("%20s %20s : % 14.6f % 14.6f % 14.6f\n" % (IC.Internals[i], IC.Internals[j], Hq[i, j], fdiff, Hq[i, j]-fdiff))
    Eigx = sorted(np.linalg.eigh(Hx)[0])
    logger.info("Hessian Eigenvalues (Cartesian):\n")
    for i in range(len(Eigx)):
        logger.info("% 10.5f %s" % (Eigx[i], "\n" if i%9 == 8 else ""))
    Eigq = sorted(np.linalg.eigh(Hq)[0])
    logger.info("Hessian Eigenvalues (Internal):\n")
    for i in range(len(Eigq)):
        logger.info("% 10.5f %s" % (Eigq[i], "\n" if i%9 == 8 else ""))
    return Hq, Hq_f
def check_internal_hess ( coords , molecule , IC , engine , dirname , verbose = 0 ) : spcalc = engine . calc ( coords , dirname ) E = spcalc [ 'energy' ] gradx = spcalc [ 'gradient' ] h = 1.0e-3 nc = len ( coords ) Hx = np . zeros ( ( nc , nc ) , dtype = float ) logger . info ( "Calculating Cartesian Hessian using finite difference on Cartesian gradients\n" ) for i in range ( nc ) : logger . info ( " coordinate %i/%i\n" % ( i + 1 , nc ) ) coords [ i ] += h gplus = engine . calc ( coords , dirname ) [ 'gradient' ] coords [ i ] -= 2 * h gminus = engine . calc ( coords , dirname ) [ 'gradient' ] coords [ i ] += h Hx [ i ] = ( gplus - gminus ) / ( 2 * h ) Hq = IC . calcHess ( coords , gradx , Hx ) q0 = IC . calculate ( coords ) Gq = IC . calcGrad ( coords , gradx ) Hq_f = np . zeros ( ( len ( q0 ) , len ( q0 ) ) , dtype = float ) logger . info ( "-=# Now checking Hessian of the energy in internal coordinates using finite difference on gradient #=-\n" ) logger . info ( "%20s %20s : %14s %14s %14s\n" % ( 'IC1 Name' , 'IC2 Name' , 'Analytic' , 'Numerical' , 'Abs-Diff' ) ) h = 1.0e-2 for i in range ( len ( q0 ) ) : dq = np . zeros_like ( q0 ) dq [ i ] += h x1 = IC . newCartesian ( coords , dq , verbose ) qplus = IC . calculate ( x1 ) gplus = engine . calc ( x1 , dirname ) [ 'gradient' ] gqplus = IC . calcGrad ( x1 , gplus ) dq [ i ] -= 2 * h x1 = IC . newCartesian ( coords , dq , verbose ) qminus = IC . calculate ( x1 ) gminus = engine . calc ( x1 , dirname ) [ 'gradient' ] gqminus = IC . calcGrad ( x1 , gminus ) fdiffg = ( gqplus - gqminus ) / ( 2 * h ) for j in range ( len ( q0 ) ) : fdiff = fdiffg [ j ] Hq_f [ i , j ] = fdiff logger . info ( "%20s %20s : % 14.6f % 14.6f % 14.6f\n" % ( IC . Internals [ i ] , IC . Internals [ j ] , Hq [ i , j ] , fdiff , Hq [ i , j ] - fdiff ) ) Eigx = sorted ( np . linalg . eigh ( Hx ) [ 0 ] ) logger . info ( "Hessian Eigenvalues (Cartesian):\n" ) for i in range ( len ( Eigx ) ) : logger . info ( "% 10.5f %s" % ( Eigx [ i ] , "\n" if i % 9 == 8 else "" ) ) Eigq = sorted ( np . linalg . eigh ( Hq ) [ 0 ] ) logger . info ( "Hessian Eigenvalues (Internal):\n" ) for i in range ( len ( Eigq ) ) : logger . info ( "% 10.5f %s" % ( Eigq [ i ] , "\n" if i % 9 == 8 else "" ) ) return Hq , Hq_f
Calculate the Cartesian Hessian using finite difference,
transform to internal coordinates, then check the internal coordinate
Hessian using finite difference.
Calculate the Cartesian Hessian using finite difference transform to internal coordinates then check the internal coordinate Hessian using finite difference .
""" Calculate the Cartesian Hessian using finite difference,
    transform to internal coordinates, then check the internal coordinate
    Hessian using finite difference. """ # Initial energy and gradient # Finite difference step # Calculate Hessian using finite difference # Internal coordinate Hessian using analytic transformation # Initial internal coordinates and gradient
coords
molecule
IC
engine
dirname
verbose