Skip to content

Add orbital tracking #12

@YHordijk

Description

@YHordijk

Often times, when reading orbital data (e.g. energies) we see that MOs can switch ordering. For example, the HOMO at the start of an IRC becomes the HOMO-1 towards the end. This will lead to incorrect data for these orbitals. We should implement an orbital tracking feature that can fix this issue.
The get_alike_orbital method of the pyfmo.Orbitals class allows us to retrieve orbitals that are similar to a given orbital. We can use this feature to fix the above problem statement.

For example:

from pyfmo import Orbitals
from tcutility import pathfunc
import os
import matplotlib.pyplot as plt

# the directory that stores the pyfrag calculations
dirs = pathfunc.get_subdirectories('../../../test/fixtures/pyfrag')
# get the files to the adf.rkf files
files = [os.path.join(d, 'adf.rkf') for d in sorted(dirs, key=lambda d: int(d.split('.')[-1]))]
# and load them
orbs = [Orbitals(f) for f in files]

# plot energy levels of the MOs based on the ordering in the adf.rkf files
plt.plot([orb.mos['HOMO'].energy for orb in orbs], label='HOMO')
plt.plot([orb.mos['HOMO-1'].energy for orb in orbs], label='HOMO-1')
plt.xlabel('IRC Step')
plt.ylabel(r'$\epsilon$ (eV)')
plt.legend()

def MO_track(start_mo, orbs):
    '''
    Function used to track a MO through a series of calculations.
    
    Args:
        start_mo: the MO object that should be tracked through the other calculations.
        orbs: the other orbitals objects that we use to compare to start_mo

    Returns:
        A list of MO objects that are tracked from start_mo
    '''
    mos = [start_mo]
    for orb in orbs[1:]:
        best = orb.get_alike_orbital(mos[-1])
        mos.append(best)
    return mos

# plot the tracked MO energies
plt.figure()
plt.plot([mo.energy for mo in MO_track(orbs[1].mos['HOMO'], orbs[1:-1])], label='HOMO')
plt.plot([mo.energy for mo in MO_track(orbs[1].mos['HOMO-1'], orbs[1:-1])], label='HOMO-1')
plt.xlabel('IRC Step')
plt.ylabel(r'$\epsilon$ (eV)')
plt.legend()
plt.show()

Figure_1
Figure_2

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workingenhancementNew feature or request

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions