Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions qstack/spahm/compute_spahm.py
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ def get_guess_orbitals_grad(mol, guess, field=None):
return e, eigenvalue_grad(mol, e, c, s1, h1), de_dfield


def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
def get_spahm_representation(mol, guess_in, xc="pbe", field=None, verbose=1):
"""Compute the ε-SPAHM molecular representation.

Reference:
Expand All @@ -128,6 +128,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
xc (str): Exchange-correlation functional name. Defaults to 'pbe'.
field (numpy ndarray, optional): 3-component uniform electric field (∇φ) in atomic units.
Defaults to None.
verbose (int): Level of verbosity

Returns:
numpy ndarray: SPAHM representation consisting of occupied orbital eigenvalues.
Expand All @@ -136,7 +137,7 @@ def get_spahm_representation(mol, guess_in, xc="pbe", field=None):
"""
guess = get_guess(guess_in)
e, _v = get_guess_orbitals(mol, guess, xc, field=field)
e1 = get_occ(e, mol.nelec, mol.spin)
e1 = get_occ(e, mol.nelec, mol.spin, verbose=verbose)
return e1


Expand Down
16 changes: 9 additions & 7 deletions qstack/spahm/guesses.py
Original file line number Diff line number Diff line change
Expand Up @@ -197,13 +197,13 @@ def get_guess(arg):
return guesses_dict[arg]


def check_nelec(nelec, nao):
def check_nelec(nelec, nao, full_shell_warning=1):
"""Validate that the number of electrons can be accommodated by available orbitals.

Args:
nelec (tuple or int): Number of electrons (alpha, beta) or total.
nao (int): Number of atomic orbitals.

full_shell_warning (int): if < 1 supresses warnings about fully occupied shells
Raises:
RuntimeError: If there are more electrons than available orbitals.

Expand All @@ -212,26 +212,27 @@ def check_nelec(nelec, nao):
"""
if np.any(np.array(nelec) > nao):
raise RuntimeError(f'Too many electrons ({nelec}) for {nao} orbitals')
elif np.any(np.array(nelec) == nao):
elif np.any(np.array(nelec) == nao) and full_shell_warning>0:
msg = f'{nelec} electrons for {nao} orbitals. Is the input intended to have a complete shell?'
warnings.warn(msg, RuntimeWarning, stacklevel=2)


def get_occ(e, nelec, spin):
def get_occ(e, nelec, spin, verbose=1):
"""Extract occupied orbital eigenvalues/energies.

Args:
e (numpy ndarray): Full array of orbital eigenvalues (1D)
or possibly arrays of larger dimensionality.
nelec (tuple): Number of (alpha, beta) electrons.
spin (int or None): Spin multiplicity. If None, assumes closed-shell.
verbose (int): Level of verbosity (< 1 supresses full-shell warnings)

Returns:
numpy ndarray: Occupied eigenvalues. Shape depends on spin:
- Closed-shell (spin=None): 1D array of occupied eigenvalues
- Open-shell: 2D array (2, nocc) for alpha and beta separately
"""
check_nelec(nelec, e.shape[0])
check_nelec(nelec, e.shape[0], full_shell_warning=verbose)
if spin is None:
nocc = nelec[0]
return e[:nocc,...]
Expand All @@ -243,20 +244,21 @@ def get_occ(e, nelec, spin):
return e1


def get_dm(v, nelec, spin):
def get_dm(v, nelec, spin, verbose=1):
"""Construct density matrix from occupied molecular orbitals.

Args:
v (numpy ndarray): 2D array of MO coefficients (eigenvectors), columns are MOs.
nelec (tuple): Number of (alpha, beta) electrons.
spin (int or None): Spin multiplicity. If None, assumes closed-shell (RHF).
verbose (int): Level of verbosity (< 1 supresses full-shell warnings)

Returns:
numpy ndarray: Density matrix in AO basis.
- Closed-shell: 2D array (nao, nao)
- Open-shell: 3D array (2, nao, nao) for alpha and beta
"""
check_nelec(nelec, len(v))
check_nelec(nelec, len(v), full_shell_warning=verbose)
if spin is None:
nocc = nelec[0]
dm = v[:,:nocc] @ v[:,:nocc].T
Expand Down
4 changes: 2 additions & 2 deletions qstack/spahm/rho/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=None, readdm=None, pri
xc (str): Exchange-correlation functional for guess. Defaults to defaults.xc.
spin (list or None): List of spin multiplicities. Defaults to None.
readdm (str, optional): Directory path to load pre-computed density matrices. Defaults to None.
printlevel (int): Verbosity level. Defaults to 0.
printlevel (int): Verbosity level. Defaults to 0. (if < 0 supresses warnings)

Returns:
list: List of density matrices (2D or 3D numpy arrays).
Expand All @@ -120,7 +120,7 @@ def mols_guess(mols, xyzlist, guess, xc=defaults.xc, spin=None, readdm=None, pri
print(xyzfile, flush=True)
if readdm is None:
_e, v = spahm.get_guess_orbitals(mol, guess, xc=xc)
dm = guesses.get_dm(v, mol.nelec, mol.spin if sp is not None else None)
dm = guesses.get_dm(v, mol.nelec, mol.spin if sp is not None else None, verbose=printlevel+1)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@briling I'm not a big fan of the increment startegy here, but seems to be the best option...
what do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think it's okay!

else:
dm = np.load(f'{readdm}/{os.path.basename(xyzfile)}.npy')
if spin and dm.ndim==2:
Expand Down
Loading