29 read vibration data from dftb calcs

TCutility/__init__.py

@@ -1,3 +1,4 @@
 ensure_list = lambda x: [x] if not isinstance(x, (list, tuple, set)) else x  # noqa: E731
+squeeze_list = lambda x: x[0] if len(x) == 1 else x  # noqa: E731
 
 from TCutility import results, constants  # noqa: F401, E402

TCutility/results/__init__.py

@@ -29,6 +29,7 @@ def read(calc_dir: str) -> Result:
     if ret.engine == 'adf':
         ret.adf = adf.get_calc_settings(ret)
         ret.properties = adf.get_properties(ret)
+        ret.level = adf.get_level_of_theory(ret)
     elif ret.engine == 'dftb':
         ret.dftb = dftb.get_calc_settings(ret)
         ret.properties = dftb.get_properties(ret)

TCutility/results/adf.py

@@ -22,9 +22,8 @@ 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'])
     reader_ams = cache.get(info.files['ams.rkf'])
     ret = Result()
 
@@ -41,6 +40,14 @@ def get_calc_settings(info: Result) -> Result:
             ret.task = words[i+1]
             break
 
+    # 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
@@ -90,12 +97,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
@@ -108,15 +135,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
@@ -129,3 +148,79 @@ def get_properties(info: Result) -> Result:
     # The ordering is not very straightfoward so this is a suggestion for the future with keys: ret.vdd.[IRREP]
 
     return ret
+
+
+def get_level_of_theory(info: Result) -> Result:
+    '''Function to get the level-of-theory from an input-file.
+
+    Args:
+        inp_path: Path to the input file. Can be a .run or .in file create for AMS
+
+    Returns:
+        :Dictionary containing information about the level-of-theory:
+
+            - **summary (str)** - a summary string that gives the level-of-theory in a human-readable format.
+            - **xc.functional (str)** - XC-functional used in the calculation.
+            - **xc.category (str)** - category the XC-functional belongs to. E.g. GGA, MetaHybrid, etc ...
+            - **xc.dispersion (str)** - the dispersion correction method used during the calculation.
+            - **xc.summary (str)** - a summary string that gives the XC-functional information in a human-readable format.
+            - **xc.empiricalScaling (str)** - which empirical scaling parameter was used. Useful for MP2 calculations.
+            - **basis.type (str)** - the size/type of the basis-set.
+            - **basis.core (str)** - the size of the frozen-core approximation.
+            - **quality (str)** - the numerical quality setting.
+    '''
+    sett = info.input
+    ret = Result()
+    # print(json.dumps(sett, indent=4))
+    xc_categories = ['GGA', 'LDA', 'MetaGGA', 'MetaHybrid', 'Model', 'LibXC', 'DoubleHybrid', 'Hybrid', 'MP2', 'HartreeFock']
+    ret.xc.functional = 'VWN'
+    ret.xc.category = 'LDA'
+    for cat in xc_categories:
+        if cat.lower() in [key.lower() for key in sett.adf.xc]:
+            ret.xc.functional = sett.adf.xc[cat]
+            ret.xc.category = cat
+
+    ret.basis.type = sett.adf.basis.type
+    ret.basis.core = sett.adf.basis.core
+    ret.quality = sett.adf.NumericalQuality or 'Normal'
+
+    ret.xc.dispersion = None
+    if 'dispersion' in [key.lower() for key in sett.adf.xc]:
+        ret.xc.dispersion = " ".join(sett.adf.xc.dispersion.split())
+
+    # the empirical scaling value is used for MP2 calculations
+    ret.xc.empirical_scaling = None
+    if 'empiricalscaling' in [key.lower() for key in sett.adf.xc]:
+        ret.xc.empiricalscaling = sett.adf.xc.empiricalscaling
+
+    # MP2 and HF are a little bit different from the usual xc categories. They are not key-value pairs but only the key
+    # we start building the ret.xc.summary string here already. This will contain the human-readable functional name
+    if ret.xc.category == 'MP2':
+        ret.xc.summary = 'MP2'
+        if ret.xc.empiricalscaling:
+            ret.xc.summary += f'-{ret.xc.empiricalscaling}'
+    elif ret.xc.category == 'HartreeFock':
+        ret.xc.summary = 'HF'
+    else:
+        ret.xc.summary = ret.xc.functional
+
+    # If dispersion was used, we want to add it to the ret.xc.summary
+    if ret.xc.dispersion:
+        if ret.xc.dispersion.lower() == 'grimme3':
+            ret.xc.summary += '-D3'
+        if ret.xc.dispersion.lower() == 'grimme3 bjdamp':
+            ret.xc.summary += '-D3(BJ)'
+        if ret.xc.dispersion.lower() == 'grimme4':
+            ret.xc.summary += '-D4'
+        if ret.xc.dispersion.lower() == 'ddsc':
+            ret.xc.summary += '-dDsC'
+        if ret.xc.dispersion.lower() == 'uff':
+            ret.xc.summary += '-dUFF'
+        if ret.xc.dispersion.lower() == 'mbd':
+            ret.xc.summary += '-MBD@rsSC'
+        if ret.xc.dispersion.lower() == 'default':
+            ret.xc.summary += '-D'
+
+    # ret.summary is simply the ret.xc.summary plus the basis set type
+    ret.summary = f'{ret.xc.summary}/{ret.basis.type}'
+    return ret

TCutility/results/ams.py

@@ -4,6 +4,7 @@
 from TCutility import ensure_list
 import os
 from datetime import datetime
+import re
 
 j = os.path.join
 
@@ -124,6 +125,9 @@ def get_ams_info(calc_dir: str) -> Result:
     else:
         ret.engine = 'adf'
 
+    # store the input of the calculation
+    ret.input = get_ams_input(reader_ams.read('General', 'user input'))
+
     # store the job id, which should be unique for the job
     ret.job_id = reader_ams.read('General', 'jobid') if ('General', 'jobid') in reader_ams else None
 
@@ -375,3 +379,152 @@ def get_history(calc_dir: str) -> Result:
                     ret[item.lower()].append(val)
 
     return ret
+
+
+def get_input_blocks():
+    '''
+    This function reads input_blocks and decomposes its content into a list of blocks and a list of non-standard blocks
+    The general format is as follows:
+
+    parentblock
+    - subblock
+    - - subsubblock
+    - subblock !nonstandard
+    parentblock
+    - subblock
+    - - subsubblock 
+    - - - subsubsubblock
+    - - - subsubsubblock !nonstandard
+
+    Each subblock has to be defined within its parent block. !nonstandard indicates that the block is a non-standard block
+    These blocks are special in that they can contain multiple of the same entry
+    '''
+    blocks = []
+    nonstandard_blocks = []
+    parent_blocks = []  # this list tracks the parent blocks of the current block
+    with open(j(os.path.split(__file__)[0], 'input_blocks')) as inpblx:
+        lines = inpblx.readlines()
+
+    for line in lines:
+        line = line.strip().lower()
+        # we can ignore some lines
+        if line == '' or line.startswith('#'):
+            continue
+
+        # block_depth indicates how many parents the block has
+        block_depth = line.count('- ')
+        # we reduce the amount of parent_blocks using block_depth
+        # if we move from a subsubblock to a subblock we remove the last-added block
+        parent_blocks = parent_blocks[:block_depth]
+
+        # remove the "- " substrings
+        block = line.split('- ')[-1].strip().lower()
+        # check if the block is non-standard. If it is, remove the !nonstandard substring
+        # and add it to the nonstandard_blocks list
+        if block.endswith('!nonstandard'):
+            block = block.split()[0]
+            nonstandard_blocks.append(parent_blocks.copy() + [block])
+        # in both standard and nonstandard cases add the block to blocks list
+        blocks.append(parent_blocks.copy() + [block])
+        # add the current block to parent_blocks for the next line
+        parent_blocks.append(block.lower())
+
+    return blocks, nonstandard_blocks
+
+
+def get_ams_input(inp: str) -> Result:
+    def get_possible_blocks():
+        # check which blocks are openable given the current open blocks
+        # we start by considering all blocks
+        possible_blocks = blocks
+        # we iterate through the open blocks and check if the first element in the possible_blocks is the same.
+        # this way we move down through the blocks
+        for open_block in open_blocks:
+            possible_blocks = [block[1:] for block in possible_blocks if len(block) > 0 and block[0].lower() == open_block.lower()]
+        # we want only the first element in the possible blocks, not the tails
+        possible_blocks = set([block[0] for block in possible_blocks if len(block) > 0])
+        return possible_blocks
+
+    sett = Result()
+
+    blocks, nonstandard_blocks = get_input_blocks()
+
+    open_blocks = ['ams']
+    for line in inp.splitlines():
+        line = line.strip()
+
+        # we remove comments from the line
+        # comments can be either at the start of a line or after a real statement
+        # so we have to search the line for comment starters and remove the part after it
+        for comment_start in ['#', '!', '::']:
+            if comment_start in line:
+                idx = line.index(comment_start)
+                line = line[:idx]
+
+        # skip empty lines
+        if not line:
+            continue
+
+        # if we encounter an end statement we can close the last openblock
+        if line.lower().startswith('end'):
+            open_blocks.pop(-1)
+            continue
+
+        # check if we are opening a block
+        # We have to store the parts of a compact block (if it is compact) 
+        # it will be a list of tuples of key-value pairs
+        compact_parts = None
+        skip = False
+        # check if the current line corresponds to a possible block
+        for possible_block in get_possible_blocks():
+            if line.lower().startswith(possible_block.lower()):
+                # a block opening can either span multiple lines or can be use with compact notation
+                # compact notation will always have = signs
+                if '=' in line:
+                    # split the key-values using some regex magic
+                    compact = line[len(possible_block):].strip()  # remove the block name
+                    compact_parts = re.findall(r"""(\S+)=['"]{0,1}([^'"\n]*)['"]{0,1}""", compact)
+                else:
+                    skip = True
+                # add the new block to the open_blocks
+                open_blocks.append(possible_block)
+
+        # if we are not in a compact block we can just skip
+        if skip:
+            continue
+
+        # get the sett to the correct level
+        # first check if the block is nonstandard
+        is_nonstandard = [block.lower() for block in open_blocks] in nonstandard_blocks
+        sett_ = sett
+        # go to the layer one above the lowest
+        for open_block in open_blocks[:-1]:
+            sett_ = sett_[open_block]
+        # at the lowest level we have to check if the block is nonstandard. If it is, we add the whole line to the settings object
+        if is_nonstandard and not sett_[open_blocks[-1]]:
+            sett_[open_blocks[-1]] = []
+        # then finally go to the lowest sett layer
+        sett_ = sett_[open_blocks[-1]]
+
+        # if we are in a compact block we just add all the key-value pairs
+        if compact_parts:
+            for key, val in compact_parts:
+                sett_[key] = val
+            # compact blocks do not have end statements, so we can remove it again from the open_blocks
+            open_blocks.pop(-1)
+        # in a normal block we split the line and add them to the settings
+        else:
+            # for nonstandard blocks we add the line to the list
+            if is_nonstandard:
+                sett_.append(line.strip())    
+            # for normal blocks we split the line and set it as a key, the rest of the line is the value
+            else:
+                sett_[line.split()[0]] = line[len(line.split()[0]):].strip()   
+
+    for engine_block in ['engine adf', 'engine dftb', 'engine band']:
+        if engine_block not in sett['ams']:
+            continue
+        sett['ams'][engine_block.split()[1]] = sett['ams'][engine_block]
+        del sett['ams'][engine_block]
+
+    return sett['ams']