#* **************************************************************************
#*
#* 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 os
import sys
import datetime
import uuid
import math as mt
import numpy as np
import CDTK.Tools.Conversion as conv
[docs]
class ElectronicState(object):
"""
This class describes different ionized states used in MD simulations including "large jumps" using
a Monte Carlo scheme.
"""
def __init__(self,**opts):
self.state = opts.get('state', None) # initial state
self.dt = float(opts.get('dte', None)) # time step for integrating electronic structure
self.t = 0.0 # initial time
self.QCE = None # quantum chemistry package
self.field = None # pulse
self.MC_dt_ratio = float(opts.get('ratio', 0.1)) # desired ratio for variable step size scheme.
self.DIR = None # Directory for outputting data
self.stop = opts.get('stop', False) # stop after 1st hop?
# logs: charge - charge of system, partial - partial charges of all atoms in the system
# events - electronic dynamics events, pulse - pulse profile
self.log = {'charge': [], 'partial': [], 'event': [], 'pulse': [], 'occupation': []}
self.pdelay=opts.get('pdelay',None)
[docs]
def initLogs(self):
"""
Initialize logs.
"""
# initialize logs
partialCharges = self.QCE.getPartialCharges()
self.log['partial'].append((self.t,) + tuple(partialCharges))
self.log['charge'].append((self.t, int(round(np.sum(partialCharges)))))
if len(self.field) > 0:
self.log['pulse'].append((self.t, self.field[0].Envelope(self.t)))
self.log['occupation'].append((self.t, self.state))
[docs]
def step(self, dtN):
"""
This functions evolves the electronic state along a full nuclear step dtN.
Input: dtN - nuclear timestep
"""
# obtain rates from Quantum chemistry
gamma = self.QCE.getGamma()
# set final time
final_t = self.t + dtN
# report on hopping
hopped = False
if self.pdelay is not None and self.t <= self.pdelay and final_t > self.pdelay:
total_cs=0.
for p in gamma:
if p[0].lower() == 'p':
total_cs+=p[2]
while not hopped:
r = np.random.random()
for p in gamma:
if p[0].lower() != 'p':
continue
if(r<=p[2]/total_cs):
self.log['event'].append((self.t,) + p)
print("Happened: (triggered):", p, " at ", self.t)
sys.stdout.flush()
hopped = True
# to stop the calculation
if self.stop:
os.system("rm stop")
# set new state, etc.
self.state = p[1]
self.QCE.set_occ(self.state)
gamma = self.QCE.getGamma()
break
else:
r=r-p[2]/total_cs
# while not hopped
# integrate to full nuclear step
while self.t < final_t:
# careful to not overstep!
dt = min(self.dt, final_t - self.t)
# obtain total rate and adjust timestep
total_rate = 0.0
for p in gamma:
total_rate += self.rate(p)
if total_rate > 0.0:
dt = min(self.MC_dt_ratio / total_rate, dt)
# random number for MC scheme
r = np.random.random()
# go through all possible processes obtained from electronic structure code
for p in gamma:
# obtain rate and hopping probability for particular process
rate = self.rate(p)
prob = dt * rate
# hop occurs for this process
if(r < prob):
# log event
self.log['event'].append((self.t,) + p)
print("Happened: ", p, " at ", self.t)
sys.stdout.flush()
hopped = True
# to stop the calculation
if self.stop:
os.system("rm stop")
# set new state, etc.
self.state = p[1]
# get new rates
self.QCE.set_occ(self.state)
gamma = self.QCE.getGamma()
break
# no hop occurs for this process
else:
r = r-prob
if hopped and self.stop:
break
self.t += dt
# update logs and write
partialCharges = self.QCE.getPartialCharges()
self.log['partial'].append((self.t,) + tuple(partialCharges))
self.log['charge'].append((self.t, int(round(np.sum(partialCharges)))))
if len(self.field) > 0:
self.log['pulse'].append((self.t, self.field[0].Envelope(self.t)))
self.log['occupation'].append((self.t, self.state))
return hopped
[docs]
def rate(self, p):
"""
Obtain rate for current process.
Input:
p = process as obtained from the interface
tupel with: 0: type (currently: 'p' photoionization, 'f' fluouresence, 'a' auger
1: final state for this process
2: rate or cross section
3: energy of corresponding electron
Output:
rate for this process
"""
# for photo-ionization we get a cross section
# else we get the rate directly
if p[0].lower() == 'p':
if len(self.field) > 0:
return p[2] * self.field[0].Envelope(self.t)
else:
return 0.
else:
return p[2]
[docs]
def write_logs(self):
"""
Wrapper-Function to write logs. Set here which logs to write.
"""
self.write_log('charge')
self.write_log('partial')
self.write_log('event')
self.write_log('pulse')
self.write_log('occupation')
[docs]
def write_log(self, logtype):
"""
Function writing logs. Called by wrapper-function 'write_logs'.
Input: logtype - string defining the log to be written.
Currently possible:
* charge: charge for full system
* event: electron dynamics events that happened
* pulse: Outputting the pulse profile
* occupation: describing the state
"""
f=open(self.DIR + '/' + logtype + '.log','w')
for t in self.log.get(logtype):
line = ' '.join(str(x) for x in t)
f.write(line + '\n')
f.flush()
f.close()