Source code for CDTK.Dynamics.InputfileSH
#* **************************************************************************
#*
#* 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 CDTK.Tools.Conversion as conv
import numpy as np
from CDTK.Tools.Inputfile import Inputfile
[docs]
class Inputfile_SurfaceHopping(Inputfile):
"""
Read and store the contents of the input file for surface-hopping method.
Returns [sections] - information parsed from input file
"""
def __init__(self,fname=None,**opts):
Inputfile.__init__(self,fname=fname,**opts)
def _inputLog(self,**opts):
"""
Print all parameters according to calculation type [LZ/FS]
"""
fout = 'log'
try:
output = file(fout,'w')
except:
raise ValueError(fout+' could not be accessed')
output.write('------------------------------------------------------------------------\n\n')
output.write('############### NEW RUN ###############\n\n')
output.write('Program name: QDTK \n\n')
output.write('Fewest-Switching Surface-Hopping Method \n\n')
output.write('########################################################################\n\n')
# parsed input parameters for running QDTK
for sec in self.sections:
output.write(str(sec)+' \n')
for sec_content in self.sections[sec]:
if sec_content == 'coord':
output.write(str(sec_content)+'\n')
for atom in self.sections[sec][sec_content]:
for xyz in atom:
output.write("%12.7f" % (xyz,))
output.write('\n')
else:
output.write(str(sec_content)+' = '+str(self.sections[sec][sec_content])+'\n')
output.write('\n')
output.close()
def _process_run_section(self,seclist):
"""
Returns the parameter from the section [run]
Basic definition parameters for the calculation:
Title - calculation description
Method_type - [Newton] newtonian dynamics
[Quantum] quantum dynamics
[Job_type] - [initial_condition] prepare initial condition
[propagation] perform dynamics calculation
[Restart] - whether to be the restart calculation
Molecule -
[Atom] - Array to specify the atoms
[Mass] - Array to specify the associated mass of each atom
Unit - atomic_unit
angstrom
Spatial_dim - Spatial dimensions
Init_condition - File containing initial conditions
TINIT - Initial time
TFINAL - Final time
DT - Time step
"""
dsection = {}
dsection['isrestart'] = False
for line in seclist:
ls_s = self._cleanstr(line).lower() # for options without = seperator
if (ls_s == 'initial_condition') or (ls_s == 'propagation'):
keyword = 'job_type'
value = ls_s
elif (ls_s == 'restart'):
keyword = 'isrestart'
value = True
else:
ls = line.split('=')
keyword = self._cleanstr(ls[0]).lower()
if keyword == 'title':
value = self._cleanstr_long(ls[1])
else:
try:
value = self._cleanstr(ls[1])
except:
value = None
dsection[keyword] = value
# Read the molecular definition
mol_fname = dsection['molecule_file']
try:
fmol = open(mol_fname, 'r')
except:
msg = 'File ' + mol_fname + ' could not be accessed'
raise ValueError(msg)
dsection['atom'] = []
dsection['mass'] = []
if dsection['job_type'] == 'initial_condition':
dsection['coord'] = []
line = fmol.readline()
while line:
if line[0] == '#': # detects comments
line = fmol.readline()
continue
elif line.split() == []: # detects empty lines
line = fmol.readline()
continue
else:
if dsection['job_type'] == 'propagation':
# read atomic name and mass
ls = line.split()
if len(ls) < 2:
raise ValueError('Error: molecule definition file. Usage: [atom] [mass]')
dsection['atom'].append(ls[0])
dsection['mass'].append(float(ls[1]))
line = fmol.readline()
elif dsection['job_type'] == 'initial_condition':
# read atomic name. mass and initial coordinate
ndim = int(dsection['spatial_dim'])
ls = line.split()
if len(ls) != (ndim+2):
raise ValueError('Error: molecule definition file. Usage: [atom] [mass] [coordinates]')
dsection['atom'].append(ls[0])
dsection['mass'].append(float(ls[1]))
for i in range(ndim):
dsection['coord'].append(float(ls[2+i]))
line = fmol.readline()
if dsection['job_type'] == 'initial_condition':
natom = len(dsection['atom'])
dsection['coord'] = np.array(dsection['coord'],dtype=float)
dsection['coord'].shape = (natom, ndim)
fmol.close()
# adjustment for internal keyword and option value consistence
dsection['method_type'] = dsection['method_type'].lower()
dsection['job_type'] = dsection['job_type'].lower()
dsection['surface_hopping_type'] = dsection['surface_hopping_type'].lower()
dsection['unit'] = dsection['unit'].lower()
return dsection
def _process_surface_hopping_section(self,seclist):
"""
Returns the parameter from the section [surface_hopping]
Parameters for the surface-hopping calculation:
init_elec_state - initial electronic state
number_elec_state(s) - number of electronic state(s)
Nstep_quantum - number of quantum steps within one classical step
coupling_matrix_method - [wavefunction_overlap]
[explicit_d_matrix]
"""
dsection = {}
# avoid input ambiguity [number_elec_state(s)]
dsection['number_elec_state'] = None
for line in seclist:
ls = line.split('=')
keyword = self._cleanstr(ls[0]).lower()
try:
value = self._cleanstr(ls[1])
except:
value = None
dsection[keyword] = value
# adjust initial electronic state for internal consistency
dsection['init_elec_state'] = int(dsection['init_elec_state']) - 1
# avoid input ambiguity [number_elec_state(s)]
if dsection['number_elec_state'] is not None:
dsection['number_elec_states'] = dsection['number_elec_state']
dsection.pop('number_elec_state')
# adjustment for internal keyword and option value consistence
dsection['coupling_matrix_method'] = dsection['coupling_matrix_method'].lower()
return dsection
def _process_electronic_struct_section(self,seclist):
"""
Returns the parameter from the section [electronic_struct]
Parameters for the electronic structure calculation:
Program - [MOLCAS]
[QCHEM]
Basis - Basis for electronic wavefunction
Gradient_method - [Analytical]
[Numerical]
Input - General input for the electronic
Energy
Energy method [Koopmans] [CASSCF]
Analytical gradient
Wavefunction
Coupling_matrix_input - Input for the wavefunction overlap calculation
non-adiabatic coupling matrix
"""
dsection = {}
for line in seclist:
ls = line.split('=')
keyword = self._cleanstr(ls[0]).lower()
try:
value = self._cleanstr(ls[1])
except:
value = None
dsection[keyword] = value
# adjustment for internal keyword and option value consistence
program_name = dsection['program']
dsection['program'] = dsection['program'].lower()
dsection['gradient_method'] = dsection['gradient_method'].lower()
dsection['energy_method'] = dsection['energy_method'].lower()
# Prepare the input files for electronic structure calculation
if dsection['program'] == 'molcas':
if dsection['energy_method'] == 'koopmans':
# Hartree-Fock-Koopmans method
pass
else:
# CI type methods
molcas_elec_struct_inp = dsection['energy_input']
molcas_gradient_inp = dsection['gradient_input']
molcas_cplmat_inp_init = dsection['coupling_matrix_input_init']
molcas_cplmat_inp = dsection['coupling_matrix_input']
if not (isinstance(molcas_elec_struct_inp,basestring) and isinstance(molcas_cplmat_inp,basestring)):
raise ValueError('Failed to locate MOLCAS input files')
elif dsection['program'] == 'qchem':
qchem_elec_struct_inp = dsection['input']
qchem_cplmat_inp = dsection['coupling_matrix_input']
if not (isinstance(qchem_elec_struct_inp,basestring) and isinstance(qchem_cplmat_inp,basestring)):
raise ValueError('Failed to locate QCHEM input files')
else:
raise ValueError('Could not find the computer program ' + program_name +' for electronic structure calculation')
return dsection
def _cleanstr(self,s):
"""
remove blanks around 1st word in string
"""
return s.split()[0]
def _cleanstr_long(self,s):
"""
remove blanks around the string
"""
sl = ''
for word in s.split():
sl = sl + word + ' '
return sl
[docs]
def check_input(self,param,keywords):
"""
Check consistency of the input parameters
"""
fout = 'check'
try:
output = file(fout,'w')
except:
raise ValueError(fout+' could not be accessed')
output.write('------------------------------------------------------------------------\n\n')
output.write('############### INPUT CHECK ###############\n\n')
output.write('Program name: QDTK \n\n')
output.write('Fewest-Switching Surface-Hopping Method \n\n')
output.write('########################################################################\n\n')
for section in keywords:
for key in keywords[section]:
try:
# ckeck whether the certain parameter is not defined
if param[section][key] is None:
output.write(section+' '+key+' NOT DEFINED!'+'\n')
else:
output.write(section+' '+key+' checked'+'\n')
except:
# certain parameter is missing
raise ValueError('Input Error: '+section+' '+key+'\n')
output.close()