#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
"""
VibronicCouplingModel
"""
# TODO
# Analysis: coherence, tr[rho^2], parameter sets, time
import os
import re
import numpy as np
import datetime
import shutil
import tarfile
TINY = 1.e-15 # smaller numbers considered to be zero
MAX_GRID_EDGE_AMPLITUDE = 1.e-3 # Maximum tolerated amplitude at final grid pts
MAX_NATURAL_WEIGHT = 1.e-4 # Maximum natural weight at highest SPF
[docs]
class VibronicCouplingModel(object):
def __init__(self, **opts):
# Operator parameters
self.w1 = opts.get('w1', 0.01)
self.w2 = opts.get('w2', 0.01)
self.gap = opts.get('gap', 0.0)
self.off_gap = opts.get('off_gap', 0.0)
# linear coupling, intra-state
self.kappa_1_1 = opts.get('kappa_1_1', 0.0) # el state index _ mode index
self.kappa_1_2 = opts.get('kappa_1_2', 0.0)
self.kappa_2_1 = opts.get('kappa_2_1', 0.0)
self.kappa_2_2 = opts.get('kappa_2_2', 0.0)
# linear coupling, inter-state
self.lambda_1 = opts.get('lambda_1', 0.0)
self.lambda_2 = opts.get('lambda_2', 0.0)
# quadratic coupling, intra-state
self.gamma_1_11 = opts.get('gamma_1_11', 0.0)
self.gamma_1_12 = opts.get('gamma_1_12', 0.0)
self.gamma_1_22 = opts.get('gamma_1_22', 0.0)
self.gamma_2_11 = opts.get('gamma_2_11', 0.0)
self.gamma_2_12 = opts.get('gamma_1_12', 0.0)
self.gamma_2_22 = opts.get('gamma_2_22', 0.0)
# quadratic coupling, inter-state
self.mu_12 = opts.get('mu_12', 0.0)
# Initial center of the gaussian wavepacket
self.q1start = opts.get('q1start', 0.0)
self.q2start = opts.get('q2start', 0.0)
# Input parameters
self.name = opts.get('name', 'vibronic_coupling') # MCTDH target dir
self.opname = opts.get('opname', 'vibronic_coupling') # MCTDH op name
self.pes = opts.get('pes', False) # Write PES file
self.pq1n = opts.get('pq1n', 100) # Primitive basis functions Q1
self.pq2n = opts.get('pq2n', 100) # Primitive basis functions Q2
self.pq1xi = opts.get('pq1xi', -5.0) # xi Q1
self.pq1xf = opts.get('pq1xf', 5.0) # xf Q1
self.pq2xi = opts.get('pq2xi', -5.0) # xi Q2
self.pq2xf = opts.get('pq2xf', 5.0) # xi Q2
self.tfinal = opts.get('tfinal', 10.0) # time [fs]
self.tout = opts.get('tout', 0.05)
self.tpsi = opts.get('tpsi', 0.1)
self.spfq1s1 = opts.get('spfq1s1', 1) # SPF multi-set functions
self.spfq1s2 = opts.get('spfq1s2', 1)
self.spfq2s1 = opts.get('spfq2s1', 1)
self.spfq2s2 = opts.get('spfq2s2', 1)
self.acoeff1 = opts.get('acoeff1',1.0)
self.acoeff2 = opts.get('acoeff2',1.0)
# Run and post run files
self.results_dir = opts.get('results_dir','')
self.keep_mctdh_data = opts.get('keep_mctdh_data', False)
self._inputfile = 'vibronic_coupling.inp'
self._rdgpop = None
# Support for additional harmonic oscillator modes
self.multid = opts.get('multid', False)
self.multid_w = opts.get('multid_w', None)
self.multid_grad = opts.get('multid_grad', None)
self.multid_curv = opts.get('multid_curv', None)
self.columntitle = 'COLUMNTITLE'
self.verbose = opts.get('verbose', False)
[docs]
def spatial_resolution(self):
"""
Returns the spatial resolution for each mode
"""
return ((self.pq1xf-self.pq1xi)/self.pq1n,\
(self.pq2xf-self.pq2xi)/self.pq2n)
[docs]
def distance_min(self):
"""
Returns the distance of the HO minima for each mode
"""
return np.array([(self.kappa_2_1 - self.kappa_1_1)/self.w1,\
(self.kappa_2_2 - self.kappa_1_2)/self.w2])
[docs]
def set_spf(self, spf):
"""
Set all SPF numbers to spf
"""
self.spfq1s1 = spf
self.spfq1s2 = spf
self.spfq2s1 = spf
self.spfq2s2 = spf
[docs]
def set_xixfn(self, x, n):
"""
Set xi to -x, xf to x for all modes. Set # DVR to n for all modes.
"""
self.pq1xi = -x
self.pq1xf = x
self.pq2xi = -x
self.pq2xf = x
self.pq1n = n
self.pq2n = n
[docs]
def write_op(self):
"""
Fill in the MCTDH operator template for the vibronic coupling
model. Write file.
"""
w_line_template = ' {w{ix:02d} = {v:.6f}\n'
kinetic_line_template = ' w{ix:02d} |{mode_ix:d} dq^2\n'
potential_line_template = ' w{ix:02d} |{mode_ix:d} q^2\n'
label_line_template = ' | L{ix:02d} '
# TODO gradient and mode-mode couplings in additional modes
w_template = ''
kinetic_template = ''
potential_template = ''
label_template = ''
if self.multid:
for i,v in enumerate(self.multid_w):
print(i)
print(v)
label_template = label_template \
+ label_line_template.format(ix=i+1)
w_template = w_template + w_line_template.format(ix=i+1, v=v)
kinetic_template = kinetic_template \
+ kinetic_line_template.format(ix=i+1,mode_ix=i+4)
potential_template = potential_template \
+ potential_line_template.format(ix=i+1,mode_ix=i+4)
t = open(self.opname + '.op.template').read()
with open(self.opname + '.op', 'w') as f:
f.write(t.format(w1 = self.w1,
w2 = self.w2,
gap = self.gap,
off_gap = self.off_gap,
kappa_1_1 = self.kappa_1_1,
kappa_1_2 = self.kappa_1_2,
kappa_2_1 = self.kappa_2_1,
kappa_2_2 = self.kappa_2_2,
lambda_1 = self.lambda_1,
lambda_2 = self.lambda_2,
gamma_1_11 = self.gamma_1_11,
gamma_1_12 = self.gamma_1_12,
gamma_1_22 = self.gamma_1_22,
gamma_2_11 = self.gamma_2_11,
gamma_2_12 = self.gamma_2_12,
gamma_2_22 = self.gamma_2_22,
mu_12 = self.mu_12,
additional_dof_w = w_template,
additional_dof_mode_labels = label_template,
additional_dof_kinetic = kinetic_template,
additional_dof_potential = potential_template
))
[docs]
def write_inp(self):
"""
Fill in the MCTDH input template for the vibronic coupling model.
Write file.
"""
spf_line_template = ' L{ix:02d} = 1,1\n'
pbs_line_template = ' L{ix:02d} HO 61 xi-xf -3.0 3.0\n'
initwf_line_template = ' L{ix:02d} HO 0.0 0.0 1.0\n'
spf_template = ''
pbs_template = ''
initwf_template = ''
if self.multid:
for i in range(self.multid_w.shape[0]):
spf_template = spf_template + spf_line_template.format(ix=i+1)
pbs_template = pbs_template + pbs_line_template.format(ix=i+1)
initwf_template = initwf_template + initwf_line_template.format(ix=i+1)
t = open('vibronic_coupling.inp.template').read()
with open(self._inputfile, 'w') as f:
f.write(t.format(name = self.name,
opname = self.opname,
pq1n = self.pq1n,
pq2n = self.pq2n,
pq1xi = self.pq1xi,
pq1xf = self.pq1xf,
pq2xi = self.pq2xi,
pq2xf = self.pq2xf,
tfinal = self.tfinal,
tout = self.tout,
tpsi = self.tpsi,
spfq1s1 = self.spfq1s1,
spfq1s2 = self.spfq1s2,
spfq2s1 = self.spfq2s1,
spfq2s2 = self.spfq2s2,
q1start = self.q1start,
q2start = self.q2start,
acoeff1real = self.acoeff1.real,
acoeff1imag = self.acoeff1.imag,
acoeff2real = self.acoeff2.real,
acoeff2imag = self.acoeff2.imag,
additional_dof_spf = spf_template,
additional_dof_pbs = pbs_template,
additional_dof_initwf = initwf_template
))
[docs]
def write_results(self, t_fs, dat, fname):
"""
Write results file for use in GNUPLOT etc
print parameters
print t - dat to fname
"""
try:
f = open(fname, 'w')
except:
print('Failed to open + ' + fname)
return
f.write('# Results from MCTDH run\n')
f.write('# Vibronic coupling model Hamiltonian\n')
f.write(self.print_parameters())
f.write('t[fs]\t' + self.columntitle + '\n')
for i in range(t_fs.size):
f.write('{0:.4f}\t{1:.4f}\n'.format(t_fs[i], dat[i]))
f.close()
[docs]
def print_parameters(self):
"""
Return formatted string with all model parameters
"""
t = '''
# Current parameter set:
# w1={w1}, w2={w2}
# kappa_1 = ({kappa_1_1}, {kappa_1_2}),
# kappa_2 = ({kappa_2_1}, {kappa_2_2}),
# gamma_1 = (({gamma_1_11}, {gamma_1_12}),({gamma_1_12}, {gamma_1_22}))
# gamma_2 = (({gamma_2_11}, {gamma_2_12}),({gamma_2_12}, {gamma_2_22}))
# gap = {gap}, off_gap = {off_gap}
# lambda_1 = {lambda_1}, lambda_2 = {lambda_2}, mu_12 = {mu_12}
# -----------------
# {multi_d}
# -----------------
# A-coeff: ({acoeff1real}, {acoeff1imag}) ; ({acoeff2real}, {acoeff2imag})
# Q1: {pq1n} xi-xf {pq1xi} {pq1xf} | SPF: {spfq1s1},{spfq1s2}
# Q2: {pq2n} xi-xf {pq2xi} {pq2xf} | SPF: {spfq2s1},{spfq2s2}
'''
if self.multid:
multi_d = 'MULTI-D'
else:
multi_d = 'No additional DOF'
return t.format(w1=self.w1,
w2=self.w2,
kappa_1_1=self.kappa_1_1,
kappa_1_2=self.kappa_1_2,
kappa_2_1=self.kappa_2_1,
kappa_2_2=self.kappa_2_2,
gamma_1_11=self.gamma_1_11,
gamma_1_12=self.gamma_1_12,
gamma_1_22=self.gamma_1_22,
gamma_2_11=self.gamma_2_11,
gamma_2_12=self.gamma_2_12,
gamma_2_22=self.gamma_2_22,
gap=self.gap,
off_gap=self.off_gap,
lambda_1=self.lambda_1,
lambda_2=self.lambda_2,
mu_12=self.mu_12,
multi_d = multi_d,
acoeff1real=self.acoeff1.real,
acoeff1imag=self.acoeff1.imag,
acoeff2real=self.acoeff2.real,
acoeff2imag=self.acoeff2.imag,
pq1n=self.pq1n,
pq2n=self.pq2n,
pq1xi=self.pq1xi,
pq1xf=self.pq1xf,
pq2xi=self.pq2xi,
pq2xf=self.pq2xf,
spfq1s1=self.spfq1s1,
spfq1s2=self.spfq1s2,
spfq2s1=self.spfq2s1,
spfq2s2=self.spfq2s2
)
[docs]
def run_mctdh(self):
"""
Run MCTDH. Make name directory if it does not exist. Overwrite name
directory if it exists.
"""
# TODO mnd
self.write_inp()
self.write_op()
os.system('mctdh85 -w ' + self._inputfile)
# Analyze rdgpop
cwd = os.getcwd()
os.chdir(cwd + '/' + self.name)
os.system('rdgpop85 2 0 > rdgpop.dat')
os.chdir(cwd)
if not rdgpop_OK(self.name + '/rdgpop.dat'):
if self.verbose:
print(self.print_parameters())
raise ValueError('WARNING: rdgpop failed')
# Analyze rdcheck
# Single SPF allowed
if self.spfq1s1==1 and self.spfq1s2==1 and self.spfq2s1==1 and self.spfq2s2==1:
print('Skipped rdcheck: 1 SPF for each mode')
else:
cwd = os.getcwd()
os.chdir(cwd + '/' + self.name)
os.system('rdcheck85 0 0 > rdcheck.dat')
os.chdir(cwd)
if not rdcheck_OK(self.name + '/rdcheck.dat'):
if self.verbose:
print(self.print_parameters())
raise ValueError('WARNING: rdcheck failed')
# Calculate coherence and tr[rho^2]
t_fs, nrm, d11, d12, d22 = get_dm_t(self.name + '/expectation')
dm = rho(d11,d12,d22)
coh=coherence(dm)
tr_sq = trace_squared(dm)
# Target directory
now = datetime.datetime.now()
rd = self.results_dir + '/{0:4d}-{1:02d}-{2:02d}_{3:02d}-{4:02d}-{5:02d}'.format(now.year, \
now.month, now.day, now.hour, now.minute, now.second)
os.makedirs(rd)
self.write_results(t_fs, coh, rd + '/coherence.dat')
self.write_results(t_fs, tr_sq, rd + '/tr_sq.dat')
print('Results written to ' + rd)
# Calculate potential energy surface if set
if self.pes:
os.system('mctdh85 -pes -w ' + self._inputfile)
# Move all MCTDH runtime data to target directory
if self.keep_mctdh_data:
with tarfile.open(rd + '/mctdh_runtime_data.tar.gz','w:gz') as tf:
for mf in os.listdir(self.name):
tf.add(self.name + '/' + mf)
################################################################################
### TOOLS ETC.
################################################################################
[docs]
def get_dm_t(fname):
"""
Get density matrix elements from expectation file
Expected format: t_fs Norm DM11 DM12 DM22
"""
try:
data = np.genfromtxt(fname)
except:
print('File does not exist: ' + fname)
return
t_fs = data[:,0]
nrm = data[:,1]
d11 = data[:,2] + 1j*data[:,3]
d12 = data[:,4] + 1j*data[:,5]
d22 = data[:,6] + 1j*data[:,7]
return t_fs, nrm, d11, d12, d22
[docs]
def rdgpop_OK(fname):
"""
Analyze output of 'rdgpop85 <n> <f>'
Return True if the density at the edge of the position grid / the final edge
of the basis grid does not exceed MAX_GRID_EDGE_AMPLITUDE
"""
try:
f = open(fname)
except:
print('Could not open ' + fname)
return None
patt1 = r'^\s*#\s*dof\s*grid'
patt2 = r'^-----'
grid_begin = []
grid_end = []
basis_begin = []
basis_end = []
line = f.readline()
while line:
if re.search(patt1,line):
line = f.readline()
while not re.search(patt2,line):
l = map(float, line.split()[2:])
grid_begin.append(l[0])
grid_end.append(l[1])
basis_begin.append(l[2])
basis_end.append(l[3])
line = f.readline()
line = f.readline()
f.close()
if np.any(np.asarray(grid_begin) > MAX_GRID_EDGE_AMPLITUDE):
print('*** grid begin')
return False
if np.any(np.asarray(grid_end) > MAX_GRID_EDGE_AMPLITUDE):
print('*** grid end')
return False
if np.any(np.asarray(basis_end) > MAX_GRID_EDGE_AMPLITUDE):
print('*** basis end')
return False
return True
[docs]
def rdcheck_OK(fname):
"""
Analyze output of 'rdcheck85 0 0'
Return True if the natural weight of the highest SPF does not exceed
MAX_NATURAL_WEIGHT
"""
try:
f = open(fname)
except:
raise ValueError('Could not open ' + fname)
patt1 = r'^ Maximum over time of lowest nat\.-weight'
patt2a = r'^\s+3' # only consider modes 1 + 2
patt2b = r'^------' # if only 2 modes given
nat_weight = []
line = f.readline()
while line:
if re.search(patt1,line):
line = f.readline() # skip the line after patt1
line = f.readline()
while not (re.search(patt2a,line) or re.search(patt2b, line)):
l = map(float, line.split()[1:])
[nat_weight.append(ll) for ll in l]
line = f.readline()
line = f.readline()
f.close()
if np.any(np.asarray(nat_weight) > MAX_NATURAL_WEIGHT):
print('*** rdcheck')
return False
return True
[docs]
def rho(dm11, dm12, dm22):
"""
Construct time dependent 2x2 density matrix from elements
"""
rho = np.zeros([dm11.size, 2, 2], dtype='complex')
for i in range(dm11.size):
rho[i,:,:] = np.array([[dm11[i], dm12[i]],\
[np.conjugate(dm12[i]),dm22[i]]])
return rho
[docs]
def coherence(rho):
"""
Calculate the coherence from a time-dependent 2x2 density matrix rho
rho[i, :, :] = rho(t[i])
Returns |rho12| / sqrt(rho11 * rho22)
"""
coh = np.zeros(rho.shape[0])
for i in range(rho.shape[0]):
if np.any(np.isclose(rho[i,:,:], 0.0)):
coh[i] = 0.0
else:
coh[i] = np.abs(rho[i,0,1]) / np.sqrt(np.abs(rho[i,0,0]) * np.abs(rho[i,1,1]))
return coh
[docs]
def trace_squared(rho):
"""
Calculate the squared trace from a time-dependent 2x2 density matrix rho
rho[i, :, :] = rho(t[i])
Returns Tr[rho^2]
"""
tr_sq = np.zeros(rho.shape[0])
for i in range(rho.shape[0]):
tt = np.trace(np.dot(rho[i],rho[i]))
if not np.allclose(tt.imag, 0.0):
raise ValueError
tr_sq[i] = tt.real
return tr_sq