#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
# Taylor expansion of the adiabatic potential energy surfaces at the Franck Condon point
# Methods for generating the coefficients and writing the MCTDH operator file
import numpy as np
import sys, getopt
import os
import re
import CDTK
import CDTK.Interfaces.GamessUSInterface as gi
import CDTK.Tools.NormalModeAnalysis as nm
import CDTK.Tools.Conversion as cv
[docs]
def cda_gradient(a,b,h):
"""
Gradient @ X0 by central difference approximation
a - f[X0 + h]
b - f[X0 - h]
h - step size
"""
return (a - b) / 2.0 / h
[docs]
def cda_curvature(a,b,c,d,e,h):
"""
Curvature @X0 by central difference approximation
a - f[X0 + 2h]
b - f[X0 + h]
c - f[X0]
d - f[X0 - h]
e - f[X0 - 2h]
h - step size
"""
return (-a + 16*b - 30*c + 16*d - e) / 12.0 / h / h
[docs]
def cda_crossterm(a,b,c,d,h):
"""
Cross term by central difference approximation
a - f[X0 + h e_i + h e_j]
b - f[X0 + h e_i - h e_j]
c - f[X0 - h e_i + h e_j]
d - f[X0 - h e_i - h e_j]
h - step size
"""
return (a - b - c + d) / 4.0 / h / h
[docs]
def calc_coeffs_diag(gs, Lx, X0, i, h, gap=None):
"""
Calculate the first and second Taylor coefficients for displacements along Lx[i]
gs - GamessUSInterface object for MCSCF calculation
Lx - Normal mode array with normal modes in rows
X0 - ground state equilibrium geometry [Angstrom]
i - Normal mode index
h - step size
gap - Energies of electronic states at the center as np array (optional)
"""
if gap is None:
gap = np.asarray(gs.mcscf(X0 * cv.an2au)[0])
a = np.asarray(gs.mcscf(nm_translate(X0, Lx[i], 2*h)*cv.an2au)[0])
b = np.asarray(gs.mcscf(nm_translate(X0, Lx[i], h)*cv.an2au)[0])
c = gap
d = np.asarray(gs.mcscf(nm_translate(X0, Lx[i], -h)*cv.an2au)[0])
e = np.asarray(gs.mcscf(nm_translate(X0, Lx[i], -2*h)*cv.an2au)[0])
return cda_gradient(b,d,h), cda_curvature(a,b,c,d,e,h)
[docs]
def calc_coeffs_diag_koopmans(gs, Lx, X0, i, h, gap=None):
"""
Calculate the first and second Taylor coefficients for displacements along Lx[i]
Koopmans level of theory
gs - GamessUSInterface object for CIS calculation
Lx - Normal mode array with normal modes in rows
X0 - ground state equilibrium geometry [Angstrom]
i - Normal mode index
h - step size
gap - Energies of electronic states at the center as np array (optional)
"""
if gap is None:
gap = np.asarray(gs.energy_HF_KOOPMANS(X0 * cv.an2au))
a = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[i], 2*h)*cv.an2au))
b = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[i], h)*cv.an2au))
c = gap
d = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[i], -h)*cv.an2au))
e = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[i], -2*h)*cv.an2au))
return cda_gradient(b,d,h), cda_curvature(a,b,c,d,e,h)
[docs]
def calc_coeffs_cross(gs, Lx, X0, i, j, h):
"""
Calculate the second order cross Taylor coefficient along Lx[i],Lx[j]
gs - GamessUSInterface object for MCSCF calculation
Lx - Normal mode array with normal modes in rows
X0 - ground state equilibrium geometry [Angstrom]
i, j - Normal mode indices
h - step size
"""
ixa = np.array([i,j])
a = np.asarray(gs.mcscf(nm_translate(X0, Lx[ixa], h*np.array([1,1]))*cv.an2au)[0])
b = np.asarray(gs.mcscf(nm_translate(X0, Lx[ixa], h*np.array([1,-1]))*cv.an2au)[0])
c = np.asarray(gs.mcscf(nm_translate(X0, Lx[ixa], h*np.array([-1,1]))*cv.an2au)[0])
d = np.asarray(gs.mcscf(nm_translate(X0, Lx[ixa], h*np.array([-1,-1]))*cv.an2au)[0])
return cda_crossterm(a,b,c,d,h)
[docs]
def calc_coeffs_cross_koopmans(gs, Lx, X0, i, j, h):
"""
Calculate the second order cross Taylor coefficient along Lx[i],Lx[j]
CI singles level of theory
gs - GamessUSInterface object for CIS calculation
Lx - Normal mode array with normal modes in rows
X0 - ground state equilibrium geometry [Angstrom]
i, j - Normal mode indices
h - step size
"""
ixa = np.array([i,j])
a = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[ixa], h*np.array([1,1]))*cv.an2au))
b = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[ixa], h*np.array([1,-1]))*cv.an2au))
c = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[ixa], h*np.array([-1,1]))*cv.an2au))
d = np.asarray(gs.energy_HF_KOOPMANS(nm_translate(X0, Lx[ixa], h*np.array([-1,-1]))*cv.an2au))
return cda_crossterm(a,b,c,d,h)
[docs]
def nm_translate(X0, Lx, dist):
"""
Calculate coordinates for translation along normal modes.
X0 - initial position
Lx - vector of cartesian normal modes
dist - vector of distances
"""
dist = np.asarray(dist)
dist = np.array(dist, copy=False, ndmin=1)
Lx = np.array(Lx, copy=False, ndmin=2)
if dist.size != Lx.shape[0]:
raise ValueError('Given distances did not match given modes.')
if Lx.shape[1] != X0.size:
raise ValueError('Dimensions of normal mode and initial position did not match.')
X = X0
for i in range(dist.size):
X = X + dist[i] * Lx[i]
return X
[docs]
def atomic_masses(atomicSymbols):
"""
Return atomic masses from periodic table for each coordinate [A.U.]
atomicSymbols - list of atomic symbols
"""
M = np.zeros(len(atomicSymbols)*3)
for i,asym in enumerate(atomicSymbols):
M[3*i:3*i+3] = cv.periodicTable[asym]['atomic_mass'] * cv.am2au
return np.asarray(M)
[docs]
def shift_energies_zero(e):
""""
Shift array of energies such that the lowest one is zero
e - array of energies
"""
return e - e.min()
[docs]
def coeffs_x_to_q_mfw(a, mu, w):
'''
Transform coefficient vector / matrix a from cartesian to mass and frequency weighted
coordinates
a.x to alpha.Q
x.g.x to Q.gamma.Q
a - object to be transformed
mu - reduced masses
w - frequencies
'''
CxQ = np.diag([np.sqrt(mu[i]*w[i]) for i in range(w.size)])
CQx = np.linalg.inv(CxQ)
if len(a.shape)==1: # vector
return np.dot(a, CQx)
elif len(a.shape)==2 and a.shape[0]==a.shape[1]: # matrix
return np.dot(CQx.T, np.dot(a, CQx))
else:
raise ValueError('a must be vector / square matrix')
[docs]
def calc_coefficients(gs, Lx, X0, h, save='', ix1=None, ix2=None, gap=None, grad=None, curv=None, verbose=False, mcscf=True):
"""
Calculate all coefficients for a Taylor series expansion at a given evaluation point
Gap, gradient (kappa_i) and curvatures (gamma_i,j) calculated in Cartesian
normal mode coordinates (cnmc)
gs - GamessUSInterface object with options set for MCSCF calculation
Lx - cartesian normal modes
X0 - evaluation point [Angstrom]
h - step size
save - save coefficients in given path after completing one outer loop (default: don't)
ix1 - pass row index (i) where the coefficients should be calculated (default: calculate all coefficients)
ix2 - pass column index (j) where the coefficients should be calculated (default: calculate all coefficients)
verbose - print runtime information if True (default: False)
gap,grad,curv - np arrays can be passed
mcscf - if TRUE, use MCSCF level of theory for calculating the PES. If FALSE, use Koopmans (default: TRUE)
"""
dof = Lx.shape[0]
if gap is None:
if mcscf:
gap = np.asarray(gs.mcscf(X0 * cv.an2au)[0])
else:
gap = np.asarray(gs.energy_HF_KOOPMANS(X0 * cv.an2au))
if verbose:
print('Finished calculating gap')
else:
if verbose:
print('Using provided values for gap')
if save:
np.save('cnmc_gap', gap)
if grad is None:
grad = np.zeros([dof, gs.nstate])
if curv is None:
curv = np.zeros([dof, dof, gs.nstate])
if ix1 is not None and ix2 is not None:
if ix1 == ix2:
if verbose:
print('Calculating diagonal coefficients: ' + str(ix1))
if mcscf:
grad[ix1,:],curv[ix1,ix1,:] = calc_coeffs_diag(gs, Lx, X0, ix1, h, gap=gap)
else:
grad[ix1,:],curv[ix1,ix1,:] = calc_coeffs_diag_koopmans(gs, Lx, X0, ix1, h, gap=gap)
if save:
np.save('cnmc_grad_{0:02d}'.format(ix1), grad[ix1,:])
np.save('cnmc_curv_{0:02d}_{1:02d}'.format(ix1, ix2), curv[ix1, ix1, :])
if verbose:
print('Created cnmc_grad_{0:02d}.npy'.format(ix1))
print('Created cnmc_curv_{0:02d}_{1:02d}'.format(ix1, ix2))
else:
if verbose:
print('Calculating cross terms: ' + str(ix1) + ', ' + str(ix2))
if mcscf:
curv[ix1,ix2,:] = calc_coeffs_cross(gs, Lx, X0, ix1, ix2, h)
else:
curv[ix1,ix2,:] = calc_coeffs_cross_koopmans(gs, Lx, X0, ix1, ix2, h)
curv[ix2,ix1,:] = curv[ix1,ix2,:]
if save:
np.save('cnmc_curv_{0:02d}_{1:02d}'.format(ix1, ix2), curv[ix1, ix2, :])
if verbose:
print('Created curv_{0:02d}_{1:02d}'.format(ix1, ix2))
else:
for i in range(dof):
if verbose:
print('Calculating diagonal coefficients: ' + str(i))
if mcscf:
grad[i,:],curv[i,i,:] = calc_coeffs_diag(gs, Lx, X0, i, h, gap=gap)
else:
grad[i,:],curv[i,i,:] = calc_coeffs_diag_koopmans(gs, Lx, X0, i, h, gap=gap)
for j in range(i+1,dof):
if verbose:
print('Calculating cross terms: ' + str(i) + ', ' + str(j))
if mcscf:
curv[i,j,:] = calc_coeffs_cross(gs, Lx, X0, i, j, h)
else:
curv[i,j,:] = calc_coeffs_cross_koopmans(gs, Lx, X0, i, j, h)
curv[j,i,:] = curv[i,j,:]
if save:
write_coefficients('TaylorAPESCoefficients{0:02d}.dat'.format(i+1),\
grad,curv,i,comment='Cartesian coordinates (x)')
return gap, grad, curv
[docs]
def write_coefficients(fname, grad, curv, ix, comment=''):
"""
Write the coefficients given in grad, curv [ix] to fname
"""
nstate = grad.shape[-1]
grad_template = ' grad{elix}_{ix:02d}\t= {v:.8f}\n'
curv_template = ' curv{elix}_{ix1:02d}_{ix2:02d}\t= {v:.8f}\n'
with open(fname, 'w') as f:
f.write(comment + '\n')
for i in range(nstate):
f.write(grad_template.format(elix=i+1,ix=ix+1,v=grad[ix,i]))
f.write(curv_template.format(elix=i+1,ix1=ix+1,ix2=ix+1,v=curv[ix,ix,i]))
for j in range(ix+1,grad.shape[0]):
f.write(curv_template.format(elix=i+1,ix1=ix+1,ix2=j+1,v=curv[ix,j,i]))
[docs]
def write_mctdh_op(opfile, freq, gap, grad, curv, title):
"""
Write an MCTDH operator file for Taylor model Hamiltonian
opfile - target file name
freq - eigenfrequencies of normal modes
gap - for each electronic state, energy gap (scalar)
grad - for each electronic state, gradient (vector)
curv - for each electronic state, curvature (matrix)
"""
dof = freq.size
nstate = gap.size
# Line templates
ef_template = ' w{ix:02d}\t= {v:.8f}\n'
ml_template = ' L{ix:02d} |'
egap_template = ' egap{elix}\t= {v:.8f}\n'
grad_template = ' grad{elix}_{ix:02d}\t= {v:.8f}\n'
curv_template = ' curv{elix}_{ix1:02d}_{ix2:02d}\t= {v:.8f}\n'
op_ke_template = ' -0.5*w{ix:02d}\t|{mix}\tdq^2\n'
op_egap_template = ' egap{elix}\t|1 S{elix}&{elix}\n'
op_grad_template = ' grad{elix}_{ix:02d}\t|1 S{elix}&{elix}\t|{mix}\tq\n'
op_curv_template = ' 0.5*curv{elix}_{ix1:02d}_{ix1:02d}\t|1 S{elix}&{elix}\t' + '|{mix1:02d}\tq^2\n'
op_cross_template = ' curv{elix}_{ix1:02d}_{ix2:02d}\t|1 S{elix}&{elix}\t' + '|{mix1:02d}\tq\t|{mix2:02d}\tq\n'
tc_string = ''
ml_string = ''
ef_string = ''
op_ke_string = ''
op_pot_string = ''
# Fill in te mplates
for i,w in enumerate(freq):
ef_string = ef_string + ef_template.format(ix=i+1, v=w)
ml_string = ml_string + ml_template.format(ix=i+1)
op_ke_string = op_ke_string + op_ke_template.format(ix=i+1, mix=i+2)
# Fill in Taylor series coefficients
for i in range(nstate):
tc_string = tc_string + egap_template.format(elix=i+1, v=gap[i])
for j in range(dof):
tc_string = tc_string + grad_template.format(elix=i+1, ix=j+1, v=grad[j,i])
for k in range(j,dof):
tc_string = tc_string + curv_template.format(elix=i+1, ix1=j+1, ix2=k+1, v=curv[j,k,i])
# Fill in gap terms
for i in range(nstate):
op_pot_string = op_pot_string + op_egap_template.format(elix=i+1)
# Fill in linear terms
for i in range(nstate):
for j in range(dof):
op_pot_string = op_pot_string + op_grad_template.format(elix=i+1, ix=j+1, mix=j+2)
# Fill in quadratic terms
for i in range(nstate):
for j in range(dof):
op_pot_string = op_pot_string + op_curv_template.format(elix=i+1, ix1=j+1, mix1=j+2)
# Fill in cross terms
for i in range(nstate):
for j in range(dof):
for k in range(j+1, dof):
op_pot_string = op_pot_string + op_cross_template.format(elix=i+1, \
ix1=j+1, ix2=k+1, mix1=j+2, mix2=k+2)
context = {'title':title,
'eigenfreqs':ef_string,
'mode-labels':ml_string,
'taylorcoeffs':tc_string,
'kinetic':op_ke_string,
'potential':op_pot_string}
with open(opfile, 'w') as f:
path = CDTK.CDTKPATH + '/CDTK/Tools/'
f.write(open(path + 'TaylorAPES_MCTDH_OP_TEMPLATE').read().format(**context))
[docs]
def parse_gamess_options(path, p, func=[]):
"""
Parse GAMESS options from a text file
Format: parameter=value
"""
try:
gso = open(path + '/gamess_options.dat').read().split('\n')
except:
return None
for g in gso:
if re.search(p,g):
g = g.split('=')
if func:
return func(g[-1])
return g[-1]
return None
[docs]
def main(argv):
"""
python TaylorAPES.py --opfile=opfile, --title=title, --step=step, --inputpath=inputpath
opfile - name of the MCTDH operator file
title - MCTDH operator title
step - step size in coefficient calculation
inputpath - path with gamess options and initial geometry
Calculate Taylor coefficients with gamess options given in inputpath
X0 (equilibrium geometry) and Hessian (cartesian) may be given in inputpath.
Else, will be calculated from X0_guess.dat (xyz initial coordinates in
inputpath). Saves MCTDH operator file in inputpath/opfile.
"""
opts,args = getopt.getopt(argv, 'o:t:s:i:', ['opfile=', 'title=', 'step=', 'inputpath='])
# Default values for command line options
path = '.'
h = 0.01
title = 'Untitled operator'
opfile = 'untitled.op'
for o,a in opts:
if o in ('-o', '--opfile'):
opfile = a
elif o in ('-t', '--title'):
title = a
elif o in ('-i', '--inputpath'):
path = a
elif o in ('-s', '--step'):
h = float(a)
gs = gi.gamess()
gs.atomicSymbols = parse_gamess_options(path,'atomicSymbols',eval)
gs.atomicNumbers = parse_gamess_options(path,'atomicNumbers',eval)
gs.ngauss = parse_gamess_options(path,'ngauss',int)
gs.gbasis = parse_gamess_options(path,'gbasis',eval)
# Optimize geometry
if os.path.isfile(path + '/X0.dat'):
print('Using geometry from ' + path + '/X0.dat')
X0 = np.loadtxt(path + '/X0.dat')
else:
X0_guess = np.loadtxt(path + '/X0_guess.dat')
X0 = gs.optimize(X0_guess * cv.an2au)
np.savetxt(path + '/X0.dat', X0)
# Normal modes
if os.path.isfile(path + '/hessianCart.tab'):
print('Using Hessian from ' + path + '/hessianCart.tab')
H = np.loadtxt(path + '/hessianCart.tab')
else:
H = gs.hessian(X0*cv.an2au)
np.savetxt(path + '/hessianCart.tab', H)
M = atomic_masses(gs.atomicSymbols)
w,L = nm.normal_modes_from_HessianCart(H,M)
Lx,mu = nm.q_mw_to_x(M,L)
N_atoms = X0.size/3
w = w[-3*N_atoms+6:]
L = L[-3*N_atoms+6:]
Lx = Lx[-3*N_atoms+6:]
mu = mu[-3*N_atoms+6:]
np.savetxt(path + '/w.dat', w)
np.savetxt(path + '/Lx.dat', Lx)
np.savetxt(path + '/mu.dat', mu)
print('Frequencies [cm-1]')
for ww in w:
print(int(ww*cv.au2ic))
# GAMESS object for MCSCF energy calculation
gs2 = gi.gamess()
gs2.atomicSymbols = parse_gamess_options(path,'atomicSymbols',eval)
gs2.atomicNumbers = parse_gamess_options(path,'atomicNumbers',eval)
gs2.icharg = parse_gamess_options(path,'icharg',int)
gs2.ngauss = parse_gamess_options(path,'ngauss',int)
gs2.gbasis = parse_gamess_options(path,'gbasis',eval)
gs2.mult = parse_gamess_options(path,'mult',int)
gs2.sz = parse_gamess_options(path,'sz',float)
gs2.nels = parse_gamess_options(path,'nels',int)
gs2.ncore = parse_gamess_options(path,'ncore',int)
gs2.nact = parse_gamess_options(path,'nact',int)
gs2.nstate = parse_gamess_options(path,'nstate',int)
gs2.wstate = parse_gamess_options(path,'wstate',eval)
print('Running GAMESS MCSCSF energy calculation with paramters:')
print('atomicSymbols=' + str(gs2.atomicSymbols))
print('atomicNumbers=' + str(gs2.atomicNumbers))
print('icharg=' + str(gs2.icharg))
print('mult=' + str(gs2.mult))
print('sz=' + str(gs2.sz))
print('ngauss=' + str(gs2.ngauss))
print('gbasis=' + str(gs2.gbasis))
print('nels=' + str(gs2.nels))
print('ncore=' + str(gs2.ncore))
print('nact=' + str(gs2.nact))
print('wstate=' + gs2.wstate)
# Calculate coefficients
print('Calculating coefficients ...')
gap, grad, curv = calc_coefficients(gs2, Lx, X0, h, save=path)
gap_check = np.array(gs2.mcscf(X0*cv.an2au)[0])
# Check if calculation in the subroutines is done with correct units etc
assert np.allclose(gap, gap_check)
# Transform to mass and frequency weighted coordinates
grad_mw = np.zeros_like(grad)
curv_mw = np.zeros_like(curv)
for i in range(grad_mw.shape[-1]):
grad_mw[:,i] = coeffs_x_to_q_mfw(grad[:,i],mu,w)
curv_mw[:,:,i] = coeffs_x_to_q_mfw(curv[:,:,i],mu,w)
gap_shifted = shift_energies_zero(gap)
write_mctdh_op(opfile,w,gap_shifted,grad_mw,curv_mw,title)
if __name__=='__main__':
main(sys.argv[1:])