Skip to content

Commit

Permalink
Merge pull request #76 from FarnazH/molecular_orbitals
Browse files Browse the repository at this point in the history 
Molecular Orbitals
  • Loading branch information
@tovrstra
tovrstra committed Apr 24, 2019
2 parents 0671fae + 732c507 commit 0d8fd39
Show file tree
Hide file tree
Showing 11 changed files with 307 additions and 376 deletions.
29 changes: 15 additions & 14 deletions iodata/formats/cp2k.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
from scipy.special import factorialk

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


__all__ = []
Expand Down Expand Up @@ -447,24 +447,26 @@ def load(lit: LineIterator) -> Dict:

# Turn orbital data into a HORTON orbital expansions
if restricted:
mo_type = 'restricted'
norb, nel = _get_norb_nel(oe_alpha)
norb_alpha = norb_beta = norb
assert nel % 2 == 0
orb_alpha = (obasis.nbasis, norb)
orb_beta = None
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
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
orb_alpha = (obasis.nbasis, norb_alpha)
orb_alpha_coeffs = np.zeros([obasis.nbasis, norb_alpha])
orb_alpha_energies = np.zeros(norb_alpha)
orb_alpha_occs = np.zeros(norb_alpha)
orb_beta = (obasis.nbasis, norb_beta)
orb_beta_coeffs = np.zeros([obasis.nbasis, norb_beta])
orb_beta_energies = np.zeros(norb_beta)
orb_beta_occs = np.zeros(norb_beta)
Expand All @@ -473,20 +475,19 @@ def load(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)

result = {
'obasis': obasis,
'orb_alpha': orb_alpha,
'orb_alpha_coeffs': orb_alpha_coeffs,
'orb_alpha_energies': orb_alpha_energies,
'orb_alpha_occs': orb_alpha_occs,
'mo': mo,
'coordinates': obasis.centers,
'numbers': np.array([number]),
'energy': energy,
'pseudo_numbers': np.array([pseudo_number]),
}
if orb_beta is not None:
result['orb_beta'] = orb_beta
result['orb_beta_coeffs'] = orb_beta_coeffs
result['orb_beta_energies'] = orb_beta_energies
result['orb_beta_occs'] = orb_beta_occs
return result
62 changes: 32 additions & 30 deletions iodata/formats/fchk.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
import numpy as np

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


__all__ = []
Expand Down Expand Up @@ -68,8 +68,8 @@ def load(lit: LineIterator) -> Dict:
-------
out : dict
Output dictionary containing ``title``, ``coordinates``, ``numbers``, ``pseudo_numbers``,
``obasis``, ``orb_alpha``, ``energy`` & ``mulliken_charges`` keys and
corresponding values. It may also contain ``npa_charges``, ``esp_charges``, ``orb_beta``,
``obasis``, ``mo``, ``energy`` & ``mulliken_charges`` keys and
corresponding values. It may also contain ``npa_charges``, ``esp_charges``,
``dm_full_mp2``, ``dm_spin_mp2``, ``dm_full_mp3``, ``dm_spin_mp3``, ``dm_full_cc``,
``dm_spin_cc``, ``dm_full_ci``, ``dm_spin_ci``, ``dm_full_scf``, ``dm_spin_scf``,
``polar``, ``dipole_moment`` & ``quadrupole_moment`` keys and their values as well.
Expand Down Expand Up @@ -169,34 +169,36 @@ def load(lit: LineIterator) -> Dict:
nbeta = fchk['Number of beta electrons']
if nalpha < 0 or nbeta < 0 or nalpha + nbeta <= 0:
lit.error('The number of electrons is not positive.')
result['orb_alpha'] = (nbasis, nbasis_indep)
result['orb_alpha_coeffs'] = np.copy(
fchk['Alpha MO coefficients'].reshape(nbasis_indep, nbasis).T)
result['orb_alpha_energies'] = np.copy(fchk['Alpha Orbital Energies'])
aoccs = np.zeros(nbasis_indep)
aoccs[:nalpha] = 1.0
result['orb_alpha_occs'] = aoccs
if nalpha < nbeta:
raise ValueError('n_alpha={0} < n_beta={1} is not valid!'.format(nalpha, nbeta))

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'])

if 'Beta Orbital Energies' in fchk:
# UHF case
result['orb_beta'] = (nbasis, nbasis_indep)
result['orb_beta_coeffs'] = np.copy(
fchk['Beta MO coefficients'].reshape(nbasis_indep, nbasis).T)
result['orb_beta_energies'] = np.copy(fchk['Beta Orbital Energies'])
boccs = np.zeros(nbasis_indep)
boccs[:nbeta] = 1.0
result['orb_beta_occs'] = boccs

elif fchk['Number of beta electrons'] != fchk['Number of alpha electrons']:
# ROHF case
result['orb_beta'] = (nbasis, nbasis_indep)
result['orb_beta_coeffs'] = fchk['Alpha MO coefficients'].reshape(nbasis_indep, nbasis).T
result['orb_beta_energies'] = fchk['Alpha Orbital Energies']
boccs = np.zeros(nbasis_indep)
boccs[:nbeta] = 1.0
result['orb_beta_occs'] = boccs

# Delete dm_full_scf because it is known to be buggy
result.pop('dm_full_scf')
# 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_occs = np.zeros(2 * nbasis_indep)
mo_occs[:nalpha] = 1.0
mo_occs[nbasis_indep: nbasis_indep + nbeta] = 1.0
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.pop('dm_full_scf')

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

# E) Load properties
result['energy'] = fchk['Total Energy']
Expand Down
92 changes: 44 additions & 48 deletions iodata/formats/molden.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
convert_conventions, HORTON2_CONVENTIONS)
from ..periodic import sym2num, num2sym
from ..overlap import compute_overlap, gob_cart_normalization
from ..utils import angstrom, LineIterator
from ..utils import angstrom, LineIterator, MolecularOrbitals


__all__ = []
Expand Down Expand Up @@ -72,8 +72,8 @@ def load(lit: LineIterator) -> Dict:
-------
out : dict
output dictionary containing ``coordinates``, ``numbers``, ``pseudo_numbers``,
``obasis``, ``orb_alpha`` & ``signs`` keys and corresponding values. It may contain
``title`` and ``orb_beta`` keys and their values as well.
``obasis``, ``mo`` & ``signs`` keys and corresponding values. It may contain
``title`` key and its corresponding value as well.
"""
result = _load_low(lit)
Expand All @@ -94,8 +94,8 @@ def _load_low(lit: LineIterator) -> Dict:
-------
out : dict
output dictionary containing ``coordinates``, ``numbers``, ``pseudo_numbers``,
``obasis``, ``orb_alpha`` & ``signs`` keys and corresponding values. It may contain
``title`` and ``orb_beta`` keys and their values as well.
``obasis``, ``mo`` & ``signs`` keys and corresponding values. It may contain
``title`` key and its corresponding value as well.
"""
pure_angmoms = set([])
Expand Down Expand Up @@ -165,24 +165,24 @@ def _load_low(lit: LineIterator) -> Dict:
obasis.centers[:] = coordinates

if coeff_beta is None:
mo_type = 'restricted'
if coeff_alpha.shape[0] != obasis.nbasis:
lit.error("Number of alpha orbital coefficients does not match the size of the basis.")
orb_alpha = (obasis.nbasis, coeff_alpha.shape[1])
orb_alpha_coeffs = coeff_alpha
orb_alpha_energies = ener_alpha
orb_alpha_occs = occ_alpha / 2
orb_beta = None
norba = norbb = coeff_alpha.shape[1]
mo_occs = occ_alpha
mo_coeffs = coeff_alpha
mo_energy = ener_alpha
else:
mo_type = 'unrestricted'
if coeff_beta.shape[0] != obasis.nbasis:
lit.error("Number of beta orbital coefficients does not match the size of the basis.")
orb_alpha = (obasis.nbasis, coeff_alpha.shape[1])
orb_alpha_coeffs = coeff_alpha
orb_alpha_energies = ener_alpha
orb_alpha_occs = occ_alpha
orb_beta = (obasis.nbasis, coeff_beta.shape[1])
orb_beta_coeffs = coeff_beta
orb_beta_energies = ener_beta
orb_beta_occs = occ_beta
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)

# filter out ghost atoms
mask = pseudo_numbers != 0
Expand All @@ -192,21 +192,13 @@ def _load_low(lit: LineIterator) -> Dict:

result = {
'coordinates': coordinates,
'orb_alpha': orb_alpha,
'orb_alpha_coeffs': orb_alpha_coeffs,
'orb_alpha_energies': orb_alpha_energies,
'orb_alpha_occs': orb_alpha_occs,
'numbers': numbers,
'obasis': obasis,
'mo': mo,
'pseudo_numbers': pseudo_numbers,
}
if title is not None:
result['title'] = title
if orb_beta is not None:
result['orb_beta'] = orb_beta
result['orb_beta_coeffs'] = orb_beta_coeffs
result['orb_beta_energies'] = orb_beta_energies
result['orb_beta_occs'] = orb_beta_occs
return result


Expand Down Expand Up @@ -545,44 +537,47 @@ def _fix_molden_from_buggy_codes(result: Dict, filename: str):
"""
obasis = result['obasis']
if _is_normalized_properly(obasis, result['orb_alpha_coeffs'],
result.get('orb_beta_coeffs')):
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:]
else:
raise ValueError('Molecular orbital type={0} not recognized'.format(result['mo'].type))
if _is_normalized_properly(obasis, coeffs_a, coeffs_b):
# 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, result['orb_alpha_coeffs'],
result.get('orb_beta_coeffs')):
if _is_normalized_properly(orca_obasis, coeffs_a, coeffs_b):
print('5:Detected typical ORCA errors in file. Fixing them...')
result['obasis'] = orca_obasis
return

# --- PSI4
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, result['orb_alpha_coeffs'], result.get('orb_beta_coeffs')):
if psi4_obasis is not None and _is_normalized_properly(psi4_obasis, coeffs_a, coeffs_b):
print('5:Detected typical PSI4 errors in file. Fixing them...')
result['obasis'] = psi4_obasis
return

# -- Turbomole
print('5:Trying to fix it as if it was a file generated by Turbomole.')
turbomole_obasis = _fix_obasis_turbomole(obasis)
if turbomole_obasis is not None and _is_normalized_properly(
turbomole_obasis, result['orb_alpha_coeffs'], result.get('orb_beta_coeffs')):
turbom_obasis = _fix_obasis_turbomole(obasis)
if turbom_obasis is not None and _is_normalized_properly(turbom_obasis, coeffs_a, coeffs_b):
print('5:Detected typical Turbomole errors in file. Fixing them...')
result['obasis'] = turbomole_obasis
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, result['orb_alpha_coeffs'],
result.get('orb_beta_coeffs')):
if _is_normalized_properly(normed_obasis, coeffs_a, coeffs_b):
print('5:Detected unnormalized contractions in file. Fixing them...')
result['obasis'] = normed_obasis
return
Expand Down Expand Up @@ -683,24 +678,25 @@ def dump(f: TextIO, data: 'IOData'):
permutation, signs = convert_conventions(obasis, CONVENTIONS)

# Print the mean-field orbitals
if hasattr(data, 'orb_beta_coeffs'):
if data.mo.type == 'unrestricted':
f.write('[MO]\n')
_dump_helper_orb(f, 'Alpha', data.orb_alpha_energies, data.orb_alpha_occs,
data.orb_alpha_coeffs[permutation] * signs.reshape(-1, 1))
_dump_helper_orb(f, 'Beta', data.orb_beta_energies, data.orb_beta_occs,
data.orb_beta_coeffs[permutation] * signs.reshape(-1, 1))
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:
f.write('[MO]\n')
_dump_helper_orb(f, 'Alpha', data.orb_alpha_energies, data.orb_alpha_occs,
data.orb_alpha_coeffs[permutation] * signs.reshape(-1, 1), 2.0)
_dump_helper_orb(f, 'Alpha', data.mo.energies, data.mo.occs,
data.mo.coeffs[permutation] * signs.reshape(-1, 1))


def _dump_helper_orb(f, spin, orb_energies, orb_occs, orb_coeffs, occ_scale=1.0):
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')
f.write(f' Spin= {spin}\n')
f.write(f' Occup= {orb_occs[ifn] * occ_scale:.17e}\n')
f.write(f' Occup= {orb_occs[ifn]:.17e}\n')
for ibasis in range(orb_coeffs.shape[0]):
# The original molden floating-point formatting is too low
# precision. Molden also reads high-precision, so we use this
Expand Down
Loading

0 comments on commit 0d8fd39

Please sign in to comment.