#* **************************************************************************
#*
#* 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 math as mt
import numpy as np
HUGE = 1.0e15-1.0
[docs]
class SurfaceHoppingLZController(object):
"""
Controller for Landau-Zener surface-hopping
"""
def __init__(self,**opts):
self.calls = 0
self.state = opts.get('state',0)
self.potUP = np.zeros(4,float) # array to keep last upper pot values
self.potDOWN = np.zeros(4,float) # array to keep last lower pot values
self.pot = np.zeros(4,float) # array to keep last current pot values
self.potUP = HUGE
self.potDOWN = -HUGE
[docs]
def update(self,v,m,dt,pot_cur,pot_down=-HUGE,pot_up=HUGE):
"""
Update state of controller and decide if a hopp takes place
v -- array with velocity of each coordinate
m -- array with mass for each coordinate
dt -- time-step between updates
pot_cur -- Current PES
pot_down-- Lower PES w.r.t. populated PES
pot_up -- Upper PES w.r.t. populated PES
Returns [hopp,v]
hopp -- -1,0,1 (hopp down, no hopping, hopp up)
v -- new velocities (old ones if hopp=0)
"""
self.calls += 1
self.update_potentials([pot_down,pot_cur,pot_up])
# Check hopp downwards
if self.is_gap_min_down():
egap,egapdt2 = gap_parameters(self.pot,self.potDOWN,dt)
if is_hopp_energy(v,self.pot[0],self.potDOWN[0],m):
p = hopping_probability(egap,egapdt2)
if is_hopp(p):
vnew = velocity_correction(v,self.pot[0],self.potDOWN[0],m)
self.potUP = self.pot.copy()
self.pot = self.potDOWN.copy()
self.potDOWN[:] = -HUGE
return -1,vnew
# Check hopp upwards
elif self.is_gap_min_up():
egap,egapdt2 = gap_parameters(self.pot,self.potUP,dt)
if is_hopp_energy(v,self.pot[0],self.potUP[0],m):
p = hopping_probability(egap,egapdt2)
if is_hopp(p):
vnew = velocity_correction(v,self.pot[0],self.potUP[0],m)
self.potDOWN = self.pot.copy()
self.pot = self.potUP.copy()
self.potUP[:] = HUGE
return 1,vnew
return 0,v
[docs]
def update_potentials(self,pots):
self.potDOWN = np.roll(self.potDOWN,1)
self.pot = np.roll(self.pot,1)
self.potUP = np.roll(self.potUP,1)
self.potDOWN[0] = pots[0]
self.pot[0] = pots[1]
self.potUP[0] = pots[2]
[docs]
def is_gap_min_down(self):
g = self.pot - self.potDOWN
if g[1] < g[0] and g[1] < g[2]: # min found
return True
else:
return False
[docs]
def is_gap_min_up(self):
g = self.potUP - self.pot
if g[1] < g[0] and g[1] < g[2]: # min found
return True
else:
return False
[docs]
def hopping_probability(egap,egapdt2):
"""
Return the Landau-Zener surface hopping probability (0,1)
egap -- energy gap between PES (absolute value)
egapdt2 -- second derivative of egap w.r.t. time
"""
return mt.exp( -0.5*mt.pi*mt.sqrt(egap**3/egapdt2) )
[docs]
def velocity_correction(ivel,ipes,fpes,m):
"""
Return new velocities after hopping
ivel -- array with initial velocities
ipes -- initial PES
fpes -- final PES
m -- array with mass of each coordinate
Simple velocity rescaling. Other conditions such as angular momentum conservation may
be included afterwards.
"""
ike = 0.5*np.dot(m*ivel,ivel) # initial kinetic energy
dpes = ipes - fpes # initial minus final potential energy
k2 = 1.0 + dpes/ike # square of correction factor
try:
k = mt.sqrt(k2)
except:
ValueError('not enough KE for surface hopping')
return k*ivel
[docs]
def is_hopp_energy(ivel,ipes,fpes,m):
"""
Return True if enough kinetic energy is available for a hopping event, else False
ivel -- array with initial velocities
ipes -- initial PES
fpes -- final PES
m -- array with mass of each coordinate
"""
ike = 0.5*np.dot(m*ivel,ivel) # initial kinetic energy
dpes = ipes - fpes # initial minus final potential energy
k2 = 1.0 + dpes/ike # square of correction factor
if k2 > 0.0:
return True
else:
return False
[docs]
def is_hopp(plz):
"""
Return True if hopping event has to occur, False otherwise
"""
r = np.random.random()
if r <= plz:
return True
else:
return False
[docs]
def gap_parameters(pot1,pot2,dt):
"""
Determine gap and gap second time derivative
pot1 -- latest four potential values on surface 1
pot2 -- latest four potential values on surface 2
dt -- time-step
"""
g = np.abs( pot1 - pot2 )
d2g = (g[0] - 2.0*g[1] + g[2])/dt**2
return g[1],d2g