Skip to content

Commit

Permalink
Merge branch '29-read-vibration-data-from-dftb-calcs' into main
Browse files Browse the repository at this point in the history
YHordijk authored Oct 1, 2023

Unverified

No user is associated with the committer email.
2 parents 378f184 + e39f5fc commit cf8b7b2
Showing 2 changed files with 56 additions and 27 deletions.
45 changes: 32 additions & 13 deletions TCutility/results/adf.py
Original file line number Diff line number Diff line change
@@ -22,14 +22,21 @@ def get_calc_settings(info: Result) -> Result:
'''

assert info.engine == 'adf', f'This function reads ADF data, not {info.engine} data'
assert 'adf.rkf' in info.files, f'Missing adf.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'
# assert 'adf.rkf' in info.files, f'Missing adf.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'

reader_adf = cache.get(info.files['adf.rkf'])
ret = Result()

# set the calculation task at a higher level
ret.task = info.input.task

# the VibrationalAnalysis task does not produce adf.rkf files
# in that case we end the reading here, since the following
# variables do not apply
if ret.task.lower() == 'vibrationalanalysis':
return ret

reader_adf = cache.get(info.files['adf.rkf'])

# determine if calculation used relativistic corrections
# if it did, variable 'escale' will be present in 'SFOs'
# if it didnt, only variable 'energy' will be present
@@ -79,12 +86,32 @@ def get_properties(info: Result) -> Result:
- **vdd.charges (list[float])** - list of Voronoi Deformation Denisty (VDD) charges in [electrons], being the difference between the final (SCF) and initial VDD charges.
'''

def read_vibrations(reader: cache.TrackKFReader) -> Result:
ret = Result()
ret.number_of_modes = reader.read('Vibrations', 'nNormalModes')
ret.frequencies = ensure_list(reader.read('Vibrations', 'Frequencies[cm-1]'))
if ('Vibrations', 'Intensities[km/mol]') in reader:
ret.intensities = ensure_list(reader.read('Vibrations', 'Intensities[km/mol]'))
ret.number_of_imag_modes = len([freq for freq in ret.frequencies if freq < 0])
ret.character = 'minimum' if ret.number_of_imag_modes == 0 else 'transitionstate'
ret.modes = []
for i in range(ret.number_of_modes):
ret.modes.append(reader.read('Vibrations', f'NoWeightNormalMode({i+1})'))
return ret


assert info.engine == 'adf', f'This function reads ADF data, not {info.engine} data'
assert 'adf.rkf' in info.files, f'Missing adf.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'

reader_adf = cache.get(info.files['adf.rkf'])

ret = Result()

if info.adf.task.lower() == 'vibrationalanalysis':
reader_ams = cache.get(info.files['ams.rkf'])
ret.vibrations = read_vibrations(reader_ams)
return ret

reader_adf = cache.get(info.files['adf.rkf'])

# read energies (given in Ha in rkf files)
ret.energy.bond = reader_adf.read('Energy', 'Bond Energy') * constants.HA2KCALMOL
ret.energy.elstat = reader_adf.read('Energy', 'elstat') * constants.HA2KCALMOL
@@ -97,15 +124,7 @@ def get_properties(info: Result) -> Result:

# vibrational information
if ('Vibrations', 'nNormalModes') in reader_adf:
ret.vibrations.number_of_modes = reader_adf.read('Vibrations', 'nNormalModes')
freqs = reader_adf.read('Vibrations', 'Frequencies[cm-1]')
ints = reader_adf.read('Vibrations', 'Intensities[km/mol]')
ret.vibrations.frequencies = freqs if isinstance(freqs, list) else [freqs]
ret.vibrations.intensities = ints if isinstance(ints, list) else [ints]
ret.vibrations.number_of_imag_modes = len([freq for freq in freqs if freq < 0])
ret.vibrations.modes = []
for i in range(ret.vibrations.number_of_modes):
ret.vibrations.modes.append(reader_adf.read('Vibrations', f'NoWeightNormalMode({i+1})'))
ret.vibrations = read_vibrations(reader_adf)

# read the Voronoi Deformation Charges Deformation (VDD) before and after SCF convergence (being "inital" and "SCF")
vdd_scf: List[float] = ensure_list(reader_adf.read('Properties', 'AtomCharge_SCF Voronoi')) # type: ignore since plams does not include typing for KFReader. List[float] is returned
38 changes: 24 additions & 14 deletions TCutility/results/dftb.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
from TCutility.results import cache, Result
from TCutility import ensure_list


def get_calc_settings(info: Result) -> Result:
@@ -9,19 +10,11 @@ def get_calc_settings(info: Result) -> Result:
assert info.engine == 'dftb', f'This function reads DFTB data, not {info.engine} data'
assert 'dftb.rkf' in info.files, f'Missing dftb.rkf file, [{", ".join([": ".join(item) for item in info.files.items()])}]'

reader_dftb = cache.get(info.files['dftb.rkf'])
ret = Result()

# set the calculation task at a higher level
ret.task = info.input.task

# read properties of the calculation
number_of_properties = int(reader_dftb.read('Properties', 'nEntries')) # type: ignore plams does not include type hints. Returns int
for i in range(1, number_of_properties + 1):
prop_type = str(reader_dftb.read('Properties', f'Type({i})')).strip()
prop_subtype = str(reader_dftb.read('Properties', f'Subtype({i})')).strip()
ret.properties[prop_type][prop_subtype] = reader_dftb.read('Properties', f'Value({i})')

return ret


@@ -46,13 +39,30 @@ def get_properties(info: Result) -> Result:
# properties are not given as in ADF. Instead you must first know which index each key is at. For example we might have:
# (Properties, SubType(k)) == Coulomb
# therefore, to get the Coulomb energy we first have to know the index k
# we therefore first loop through every property entry and store the properties in order in a list
# we therefore first loop through every subtype entry and store the substypes in order in a list
# we also store the main type (e.g. 'Energy', 'Gradient', ...) and the value
number_of_properties = reader_dftb.read('Properties', 'nEntries')
properties = []
subtypes = []
types = []
values = []
for i in range(1, number_of_properties+1):
properties.append(reader_dftb.read('Properties', f'Subtype({i})').strip())
print(properties)
subtypes.append(reader_dftb.read('Properties', f'Subtype({i})').strip())
types.append(reader_dftb.read('Properties', f'Type({i})').strip())
values.append(reader_dftb.read('Properties', f'Value({i})'))

# then simply add the properties to ret
for typ, subtyp, value in zip(types, subtypes, values):
ret[typ.replace(' ', '_')][subtyp] = value

# we also read vibrations
ret.vibrations.number_of_modes = reader_dftb.read('Vibrations', 'nNormalModes')
ret.vibrations.frequencies = ensure_list(reader_dftb.read('Vibrations', 'Frequencies[cm-1]'))
if ('Vibrations', 'Intensities[km/mol]') in reader_dftb:
ret.vibrations.intensities = ensure_list(reader_dftb.read('Vibrations', 'Intensities[km/mol]'))
ret.vibrations.number_of_imag_modes = len([freq for freq in ret.vibrations.frequencies if freq < 0])
ret.vibrations.character = 'minimum' if ret.vibrations.number_of_imag_modes == 0 else 'transitionstate'
ret.vibrations.modes = []
for i in range(ret.vibrations.number_of_modes):
ret.vibrations.modes.append(reader_dftb.read('Vibrations', f'NoWeightNormalMode({i+1})'))

# read energies (given in Ha i
return ret

0 comments on commit cf8b7b2

Please sign in to comment.