#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
#! /usr/bin/env python
import sys
import string
import os
import math
import numpy as np
import CDTK.Tools.Conversion as conv
PeriodTable = {'H':1, 'He':2 ,
'Li':3, 'Be':4, 'B':5, 'C':6, 'N':7, 'O':8, 'F':9, 'Ne':10,
}
MOLCAS_MO_OUT = 'molcas_mo.out'
[docs]
class ObservableEH(object):
"""
Compute motion of the electronic hole in the molecule
"""
def __init__(self,**opts):
"""
Init an ObservableEH object
- atom_symbols: symbols of the chemical elements
default = None
- simbox: SimulationBox object
default = None
- external_AO: function returning the coefficients of atomic orbitals
in molecular orbital of current electronic state
default = None
- period_table: periodic table of chemical elements
- _ONAME: name of the observable [ob]servable_[e]lectronic[h]ole
"""
self.atom_symbols = opts.get('atom_symbols',None)
self.simbox = opts.get('simbox',None)
self.external_AO = opts.get('external_AO',None)
self.period_table = PeriodTable
self._ONAME = 'ob_eh'
[docs]
def getMapAO_Atom(self,**opts):
"""
Map the atomic orbitals to each atom
"""
try:
StateF = open(MOLCAS_MO_OUT,'r')
except:
msg = 'File '+MOLCAS_MO_OUT+' could not be accessed for reading'
raise ValueError(msg)
basis_set = None
Line=StateF.readline()
while Line:
if Line.find('Basis set label:') >= 0:
basis_set = Line.split('.')[1].upper()
break
Line=StateF.readline()
StateF.close()
if basis_set is None:
raise ValueError('Basis set could not be accessed for reading')
# construct the map for each basis set
# [3-21G] [6-31G] [6-31G*]
mapAO_Atom = []
AO_c = 0
for atom in self.atom_symbols:
AtomNum = self.period_table[atom] # atomic number of element
if basis_set == 'STO-3G':
if AtomNum >= 1 and AtomNum <= 2: # first period
# 1[1s]
AO_c += 1
mapAO_Atom.append(AO_c)
elif AtomNum >= 3 and AtomNum <= 8: # second period
# 1[1s] 1[2s] 3[2p]
AO_c += 5
mapAO_Atom.append(AO_c)
elif basis_set == '3-21G':
if AtomNum >= 1 and AtomNum <= 2: # first period
# 2[1s]
AO_c += 2
mapAO_Atom.append(AO_c)
elif AtomNum >= 3 and AtomNum <= 8: # second period
# 1[1s] 2[2s] 3x2[2p]
AO_c += 9
mapAO_Atom.append(AO_c)
else:
raise ValueError('Atomic number of element not defined')
elif basis_set == '6-31G':
if AtomNum >= 1 and AtomNum <= 2: # first period
# 2[1s]
AO_c += 2
mapAO_Atom.append(AO_c)
elif AtomNum >= 3 and AtomNum <= 8: # second period
# 1[1s] 2[2s] 3x2[2p]
AO_c += 9
mapAO_Atom.append(AO_c)
else:
raise ValueError('Atomic number of element not defined')
elif basis_set == '6-31G*':
if AtomNum >= 1 and AtomNum <= 2: # first period
# 2[1s]
AO_c += 2
mapAO_Atom.append(AO_c)
elif AtomNum >= 3 and AtomNum <= 8: # second period
# 1[1s] 2[2s] 3x2[[2p]+[d]]
AO_c += 15
mapAO_Atom.append(AO_c)
else:
raise ValueError('Atomic number of element not defined')
else:
raise ValueError('Basis set not defined')
return mapAO_Atom, AO_c
[docs]
def getElecHolePosition(self, a_x, state_current, **opts):
"""
Return the coordinate of the electronic hole
Inputs:
- a_x: molecular coordinates in bohr
- state_current: current electronic state
"""
X = a_x.copy()
natoms = X.size/3
X.shape = (natoms,3)
# run Molcas calculation for the molecular orbitals
if self.external_AO is None:
raise ValueError('external atomic orbital coefficients not defined')
# AO coefficients for Koopmans MO of [state_current]
self.AO_coeff = self.external_AO(a_x, state_current)
self.AO_coeff /= np.linalg.norm(self.AO_coeff) # normalization
# construct the AO to atom map
mapAO_Atom, NORB = self.getMapAO_Atom()
if len(self.AO_coeff) != NORB:
raise ValueError('Atomic orbital coefficients incorrect')
# The position of electronic hole
X_eh = np.zeros(3,float)
AO_previous_atom = 0
for i, AO_current_atom in enumerate(mapAO_Atom):
rho_AO_atom = np.linalg.norm(self.AO_coeff[AO_previous_atom:AO_current_atom])**2
X_eh += X[i] * rho_AO_atom
AO_previous_atom = AO_current_atom # AO of atom updated
# save the molecular orbital file
MOLDEN_MO = 'molcas_mo.scf.molden'
TIME = str(self.simbox.TIME)
STATE = str(state_current)
FileMO_TIME = MOLDEN_MO + '_' + TIME +'au_' + STATE
EXEC = 'cp ' + MOLDEN_MO + ' ' + FileMO_TIME
status = os.system(EXEC)
return X_eh
def _update_logfile(self):
if self.simbox._stepnum == 0:
f = self.simbox.logfile
f.write('Observable [electronic hole position] computed\n')
fN_vtf = self.simbox._NAME+'_h.vtf'
self.trajXeh2VTF(fileName=fN_vtf)
[docs]
def trajXeh2VTF(self,**opts):
"""
Generate a trajectory with position of the electronic hole in .vtf format
Useful for plots with VMD
"""
if self.simbox.internals: return
fileName = opts.get('fileName',self.simbox._NAME+'.vtf')
if not self.simbox.isrestart:
out = file(fileName,'w')
elif self.simbox.isrestart:
out = file(fileName,'w')
# symbols for the atoms
ncoords = self.simbox.X.size
nat = ncoords/self.simbox.ndim
for i,s in enumerate(self.simbox.atomSymbols):
out.write( 'atom %i radius 1.0 name %s\n' % (i,s) )
# symbol for the electronic hole
out.write( 'atom %i radius 1.0 name %s\n' % (i+1, 'X'))
out.write('\n')
# coordinates for the atoms
for i,geom in enumerate(self.simbox.TrajX):
out.write('timestep\n')
geom.shape = (nat,self.simbox.ndim)
for atom in geom:
for coord in atom:
out.write(' %12.7f' % (coord*conv.au2an))
out.write('\n')
geom.shape = (ncoords,)
# coordinates for the electronic hole
neh = 1 # one electronic hole
geom_eh = self.simbox.TrajXeh[i]
for coord in geom_eh:
out.write(' %12.7f' % (coord*conv.au2an))
out.write('\n')
out.write('\n')