import numpy as np
from math import sqrt
from scm import plams
from pyorbb.nested_dict import NestedDict
from pyorbb.orbitals import fragments
ensure_list = lambda x: [x] if not isinstance(x, (list, tuple, set)) else list(x) # noqa: E731
def _first_principal_numbers(nfrozen_cores):
aufbau = [
"S",
"S",
"P", "P", "P",
"S",
"P", "P", "P",
"D", "D", "D", "D", "D",
"S",
"P", "P", "P",
"D", "D", "D", "D", "D",
"F", "F", "F", "F", "F", "F", "F",
"S",
"P", "P", "P",
"D", "D", "D", "D", "D",
"S",
"F", "F", "F", "F", "F", "F", "F",
"P", "P", "P",
"D", "D", "D", "D", "D",
"S",
"P", "P", "P",
]
S = 1
P = 2
D = 3
F = 4
for typ in aufbau[:nfrozen_cores]:
if typ == 'S':
S += 1
if typ == 'P':
P += 1/3
if typ == 'D':
D += 1/5
if typ == 'F':
F += 1/7
return {'S': round(S), 'P': round(P), 'D': round(D), 'F': round(F)}
def _get_fragoccupations(reader: plams.KFReader) -> dict:
'''
Read the fragment occupations from a calculation.
Returns:
A dictionary containing the alpha and beta spin occupations for each fragment and irrep.
'''
inp = reader.read('General', 'engine input')
if 'fragoccupations' not in inp.lower():
return {}
lines = []
read = False
for line in inp.splitlines():
line = line.strip()
if line.lower() == 'fragoccupations':
read = True
continue
if read and line.lower() == 'end':
break
if read:
lines.append(line)
lines_lower = [line.lower() for line in lines]
indices = [0]
while 'subend' in lines_lower[indices[-1]+1:]:
index = lines_lower[indices[-1]+1:].index('subend')
indices.append(index + indices[-1] + 2)
blocks = []
for start, end in zip(indices, indices[1:]):
blocks.append(lines[start:end-1])
data = {}
for block in blocks:
frag = block[0]
data[frag] = {}
for occ in block[1:]:
irrep, rest = occ.split(' ', 1)
a, b = rest.split('//')
Na = sum(int(float(part)) for part in a.split())
Nb = sum(int(float(part)) for part in b.split())
data[frag][irrep] = (Na, Nb)
return data
def _get_calc_info(reader: plams.KFReader) -> dict:
'''
Function to read useful info about orbitals from kf reader
'''
ret = NestedDict()
ret.set('engine', 'ADF')
# determine if calculation used relativistic corrections
# if it did, variable 'escale' will be present in 'FMOs'
# if it didnt, only variable 'energy' will be present
ret.set('relativistic', ('SFOs', 'escale') in reader)
ret.set('symlabels', reader.read('Symmetry', 'symlab').strip().split())
ret.set('symmetry', reader.read('Symmetry', 'grouplabel'))
# determine if FMOs are unrestricted or not
ret.set('unrestricted_fmos', ('SFOs', 'energy_B') in reader)
ret.set('fmo_spins', ['A', 'B'] if ret['unrestricted_fmos'] else ['AB'])
# determine if MOs are unrestricted or not
ret.set('unrestricted_mos', (ret['symlabels'][0], 'eps_B') in reader)
ret.set('mo_spins', ['A', 'B'] if ret['unrestricted_mos'] else ['AB'])
frag_data = fragments.get_fragments_data(reader)
# determine if the calculation used regions or not
ret.set('used_regions', not frag_data['used_atomic_fragments'])
ret.set('fragments', frag_data['fragment_names'])
# determine the spin polarization of the complex and fragments
# if we were given fragoccupations we use those
spin_pols = _get_fragoccupations(reader)
for frag in ret['fragments']:
spin_pols.setdefault(frag, {})
# otherwise, if unrestricted, we use the occupations
if ret['unrestricted_fmos']:
frag_index = np.array(reader.read('SFOs', 'fragment'))
subspecies = np.array([subsp.split(':')[0] for subsp in reader.read('SFOs', 'subspecies').split()])
occs_A = np.array(reader.read('SFOs', 'occupation'))
occs_B = np.array(reader.read('SFOs', 'occupation_B'))
for i, frag in enumerate(ret['fragments'], start=1):
subspecies_of_frag = subspecies[frag_index == i]
for subsp in np.unique(subspecies_of_frag):
spin_pols[frag][subsp] = (sum(occs_A[np.logical_and(frag_index == i, subspecies == subsp)]), sum(occs_B[np.logical_and(frag_index == i, subspecies == subsp)]))
ret.set('fmo_spinpolarizations', spin_pols)
# determine if we have access to effective orbital energies
if ('SFOs', 'effective_energy') in reader or 'SFO_Fock_A' in reader or 'SFO_Fock' in reader:
ret.set('has_effective_energy', True)
else:
ret.set('has_effective_energy', False)
return ret
def _square_matrix(S: np.ndarray) -> np.ndarray:
'''
Convert a flattened lower-echelon type matrix into its square matrix.
This is useful for reading data from AMS calculations as they often
store symmetric square matrices in this form, e.g. overlap and Fock matrices.
'''
# lists are easier to work with in this case
S = np.atleast_1d(S).tolist()
size = len(S)
n = int(sqrt(.25 + 2*size) - .5) # number of rows and columns
Srows = [] # this will hold the first m elements of the row
for i in range(n):
# start index will be the number of elements before this row
min_idx = i * (i+1) // 2
# stop index will be the number of elements of the next row
max_idx = (i+1) * (i+2) // 2
Srows.append(S[min_idx:max_idx])
# then we go through rows again and add the remaining (n-m) terms
Srowsfixed = []
for i, row in enumerate(Srows):
Srowsfixed.append(row + [row2[i] for row2 in Srows[i+1:]])
return Srowsfixed
[docs]
def read_data(reader: plams.KFReader, SCF0_reader: plams.KFReader = None, output: str = None):
## HELPER FUNCTIONS
def _read_spin_indep(section, variable, spin, R=reader):
# this function is used to read a section%variable from
# an rkf file that can contain the spin state in its variable name
if spin in ['A', 'AB']:
spin_suffix = '_A'
else:
spin_suffix = '_B'
if not R:
return
# check for section%variable_{spin}
if (section, variable + spin_suffix) in R:
return R.read(section, variable + spin_suffix)
# check for section_{spin}%variable
if (section + spin_suffix, variable) in R:
return R.read(section + spin_suffix, variable)
# check for default case section%variable
if (section, variable) in R:
return R.read(section, variable)
def _compose_vector(data, spins):
# concatenate vectors of all spins
return np.hstack([data[spin] for spin in spins])
def _read_effective_energy(spin, scf0=False):
# function used to read site energies from rkf file with specific spin
# also can read from scf0 kfreader if specified
R = SCF0_reader if scf0 else reader
effective_energy = _read_spin_indep('SFOs', 'effective_energy', fmo_spin, R)
if effective_energy:
return np.atleast_1d(effective_energy)
effective_energy = []
for symlabel in ret['calc_info']['symlabels']:
# This is to correct for symlable being split up into :1, :2, etc., e.g. 1E:1 becomes 1E
if ':' in symlabel and symlabel.split(':')[1].isdigit():
Fock_symlabel = symlabel.split(':')[0]
else:
Fock_symlabel = symlabel
if ret['calc_info']['unrestricted_fmos']:
fock_A = _read_spin_indep('SFO_Fock_A', Fock_symlabel, fmo_spin, R)
fock_B = _read_spin_indep('SFO_Fock_B', Fock_symlabel, fmo_spin, R)
fmats = [fock_A, fock_B]
else:
fock = _read_spin_indep('SFO_Fock', Fock_symlabel, fmo_spin, R)
fmats = [fock]
# the fmat is given as a vector instead of a matrix
# this algorithm is used to read the correct element on the diagonal
# fmat is given as a diagonal matrix and not a full matrix
for fmat in fmats:
if fmat is None:
return
idx = 0
loop = 2
while idx <= len(fmat):
effective_energy.append(fmat[idx])
idx += loop
loop += 1
return np.atleast_1d(effective_energy)
def _read_kinetic_energy():
# the kinetic energy is read from the output file
# and not the rkf file
if not output:
return
with open(output) as outp:
lines = outp.readlines()
read = False
store = []
for line in lines:
if '----------------------------------' in line:
continue
if read and 'Total :' in line:
break
if read:
store.extend(line.strip().split())
if 'Orbital Kinetic Energies (hartree)' in line:
read = True
continue
Ekin = {}
curr_irrep = None
for part in store:
try:
float(part)
is_float = True
except ValueError:
is_float = False
if not is_float and part not in ['(equivalent', 'subspecies)', '----']:
Ekin[part] = []
curr_irrep = part
elif part not in ['(equivalent', 'subspecies)', '----']:
Ekin[curr_irrep].append(float(part))
if part == 'subspecies)':
Ekin[curr_irrep] = Ekin[curr_irrep.split(':')[0] + ':1']
return Ekin
def _compose_matrix(data, spins):
# form a block-diagonal matrix from its corresponding blocks
blocks = [np.array(data[symlabel][spin]) for symlabel in ret['calc_info']['symlabels'] for spin in spins]
shapes = [block.shape for block in blocks]
total_shape = sum(shape[0] for shape in shapes), sum(shape[1] for shape in shapes)
out = np.zeros(total_shape)
current_start_index = 0
for block in blocks:
out[current_start_index:current_start_index + block.shape[0], current_start_index:current_start_index + block.shape[1]] = block
current_start_index += block.shape[0]
return out
## MAIN FUNCTION
ret = NestedDict()
ret.set('calc_info', _get_calc_info(reader))
ret.set('fragment_data', fragments.get_fragments_data(reader))
molecules = ret['fragment_data']['fragment_molecules']
molecules['complex'] = ret['fragment_data']['complex_molecule']
ret.set('molecules', molecules)
# # if we used atomic basis we always have zero spinpol
# for frag in ret['fragment_data']['fragment_names']:
# ret['calc_info']['fmo_spinpolarizations'][frag] = {}
ret.set('FMOs', 'number', reader.read('SFOs', 'number'))
# the name of the fragment
ret.set('FMOs', 'fragment_index', np.atleast_1d(reader.read('SFOs', 'fragment')))
# we use the fragment data to map the fmo fragment index to the fragment name
ret.set('FMOs', 'fragment_types', [ret['fragment_data']['fmo_fragtype_to_fragname_map'][i] for i in ret['FMOs']['fragment_index']])
# the symmlabel of the FMO
ret.set('FMOs', 'subspecies', np.atleast_1d(reader.read('SFOs', 'subspecies').split()))
# index of the FMO in its symmlabel
ret.set('FMOs', 'symmetry_index', np.atleast_1d(reader.read('SFOs', 'isfo')) - 1)
# the index of the FMO in its symlabel
ret.set('FMOs', 'ifo', np.atleast_1d(reader.read('SFOs', 'ifo')) - 1)
ret.set('FMOs', 'spin', [spin for spin in ret['calc_info']['fmo_spins'] for _ in range(ret['FMOs']['number'])])
ret.set('FMOs', 'subspecies_fixed', [])
# some symmetry species can have a subspecies
# for example, C(3V) symmetry has the E1:1 and E1:2 symmetry species
# however, ADF only reports for one of the (general E1 label)
if ret['calc_info']['used_regions']:
subspecies_visited_symm_index = {}
for subsp, ifmo in zip(ret['FMOs']['subspecies'], ret['FMOs']['symmetry_index']):
subspecies_visited_symm_index.setdefault(subsp, [])
if ifmo in subspecies_visited_symm_index[subsp]:
n = int(subsp.split(':')[1])
subsp = subsp.split(':')[0] + ':' + str(n + 1)
subspecies_visited_symm_index.setdefault(subsp, [])
subspecies_visited_symm_index[subsp].append(int(ifmo))
ret['FMOs']['subspecies_fixed'].append(subsp)
else:
# if we have atomic fragments we might have FMOs with the same name for the same fragment
# we should give these unique names as well
# first get the total number of FMOs that have the same subspecies and fragment
total_counts = {}
for subsp, ifo, frag in zip(ret['FMOs']['subspecies'], ret['FMOs']['ifo'], ret['FMOs']['fragment_types']):
total_counts.setdefault(str(frag), {})
total_counts[str(frag)].setdefault(str(subsp), {})
total_counts[str(frag)][str(subsp)].setdefault(int(ifo), 0)
total_counts[str(frag)][str(subsp)][int(ifo)] += 1
counts = {}
# then loop again and set the fixed subspecies
for subsp, ifo, frag in zip(ret['FMOs']['subspecies'], ret['FMOs']['ifo'], ret['FMOs']['fragment_types']):
# if there is only one FMO with this subspecies for this fragment
# we simply set the subspecies as its corrected name
if total_counts[frag][subsp][ifo] == 1:
ret['FMOs']['subspecies_fixed'].append(subsp)
continue
# otherwise we count which FMO we are at and use that
# to name the FMO subspecies
counts.setdefault(str(frag), {})
counts[str(frag)].setdefault(str(subsp), {})
counts[str(frag)][str(subsp)].setdefault(int(ifo), 0)
counts[str(frag)][str(subsp)][int(ifo)] += 1
ret['FMOs']['subspecies_fixed'].append(f'{subsp}/{counts[frag][subsp][ifo]}')
# read basic information about the fmos here
for fmo_spin in ret['calc_info']['fmo_spins']:
ret.set('FMOs', 'energy', fmo_spin, np.atleast_1d(_read_spin_indep('SFOs', 'escale', fmo_spin)))
ret.set('FMOs', 'occupation', fmo_spin, np.atleast_1d(_read_spin_indep('SFOs', 'occupation', fmo_spin)))
# the order in terms of the energy of the FMO
ret.set('FMOs', 'order', fmo_spin, np.argsort(ret['FMOs']['energy'][fmo_spin]))
s = _read_effective_energy(fmo_spin, False)
if s is not None:
ret.set('FMOs', 'effective_energy', fmo_spin, s)
if SCF0_reader:
s = _read_effective_energy(fmo_spin, True)
if s is not None:
ret.set('FMOs', 'effective_energy_SCF0', fmo_spin, s)
for symlabel in ret['calc_info']['symlabels']:
energy_by_symlabel = ret['FMOs']['energy'][fmo_spin][ret['FMOs']['subspecies_fixed'] == symlabel]
ret.set('FMOs', 'order_by_symlabel', symlabel, fmo_spin, np.argsort(energy_by_symlabel))
if 'effective_energy' in ret['FMOs']:
ret.set('FMOs', 'effective_energy', 'total', _compose_vector(ret['FMOs']['effective_energy'], ret['calc_info']['fmo_spins']))
if 'effective_energy_SCF0' in ret['FMOs']:
ret.set('FMOs', 'effective_energy_SCF0', 'total', _compose_vector(ret['FMOs']['effective_energy_SCF0'], ret['calc_info']['fmo_spins']))
ret.set('FMOs', 'energy', 'total', _compose_vector(ret['FMOs']['energy'], ret['calc_info']['fmo_spins']))
ret.set('FMOs', 'occupation', 'total', _compose_vector(ret['FMOs']['occupation'], ret['calc_info']['fmo_spins']))
ret.set('FMOs', 'order', 'total', np.argsort(ret['FMOs']['energy']['total']))
# construct the names of the FMOs as they would appear in ADFLevels
for spin in ret['calc_info']['fmo_spins']:
if spin == 'AB':
ret.set('FMOs', 'adf_names', spin, [f'{index + 1}{symlabel}' for index, symlabel in zip(ret['FMOs']['ifo'], ret['FMOs']['subspecies_fixed'])])
else:
ret.set('FMOs', 'adf_names', spin, [f'{index + 1}{symlabel}_{spin}' for index, symlabel in zip(ret['FMOs']['ifo'], ret['FMOs']['subspecies_fixed'])])
ret.set('FMOs', 'adf_names', 'total', _compose_vector(ret['FMOs']['adf_names'], ret['calc_info']['fmo_spins']))
# construct here the unique names, e.g. ``NH3(4E1:1)`` that contains both the FMO orbital name and its fragment name
# ret.set('FMOs', 'unique_names', {spin: [f'{frag}({name})' for frag, name in zip(ret['FMOs']['fragment_unique'][spin], ret['FMOs']['adf_names'][spin])] for spin in ret['calc_info']['fmo_spins']})
# ret.set('FMOs', 'unique_names', 'total', _compose_vector(ret['FMOs']['unique_names'], ret['calc_info']['fmo_spins']))
# read the matrix data such as overlaps, coefficients, etc.
# we correct for the number of frozen cores later
ret.set('MOs', 'nfrozencores', {symlabel: ncbs for symlabel, ncbs in zip(ret['calc_info']['symlabels'], ensure_list(reader.read('Symmetry', 'ncbs')))})
ret.set('MOs', 'nfrozencores', 'total', sum(ret['MOs']['nfrozencores'].values()))
charge_diff = np.atleast_1d(reader.read('Geometry', 'atomtype total charge')) - np.atleast_1d(reader.read('Geometry', 'atomtype effective charge'))
charge_diff = {atomtype: charge_diff[i] for i, atomtype in enumerate(reader.read('Geometry', 'atomtype').split())}
# also get the frozencores per atom
ret.set('MOs', 'nfrozencores_atom', {atomtype: int(diff/2) for atomtype, diff in charge_diff.items()})
# if a fragment has atomic symmetry we need to fix the principle quantum numbers
# first check if we have atomic fragments, should have S, P, D, or F subspecies names
atomic_fmo_idxs = [i for i, subsp in enumerate(ret['FMOs']['subspecies_fixed']) if subsp.split(':')[0] in 'SPDF']
frag_types = np.atleast_1d(reader.read('SFOs', 'fragtype').split())
# get the fragment names of the atomic fragments
atomic_frags = np.unique([frag_types[i] for i in atomic_fmo_idxs])
# then get the atoms that correspond
frag_atomtype_index = np.atleast_1d(reader.read('Geometry', 'fragment and atomtype index')).reshape(2, -1)
atom_types = {}
for frag_type in atomic_frags:
frag_index = list(np.atleast_1d(reader.read('Geometry', 'fragmenttype').split())).index(frag_type) + 1
atomtype_index = frag_atomtype_index[0].tolist().index(frag_index)
atom_type = list(np.atleast_1d(reader.read('Geometry', 'atomtype').split()))[frag_atomtype_index[1][atomtype_index] - 1]
atom_types[frag_type] = atom_type
first_principal = {frag: _first_principal_numbers(ret['MOs']['nfrozencores_atom'][atom]) for frag, atom in atom_types.items()}
for spin in ret['calc_info']['fmo_spins']:
if spin == 'AB':
ret.set('FMOs', 'adf_names_fixed_principal', spin, [f'{index + first_principal[frag][symlabel.split(":")[0].split("/")[0]]}{symlabel}' if frag in atomic_frags else f'{index + 1}{symlabel}' for index, frag, symlabel in zip(ret['FMOs']['ifo'], frag_types, ret['FMOs']['subspecies_fixed'])])
else:
ret.set('FMOs', 'adf_names_fixed_principal', spin, [f'{index + first_principal[frag][symlabel.split(":")[0].split("/")[0]]}{symlabel}_{spin}' if frag in atomic_frags else f'{index + 1}{symlabel}_{spin}' for index, frag, symlabel in zip(ret['FMOs']['ifo'], frag_types, ret['FMOs']['subspecies_fixed'])])
for symlabel in ret['calc_info']['symlabels']:
for fmo_spin in ret['calc_info']['fmo_spins']:
S = _read_spin_indep(symlabel, 'S-CoreSFO', fmo_spin)
S = _square_matrix(S)
ret.set('matrices', 'overlap', symlabel, fmo_spin, S)
for fmo_spin in ret['calc_info']['fmo_spins']:
F = _read_spin_indep('SFO_Fock', symlabel.split(':')[0], fmo_spin)
if F:
F = _square_matrix(F)
ret.set('matrices', 'fock', symlabel, fmo_spin, F)
for mo_spin in ret['calc_info']['mo_spins']:
nmo = _read_spin_indep(symlabel, 'nmo', mo_spin)
ret.set('MOs', 'number', symlabel, mo_spin, nmo)
ret.set('MOs', 'energy', symlabel, mo_spin, np.atleast_1d(_read_spin_indep(symlabel, 'escale', mo_spin)))
occupation = np.atleast_1d(_read_spin_indep(symlabel, 'froc', mo_spin))
ret.set('MOs', 'occupation', symlabel, mo_spin, occupation)
coefficients = np.atleast_2d(_read_spin_indep(symlabel, 'Eig-CoreSFO', mo_spin))
coefficients = coefficients.reshape(nmo, -1)
ret.set('matrices', 'coefficients', symlabel, mo_spin, coefficients)
# perform mulliken analysis here
# contribution = C * (C @ S)
# population = O * contribution
# see: https://github.com/TheoChem-VU/PyOrbb/issues/28
if ret['calc_info']['fmo_spins'] == ret['calc_info']['mo_spins']:
S = ret['matrices']['overlap'][symlabel][mo_spin]
else:
S = ret['matrices']['overlap'][symlabel]['AB']
contr = coefficients * (coefficients @ S)
contr_pos = abs(contr)
contr_pos = np.maximum(0, contr)
contr_artifact = (np.sum(contr_pos, axis=1, keepdims=True) + np.sum(contr_pos, axis=0, keepdims=True)) / 2
contr_normed = contr_pos / contr_artifact
# contr_normed = (contr.T / np.sum(np.maximum(0, contr), axis=0)).T
# import matplotlib.pyplot as plt
# plt.imshow(contr_artifact)
# plt.show()
# plt.imshow(contr)
# print(contr)
# plt.figure()
# plt.imshow(contr_normed)
# print(contr_normed)
# plt.show()
ret.set('matrices', 'mulliken_contribution', symlabel, mo_spin, contr)
ret.set('matrices', 'mulliken_contribution_normalized', symlabel, mo_spin, contr_normed)
ret.set('matrices', 'mulliken_population', symlabel, mo_spin, np.atleast_2d(occupation).T * contr)
ret.set('MOs', 'energy', 'total', np.hstack([_compose_vector(ret['MOs']['energy'][symlabel], ret['calc_info']['mo_spins']) for symlabel in ret['calc_info']['symlabels']]))
ret.set('MOs', 'occupation', 'total', np.hstack([_compose_vector(ret['MOs']['occupation'][symlabel], ret['calc_info']['mo_spins']) for symlabel in ret['calc_info']['symlabels']]))
ret.set('MOs', 'order', 'total', np.argsort(ret['MOs']['energy']['total']))
ret.set('MOs', 'number', 'total', len(ret['MOs']['energy']['total']))
ret.set('MOs', 'spin', [spin for spin in ret['calc_info']['mo_spins'] for _ in range(ret['MOs']['number']['total'])])
K = _read_kinetic_energy()
if K is not None:
ret.set('MOs', 'kinetic_energy', K)
if 'fock' in ret['matrices']:
ret.set('matrices', 'fock', 'total', _compose_matrix(ret['matrices']['fock'], ret['calc_info']['fmo_spins']))
ret.set('matrices', 'overlap', 'total', _compose_matrix(ret['matrices']['overlap'], ret['calc_info']['fmo_spins']))
ret.set('matrices', 'coefficients', 'total', _compose_matrix(ret['matrices']['coefficients'], ret['calc_info']['mo_spins']))
ret.set('matrices', 'mulliken_contribution', 'total', _compose_matrix(ret['matrices']['mulliken_contribution'], ret['calc_info']['mo_spins']))
ret.set('matrices', 'mulliken_population', 'total', _compose_matrix(ret['matrices']['mulliken_population'], ret['calc_info']['mo_spins']))
# gross population is the vertical marginal of the Mulliken population matrix
ret.set('FMOs', 'symlabel', [])
ret.set('MOs', 'symlabel', [])
ret.set('MOs', 'symmetry_index', [])
for mo_spin in ret['calc_info']['mo_spins']:
ret.set('FMOs', 'gross_population', mo_spin, [])
ret.set('FMOs', 'approx_effective_energy', mo_spin, [])
for symlabel in ret['calc_info']['symlabels']:
norb = ret['MOs']['number'][symlabel][ret['calc_info']['mo_spins'][0]]
ret['FMOs']['symlabel'].extend([symlabel] * norb)
ret['MOs']['symlabel'].extend([symlabel] * norb)
ret['MOs']['symmetry_index'].extend(range(norb))
for mo_spin in ret['calc_info']['mo_spins']:
gp = ret['matrices']['mulliken_population'][symlabel][mo_spin]
gp = np.sum(gp, axis=0)[ret['MOs']['nfrozencores'][symlabel]:]
ret['FMOs']['gross_population'][mo_spin].extend(gp.tolist())
# we also approximate the effective energies here
contr = ret['matrices']['mulliken_contribution'][symlabel][mo_spin]
mo_energy = ret['MOs']['energy'][symlabel][mo_spin]
approx_site = mo_energy @ contr
ret['FMOs']['approx_effective_energy'][mo_spin].extend(approx_site.tolist())
ret.set('FMOs', 'gross_population', 'total', _compose_vector(ret['FMOs']['gross_population'], ret['calc_info']['mo_spins']))
ret.set('FMOs', 'approx_effective_energy', 'total', _compose_vector(ret['FMOs']['approx_effective_energy'], ret['calc_info']['mo_spins']))
return ret
