#!/usr/bin/env python
# coding: utf-8
# In[2]:
"""Defines the ORCAWrapper Class"""
# Set up imports
import ase
from ase.calculators.orca import ORCA
from ase.calculators.calculator import PropertyNotImplementedError
from ase.units import Hartree
from os import environ, path, makedirs, remove
import numpy as np
[docs]class ORCAWrapper():
"""Sets up and runs the ORCA Calculator (via ASE) for DFT and TDDFT calculations"""
def __init__(self):
"""Sets up instance attributes for ORCAWrapper """
# Set up some defaults
self.nprocs=24
[docs] def setup(self,nprocs=None,orca_cmd=None,ORCA_top=None):
"""Sets up the internal variables of the ORCAWrapper class, including run command"""
try:
orca_cmd = environ["ASE_ORCA_COMMAND"]
except KeyError:
if orca_cmd is None:
orca_cmd = "~/orca5/orca"
environ["ASE_ORCA_COMMAND"]=f'{orca_cmd} PREFIX.inp --oversubscribe >> PREFIX.out 2> PREFIX.err'
if nprocs is not None:
self.nprocs = nprocs
[docs] def cleanup(self,seed):
"""Cleans up temporary files created by a ORCA run that are of no further use"""
# Cleanup temporary files - TODO
for dir in "./",seed+"/":
for s in ("gbw","densities","err"):
if path.isfile(dir+seed+"."+s):
remove(dir+seed+"."+s)
def unpack_params(self,calc_params):
if 'basis' in calc_params:
basis = calc_params['basis']
else:
raise Exception("Basis not specified in calc_params")
if 'func' in calc_params:
xc = calc_params['func']
else:
raise Exception("XC Func not specified in calc_params")
if 'target' in calc_params:
target = calc_params['target']
else:
raise Exception("Target not specified in calc_params")
if 'disp' in calc_params:
disp = calc_params['disp']
else:
disp = None
return basis,xc,target,disp
def _cosmo_seed(self,solvent):
# See here for a list of solvents:
# http://
if solvent=='meth':
return 'methanol'
elif solvent=='etoh':
return 'ethanol'
elif solvent=='acet':
return 'acetntrl'
elif solvent=='dich':
return 'dcm'
elif solvent=='cycl':
return 'cychexan'
elif solvent=='watr':
return 'water'
elif solvent=='diox':
return 'dioxane'
else:
return solvent
def check_func(self,func):
fullfunc = func
return fullfunc
def _add_tddft(self,calc,nroots,target=None):
calc.parameters['orcablocks'] += f"\n%tddft\n nroots {nroots}\n tda false"
if target is not None:
calc.parameters['orcablocks'] += f"\n Iroot {target}\n end"
[docs] def singlepoint(self,model,label,calc_params={},solvent=None,charge=0,spin=0,
forces=False,continuation=False,readonly=False,calconly=False,
cleanup=True):
"""Runs a singlepoint calculation with the ORCA ASE calculator"""
basis, xc, target, disp = self.unpack_params(calc_params)
dispstr = 'D3ZERO' if disp else ''
extra = ""
calc_orca = ORCA(label=label,orcasimpleinput=f"{xc} {dispstr} {basis} {extra}",
orcablocks=f'%pal nprocs {self.nprocs} end');
if (target is not None) and (target != 0):
self._add_tddft(calc_orca,target,target)
#if solvent is not None:
# self._add_cosmo(calc_orca,solvent)
# Set up spin
calc_orca.set(mult=int(2*spin+1))
#calc_orca.set(dft=self.dftblock_singlepoint)
#if continuation:
# calc_orca.set(restart_kw='restart')
# Set up charge
calc_orca.set(charge=charge)
model.calc = calc_orca
if readonly:
calc_orca.atoms = model
calc_orca.read_results() # skip calculation
if calconly:
return calc_orca
if not forces:
energy = model.get_potential_energy()
if cleanup:
self.cleanup(label)
return energy, calc_orca
else:
forces = model.get_forces()
energy = model.get_potential_energy()
if cleanup:
self.cleanup(label)
return energy, forces, calc_orca
[docs] def geom_opt(self,model_opt,label,calc_params={},driver_tol='default',
solvent=None,continuation=False,charge=0,readonly=False,
calconly=False,cleanup=True):
"""Runs a Geometry Optimisation calculation with the ORCA ASE calculator"""
basis, xc, target, disp = self.unpack_params(calc_params)
dispstr = 'D3ZERO' if disp else ''
calc_orca = ORCA(label=label,orcasimpleinput=f"{xc} {dispstr} {basis} opt",
orcablocks=f'%pal nprocs {self.nprocs} end');
if (target is not None) and (target!=0):
self._add_tddft(calc_orca,target,target)
#if solvent is not None:
# self._add_cosmo(calc_orca,solvent)
#if continuation:
# calc_orca.set(restart_kw='restart')
calc_orca.set(charge=charge)
#calc_orca.set(print='default')
calc_orca.set(task='opt')
#calc_orca.set(memory='2000 mb')
model_opt.calc = calc_orca
if calconly:
return calc_orca
if readonly:
calc_orca.atoms = model_opt
calc_orca.read_results() # skip calculation
forces = model_opt.get_forces()
model_opt.positions = calc_orca.calc.atoms.positions
energy = model_opt.get_potential_energy()
if cleanup:
self.cleanup(label)
return energy, forces, model_opt.positions
[docs] def freq(self,model_opt,label,calc_params={},solvent=None,charge=0,
temp=300,writeonly=False,readonly=False,continuation=False,
cleanup=True):
"""Runs a Vibrational Frequency calculation with the ORCA ASE calculator"""
raise PropertyNotImplementedError("No freq implementation yet")
basis, xc, target, disp = self.unpack_params(calc_params)
calc_orca = ORCA(label=label,basis=basis,xc=self.check_func(xc))
self.freq_block['temp'] = f'1 {temp}'
#calc_orca.set(freq=self.freq_block)
if (target is not None) and (target != 0):
self._add_tddft(calc_orca,target,target)
#if solvent is not None:
# self._add_cosmo(calc_orca,solvent)
#if continuation:
# calc_orca.set(restart_kw='restart')
#if disp:
# self.dftblock_freq['disp'] = {'vdw': '3'}
#calc_orca.set(dft=self.dftblock_freq)
calc_orca.set(charge=charge)
#calc_orca.set(print='default')
#calc_orca.set(memory='2000 mb')
calc_orca.set(task='freq')
#if disp:
# calc_orca.set(disp={'vdw': '3'})
model_opt.calc = calc_orca
if readonly:
calc_orca.atoms = model_opt
print("Reading results")
calc_orca.read_results() # skip calculation
forces = model_opt.get_potential_energy() # Run frequencies
self.read_freq(calc_orca)
if cleanup:
self.cleanup(label)
[docs] def run_qmd(self,model,steplabel,calc_params,qmd_steps,qmd_timestep,superstep,temp,
solvent=None,charge=0,continuation=False,readonly=False,
constraints=None,dynamics=None,cleanup=True):
"""Runs a Quantum Molecular Dynamics calculation with the ORCA ASE calculator"""
raise PropertyNotImplementedError("No qmd implementation yet")
basis, xc, target, disp = self.unpack_params(calc_params)
calc_orca = ORCA(label=steplabel,basis=basis,xc=self.check_func(xc))
qmd_block = {'nstep_nucl': qmd_steps*(superstep+1),
'dt_nucl': qmd_timestep,
'targ_temp': temp,
'com_step': 10,
'thermostat': 'svr 100.0',
'print_xyz': 1}
calc_orca.set(qmd=qmd_block)
if continuation:
calc_orca.set(restart_kw='restart')
if (target is not None) and (target != 0):
self._add_tddft(calc_orca,target,target)
if solvent is not None:
self._add_cosmo(calc_orca,solvent)
if disp:
self.dftblock_qmd['disp'] = {'vdw': '3'}
calc_orca.set(dft=self.dftblock_qmd)
calc_orca.set(print='low')
if constraints is not None:
calc_orca.set(constraints=constraints)
calc_orca.set(charge=charge)
calc_orca.set(memory='2000 mb')
calc_orca.set(task='qmd')
if disp:
calc_orca.set(disp={'vdw': '3'})
model.calc = calc_orca
if readonly:
print("Reading results")
calc_orca.atoms = model
calc_orca.read_results()
# these should not trigger re-runs
energy = model.get_potential_energy()
forces = model.get_forces()
else:
forces = model.get_forces()
energy = model.calc.results["energy"]
model.positions = calc_orca.calc.atoms.positions
model.calc.atoms = calc_orca.calc.atoms
# after call to QMD, atoms /= calc_orca.atoms due to MD motion
#next_model = calc_orca.get_atoms()
next_model = model
if cleanup:
self.cleanup(steplabel)
print(superstep,calc_orca.label,energy,forces[0],model.positions[0])
return energy,forces,model.positions
[docs] def excitations(self,model,label,calc_params={},nroots=1,solvent=None,charge=0,
writeonly=False,readonly=False,continuation=False,cleanup=True,
plot_homo=None,plot_lumo=None,plot_trans_den=None):
"""Calculates TDDFT excitations with the ORCA ASE calculator"""
raise PropertyNotImplementedError("No excitations implementation yet")
# Set up calculator
basis, xc, target, disp = self.unpack_params(calc_params)
calc_orca = ORCA(label=label,basis=basis,xc=self.check_func(xc))
self._add_tddft(calc_orca,nroots,target)
if solvent is not None:
self._add_cosmo(calc_orca,solvent)
calc_orca.set(dft=self.dftblock_excitations)
calc_orca.set(charge=charge)
calc_orca.set(task='energy')
calc_orca.set(memory='2000 mb')
if disp:
calc_orca.set(disp={'vdw': '3'})
if continuation:
calc_orca.set(restart_kw='restart')
model.calc = calc_orca
if writeonly:
calc_orca.write_input(atoms=model)
return 0
if readonly:
calc_orca.read_results()
print("Reading excitations")
s_excit = self.read_excitations(calc_orca)
#model.set_array('singlet_excitations',s_excit)
energy = calc_orca.get_property('energy',atoms=model,allow_calculation=False)
else:
energy = model.get_potential_energy()
print("Reading excitations")
s_excit = self.read_excitations(calc_orca)
if plot_trans_den is not None:
for iexc in range(min(plot_trans_den,nroots)):
self._add_dplot_trans_den(calc_orca,iexc+1)
calc_orca.calculate(model)
if plot_homo is not None:
if isinstance(plothomo,int):
irange = range(plot_homo)
elif isinstance(plothomo,list):
irange = plot_homo
else:
raise Exception(f'Unrecognised type for plot_homo: {type(plot_homo)}')
for ih in irange:
self._add_dplot_orb(calc_orca,ih+1)
calc_orca.calculate(model)
if plot_lumo is not None:
if isinstance(plotlumo,int):
irange = range(plotlumo)
elif isinstance(plotlumo,list):
irange = plot_lumo
else:
raise Exception(f'Unrecognised type for plot_lumo: {type(plot_lumo)}')
for il in range(plot_lumo):
self._add_dplot_orb(calc_orca,il+1)
calc_orca.calculate(model)
if cleanup:
self.cleanup(label)
return s_excit, energy
[docs] def read_excitations(self,calc):
"""Read Excitations from ORCA calculator."""
filename = calc.label+'.nwo'
file = open(filename, 'r')
lines = file.readlines()
file.close()
s_excit = []
trans_lines = []
for i, line in enumerate(lines):
if line.find('ORCA TDDFT Module') >= 0:
# In case of multiple copies of TDDFT output in file,
# return only the last one
s_excit = []
trans_lines = []
for j, line2 in enumerate(lines[i:]):
if line2.find('No. of roots') >= 0:
word = line2.split()
nroots = word[4]
if line2.find('Alpha electrons') >= 0:
word = line2.split()
nalpha = int(word[3])
if line2.find('Root') >= 0 and line2.find('singlet a') >= 0:
word = line2.split()
root = int(word[1])
energy = float(word[6])
word = lines[i+j+5].split()
osc = float(word[3])
tl = []
for k, line3 in enumerate(lines[i+j+10:]):
if '----------' in line3 or len(line3.split())==0:
break
omo = int(line3.split()[1])
umo = int(line3.split()[5])
weight = float(line3.split()[7])
xy = line3.split()[8]
tl.append([nalpha-omo,umo-nalpha-1,weight,xy])
tl = sorted(tl,key = lambda x: x[2]**2,reverse=True)
tl.insert(0, (root,energy,osc))
s_excit.append((root,energy,osc))
trans_lines.append(tl)
if root == nroots:
break
calc.results['excitations'] = np.array(s_excit)
calc.results['transition_origins'] = trans_lines
return s_excit
[docs] def read_freq(self,calc):
"""Read Vibrational Frequencies and Normal Modes from results of ORCA calculator."""
file = open(calc.label+'.nwo', 'r')
lines = file.readlines()
file.close()
nat = len(calc.atoms)
nRa = 3*nat
freq = []
for i, line in enumerate(lines):
if line.find('NORMAL MODE EIGENVECTORS IN CARTESIAN COORDINATES') >= 0:
# In case of multiple copies of frequency output in file,
# return only the last one
freq = np.zeros(nRa)
nmode = np.zeros((nRa,nRa))
maxRa = 0
minRa = 0
Rb = -1
for j, line2 in enumerate(lines[i:i+nRa*nRa/6+nRa/6*10]):
words = line2.split()
if len(words)>0:
if words==[str(e) for e in range(maxRa+1,maxRa+len(words)+1)]:
minRa = maxRa + 1
maxRa = minRa + len(words) - 1
Rb = 0
if words[0] ==str(Rb+1) and len(words)==maxRa-minRa+2:
nmode[Rb,minRa-1:maxRa] = [float(w) for w in words[1:]]
Rb = Rb + 1
if line2.find('Frequency') >= 0:
new_freqs = [float(s) for s in words[1:]]
freq[minRa-1:maxRa]= new_freqs
if line.find('Derivative Dipole Moments') >= 0:
ddip = np.zeros((nRa,3))
Rb = 0
for j, line2 in enumerate(lines[i:i+nRa+10]):
words = line2.split()
if len(words)>0:
if words[0]==str(Rb+1):
ddip[Rb,1:3] = [float(w) for w in words[4:6]]
Rb = Rb + 1
if line.find('Infra Red Intensities') >= 0:
intense = np.zeros((nRa))
Rb = 0
for j, line2 in enumerate(lines[i:i+nRa+10]):
words = line2.split()
if len(words)>0:
if words[0]==str(Rb+1):
intense[Rb] = float(words[3])
Rb = Rb + 1
print(nmode)
print(freq)
print(ddip)
print(intense)
calc.results['frequencies'] = freq
calc.results['normal modes'] = nmode
calc.results['IR intensities'] = intense
calc.results['derivative dipole moments'] = ddip
return freq
# In[ ]: