Skip to content

Commit

Permalink
Merge pull request theochem#99 from tovrstra/moprops-namedtuple
Browse files Browse the repository at this point in the history 
Add properties to MolecularOrbitals class
  • Loading branch information
@tovrstra
tovrstra authored May 13, 2019
2 parents e71e34e + b2766ac commit ddeaa47
Show file tree
Hide file tree
Showing 14 changed files with 447 additions and 326 deletions.
25 changes: 12 additions & 13 deletions iodata/formats/cp2k.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@
from scipy.special import factorialk

from ..basis import angmom_sti, MolecularBasis, Shell, HORTON2_CONVENTIONS
from ..utils import LineIterator, MolecularOrbitals
from ..orbitals import MolecularOrbitals
from ..utils import LineIterator


__all__ = []
Expand Down Expand Up @@ -447,20 +448,17 @@ def load_one(lit: LineIterator) -> dict:

# Turn orbital data into a MolecularOrbitals object.
if restricted:
mo_type = 'restricted'
norb, nel = _get_norb_nel(oe_alpha)
norb_alpha = norb_beta = norb
assert nel % 2 == 0
orb_alpha_coeffs = np.zeros([obasis.nbasis, norb])
orb_alpha_energies = np.zeros(norb)
orb_alpha_occs = np.zeros(norb)
_fill_orbitals(orb_alpha_coeffs, orb_alpha_energies, orb_alpha_occs,
oe_alpha, coeffs_alpha, obasis, restricted)
mo_occs = orb_alpha_occs
mo_coeffs = orb_alpha_coeffs
mo_energy = orb_alpha_energies
mo = MolecularOrbitals(
'restricted', norb, norb, 2 * orb_alpha_occs, orb_alpha_coeffs,
orb_alpha_energies, None)
else:
mo_type = 'unrestricted'
norb_alpha = _get_norb_nel(oe_alpha)[0]
norb_beta = _get_norb_nel(oe_beta)[0]
assert norb_alpha == norb_beta
Expand All @@ -475,12 +473,13 @@ def load_one(lit: LineIterator) -> dict:
_fill_orbitals(orb_beta_coeffs, orb_beta_energies, orb_beta_occs,
oe_beta, coeffs_beta, obasis, restricted)

mo_occs = np.concatenate((orb_alpha_occs, orb_beta_occs), axis=0)
mo_energy = np.concatenate((orb_alpha_energies, orb_beta_energies), axis=0)
mo_coeffs = np.concatenate((orb_alpha_coeffs, orb_beta_coeffs), axis=1)

# create a MO namedtuple
mo = MolecularOrbitals(mo_type, norb_alpha, norb_beta, mo_occs, mo_coeffs, None, mo_energy)
mo = MolecularOrbitals(
'unrestricted', norb_alpha, norb_beta,
np.concatenate((orb_alpha_occs, orb_beta_occs), axis=0),
np.concatenate((orb_alpha_coeffs, orb_beta_coeffs), axis=1),
np.concatenate((orb_alpha_energies, orb_beta_energies), axis=0),
None,
)

result = {
'obasis': obasis,
Expand Down
16 changes: 7 additions & 9 deletions iodata/formats/fchk.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,8 @@
import numpy as np

from ..basis import MolecularBasis, Shell, HORTON2_CONVENTIONS
from ..utils import LineIterator, MolecularOrbitals, amu
from ..orbitals import MolecularOrbitals
from ..utils import LineIterator, amu


__all__ = []
Expand Down Expand Up @@ -190,31 +191,28 @@ def load_one(lit: LineIterator) -> dict:

norba = fchk['Alpha Orbital Energies'].shape[0]
mo_coeffs = np.copy(fchk['Alpha MO coefficients'].reshape(nbasis_indep, nbasis).T)
mo_energy = np.copy(fchk['Alpha Orbital Energies'])
mo_energies = np.copy(fchk['Alpha Orbital Energies'])

if 'Beta Orbital Energies' in fchk:
# unrestricted
mo_type = 'unrestricted'
norbb = fchk['Beta Orbital Energies'].shape[0]
mo_coeffs_b = np.copy(fchk['Beta MO coefficients'].reshape(nbasis_indep, nbasis).T)
mo_coeffs = np.concatenate((mo_coeffs, mo_coeffs_b), axis=1)
mo_energy = np.concatenate((mo_energy, np.copy(fchk['Beta Orbital Energies'])), axis=0)
mo_energies = np.concatenate((mo_energies, np.copy(fchk['Beta Orbital Energies'])), axis=0)
mo_occs = np.zeros(2 * nbasis_indep)
mo_occs[:nalpha] = 1.0
mo_occs[nbasis_indep: nbasis_indep + nbeta] = 1.0
mo = MolecularOrbitals('unrestricted', norba, norbb, mo_occs, mo_coeffs, mo_energies, None)
else:
# restricted closed-shell and open-shell
mo_type = 'restricted'
norbb = norba
mo_occs = np.zeros(nbasis_indep)
mo_occs[:nalpha] = 1.0
mo_occs[:nbeta] = 2.0
if nalpha != nbeta:
# delete dm_full_scf because it is known to be buggy
result['one_rdms'].pop('scf')

# create a MO namedtuple
result['mo'] = MolecularOrbitals(mo_type, norba, norbb, mo_occs, mo_coeffs, None, mo_energy)
mo = MolecularOrbitals('restricted', norba, norba, mo_occs, mo_coeffs, mo_energies, None)
result['mo'] = mo

# E) Load properties
if 'Polarizability' in fchk:
Expand Down
164 changes: 86 additions & 78 deletions iodata/formats/molden.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,9 @@
convert_conventions, HORTON2_CONVENTIONS)
from ..iodata import IOData
from ..periodic import sym2num, num2sym
from ..orbitals import MolecularOrbitals
from ..overlap import compute_overlap, gob_cart_normalization
from ..utils import angstrom, LineIterator, MolecularOrbitals
from ..utils import angstrom, LineIterator


__all__ = []
Expand Down Expand Up @@ -103,12 +104,12 @@ def _load_low(lit: LineIterator) -> dict:
atnums = None
atcoords = None
obasis = None
coeff_alpha = None
ener_alpha = None
occ_alpha = None
coeff_beta = None
ener_beta = None
occ_beta = None
coeffsa = None
energiesa = None
occsa = None
coeffsb = None
energiesb = None
occsb = None
title = None

line = next(lit)
Expand Down Expand Up @@ -153,8 +154,8 @@ def _load_low(lit: LineIterator) -> dict:
# molecular-orbital coefficients.
elif line == '[mo]':
data_alpha, data_beta = _load_helper_coeffs(lit)
coeff_alpha, ener_alpha, occ_alpha = data_alpha
coeff_beta, ener_beta, occ_beta = data_beta
occsa, coeffsa, energiesa, irrepsa = data_alpha
occsb, coeffsb, energiesb, irrepsb = data_beta

# Assign pure and Cartesian correctly. This needs to be done after reading
# because the tags for pure functions may come after the basis set.
Expand All @@ -163,25 +164,21 @@ def _load_low(lit: LineIterator) -> dict:
if shell.angmoms[0] in pure_angmoms:
shell.kinds[0] = 'p'

if coeff_beta is None:
mo_type = 'restricted'
if coeff_alpha.shape[0] != obasis.nbasis:
if coeffsb is None:
if coeffsa.shape[0] != obasis.nbasis:
lit.error("Number of alpha orbital coefficients does not match the size of the basis.")
norba = norbb = coeff_alpha.shape[1]
mo_occs = occ_alpha
mo_coeffs = coeff_alpha
mo_energy = ener_alpha
mo = MolecularOrbitals(
'restricted', coeffsa.shape[1], coeffsa.shape[1],
occsa, coeffsa, energiesa, irrepsa)
else:
mo_type = 'unrestricted'
if coeff_beta.shape[0] != obasis.nbasis:
if coeffsb.shape[0] != obasis.nbasis:
lit.error("Number of beta orbital coefficients does not match the size of the basis.")
norba = coeff_alpha.shape[1]
norbb = coeff_beta.shape[1]
mo_occs = np.concatenate((occ_alpha, occ_beta), axis=0)
mo_energy = np.concatenate((ener_alpha, ener_beta), axis=0)
mo_coeffs = np.concatenate((coeff_alpha, coeff_beta), axis=1)
# create a MO namedtuple
mo = MolecularOrbitals(mo_type, norba, norbb, mo_occs, mo_coeffs, None, mo_energy)
mo = MolecularOrbitals(
'unrestricted', coeffsa.shape[1], coeffsb.shape[1],
np.concatenate((occsa, occsb), axis=0),
np.concatenate((coeffsa, coeffsb), axis=1),
np.concatenate((energiesa, energiesb), axis=0),
irrepsa + irrepsb)

result = {
'atcoords': atcoords,
Expand Down Expand Up @@ -252,16 +249,18 @@ def _load_helper_obasis(lit: LineIterator) -> MolecularBasis:

def _load_helper_coeffs(lit: LineIterator) -> Tuple:
"""Load the orbital coefficients."""
coeff_alpha = []
ener_alpha = []
occ_alpha = []
coeff_beta = []
ener_beta = []
occ_beta = []
occsa = []
coeffsa = []
energiesa = []
irrepsa = []
occsb = []
coeffsb = []
energiesb = []
irrepsb = []

while True:
try:
line = next(lit).lower().strip()
line = next(lit).strip()
except StopIteration:
# We have no proper way to check when a Molden file has ended, so
# we must anticipate for it here.
Expand All @@ -284,18 +283,21 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
break
key, value = line.split('=')
info[key.strip().lower()] = value
energy = float(info['ene'])
occ = float(info['occup'])
col = []
energy = float(info['ene'])
irrep = info.get('sym', '??').strip()
# store column of coefficients, i.e. one orbital, energy and occ
if info['spin'].strip().lower() == 'alpha':
coeff_alpha.append(col)
ener_alpha.append(energy)
occ_alpha.append(occ)
occsa.append(occ)
coeffsa.append(col)
energiesa.append(energy)
irrepsa.append(irrep)
else:
coeff_beta.append(col)
ener_beta.append(energy)
occ_beta.append(occ)
occsb.append(occ)
coeffsb.append(col)
energiesb.append(energy)
irrepsb.append(irrep)
for line in lit:
words = line.split()
if len(words) != 2 or not words[0].isdigit():
Expand All @@ -305,18 +307,18 @@ def _load_helper_coeffs(lit: LineIterator) -> Tuple:
break
col.append(float(words[1]))

coeff_alpha = np.array(coeff_alpha).T
ener_alpha = np.array(ener_alpha)
occ_alpha = np.array(occ_alpha)
if not coeff_beta:
coeff_beta = None
ener_beta = None
occ_beta = None
coeffsa = np.array(coeffsa).T
energiesa = np.array(energiesa)
occsa = np.array(occsa)
if not coeffsb:
coeffsb = None
energiesb = None
occsb = None
else:
coeff_beta = np.array(coeff_beta).T
ener_beta = np.array(ener_beta)
occ_beta = np.array(occ_beta)
return (coeff_alpha, ener_alpha, occ_alpha), (coeff_beta, ener_beta, occ_beta)
coeffsb = np.array(coeffsb).T
energiesb = np.array(energiesb)
occsb = np.array(occsb)
return (occsa, coeffsa, energiesa, irrepsa), (occsb, coeffsb, energiesb, irrepsb)


def _is_normalized_properly(obasis: MolecularBasis, atcoords: np.ndarray,
Expand Down Expand Up @@ -522,23 +524,25 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
"""
obasis = result['obasis']
atcoords = result['atcoords']
if result['mo'].type == 'restricted':
coeffs_a = result['mo'].coeffs
coeffs_b = None
elif result['mo'].type == 'unrestricted':
coeffs_a = result['mo'].coeffs[:, :result['mo'].norba]
coeffs_b = result['mo'].coeffs[:, result['mo'].norba:]
if result['mo'].kind == 'restricted':
coeffsa = result['mo'].coeffs
# Skip testing coeffsb if it is the same array as coeffsa.
coeffsb = None
elif result['mo'].kind == 'unrestricted':
coeffsa = result['mo'].coeffsa
coeffsb = result['mo'].coeffsb
else:
raise ValueError('Molecular orbital type={0} not recognized'.format(result['mo'].type))
if _is_normalized_properly(obasis, atcoords, coeffs_a, coeffs_b):
raise ValueError('Molecular orbital kind={0} not recognized'.format(result['mo'].kind))

if _is_normalized_properly(obasis, atcoords, coeffsa, coeffsb):
# The file is good. No need to change obasis.
return
print('5:Detected incorrect normalization of orbitals loaded from a Molden or MKL file.')

# --- ORCA
print('5:Trying to fix it as if it was a file generated by ORCA.')
orca_obasis = _fix_obasis_orca(obasis)
if _is_normalized_properly(orca_obasis, atcoords, coeffs_a, coeffs_b):
if _is_normalized_properly(orca_obasis, atcoords, coeffsa, coeffsb):
print('5:Detected typical ORCA errors in file. Fixing them...')
result['obasis'] = orca_obasis
return
Expand All @@ -547,7 +551,7 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
print('5:Trying to fix it as if it was a file generated by PSI4 (pre 1.0).')
psi4_obasis = _fix_obasis_psi4(obasis)
if psi4_obasis is not None and \
_is_normalized_properly(psi4_obasis, atcoords, coeffs_a, coeffs_b):
_is_normalized_properly(psi4_obasis, atcoords, coeffsa, coeffsb):
print('5:Detected typical PSI4 errors in file. Fixing them...')
result['obasis'] = psi4_obasis
return
Expand All @@ -556,15 +560,15 @@ def _fix_molden_from_buggy_codes(result: dict, filename: str):
print('5:Trying to fix it as if it was a file generated by Turbomole.')
turbom_obasis = _fix_obasis_turbomole(obasis)
if turbom_obasis is not None and \
_is_normalized_properly(turbom_obasis, atcoords, coeffs_a, coeffs_b):
_is_normalized_properly(turbom_obasis, atcoords, coeffsa, coeffsb):
print('5:Detected typical Turbomole errors in file. Fixing them...')
result['obasis'] = turbom_obasis
return

# --- Renormalized contractions
print('5:Last resort: trying by renormalizing all contractions')
normed_obasis = _fix_obasis_normalize_contractions(obasis)
if _is_normalized_properly(normed_obasis, atcoords, coeffs_a, coeffs_b):
if _is_normalized_properly(normed_obasis, atcoords, coeffsa, coeffsb):
print('5:Detected unnormalized contractions in file. Fixing them...')
result['obasis'] = normed_obasis
return
Expand Down Expand Up @@ -663,28 +667,32 @@ def dump_one(f: TextIO, data: IOData):
permutation, signs = convert_conventions(obasis, CONVENTIONS)

# Print the mean-field orbitals
if data.mo.type == 'unrestricted':
if data.mo.kind == 'unrestricted':
f.write('[MO]\n')
norba = data.mo.norba
_dump_helper_orb(f, 'Alpha', data.mo.energies[:norba], data.mo.occs[:norba],
data.mo.coeffs[:, :norba][permutation] * signs.reshape(-1, 1))
_dump_helper_orb(f, 'Beta', data.mo.energies[norba:], data.mo.occs[norba:],
data.mo.coeffs[:, norba:][permutation] * signs.reshape(-1, 1))
else:
_dump_helper_orb(f, 'Alpha', data.mo.occsa,
data.mo.coeffsa[permutation] * signs.reshape(-1, 1),
data.mo.energiesa, data.mo.irrepsa)
_dump_helper_orb(f, 'Beta', data.mo.occsb,
data.mo.coeffsb[permutation] * signs.reshape(-1, 1),
data.mo.energiesb, data.mo.irrepsb)
elif data.mo.kind == 'restricted':
f.write('[MO]\n')
_dump_helper_orb(f, 'Alpha', data.mo.energies, data.mo.occs,
data.mo.coeffs[permutation] * signs.reshape(-1, 1))
_dump_helper_orb(f, 'Alpha', data.mo.occs,
data.mo.coeffs[permutation] * signs.reshape(-1, 1),
data.mo.energies, data.mo.irreps)
else:
raise NotImplementedError


def _dump_helper_orb(f, spin, orb_energies, orb_occs, orb_coeffs):
for ifn in range(orb_coeffs.shape[1]):
f.write(f' Ene= {orb_energies[ifn]:.17e}\n')
f.write(' Sym= 1a\n')
def _dump_helper_orb(f, spin, occs, coeffs, energies, irreps):
for ifn in range(coeffs.shape[1]):
f.write(f' Ene= {energies[ifn]:.17e}\n')
f.write(f' Sym= {irreps[ifn]}\n')
f.write(f' Spin= {spin}\n')
f.write(f' Occup= {orb_occs[ifn]:.17e}\n')
for ibasis in range(orb_coeffs.shape[0]):
f.write(f' Occup= {occs[ifn]:.17e}\n')
for ibasis in range(coeffs.shape[0]):
# The original molden floating-point formatting is too low
# precision. Molden also reads high-precision, so we use this
# instead.
# f.write('{:4d} {:10.6f}\n'.format(ibasis + 1, orb_coeffs[ibasis, ifn]))
f.write('{:4d} {:.17e}\n'.format(ibasis + 1, orb_coeffs[ibasis, ifn]))
f.write('{:4d} {:.17e}\n'.format(ibasis + 1, coeffs[ibasis, ifn]))
Loading

0 comments on commit ddeaa47

Please sign in to comment.