from __future__ import print_function, division
from pynao.m_coulomb_am import coulomb_am
import numpy as np
[docs]def comp_coulomb_pack(sv, ao_log=None, funct=coulomb_am, dtype=np.float64, **kw):
"""
Computes the matrix elements given by funct, for instance coulomb interaction
Input parameters:
sv: (System Variables), this must have arrays of coordinates and species, etc
ao_log : description of functions (either orbitals or product basis functions)
Output parameters:
res: 1D np.array dtype
matrix elements for the whole system in packed form (lower triangular part)
norbs: int
Dimension of the matrix (norbs, norbs)
"""
from pynao.m_ao_matelem import ao_matelem_c
from pynao.m_pack2den import cp_block_pack_u
aome = ao_matelem_c(sv.ao_log.rr, sv.ao_log.pp)
me = ao_matelem_c(sv.ao_log) if ao_log is None else aome.init_one_set(ao_log)
atom2s = np.zeros((sv.natm+1), dtype=np.int64)
for atom,sp in enumerate(sv.atom2sp):
atom2s[atom+1]=atom2s[atom]+me.ao1.sp2norbs[sp]
norbs = atom2s[-1]
res = np.zeros(norbs*(norbs+1)//2, dtype=dtype)
for atom1, [sp1, rv1, s1, f1] in enumerate(zip(sv.atom2sp, sv.atom2coord,
atom2s,atom2s[1:])):
#print("atom1 = {0}, rv1 = {1}".format(atom1, rv1))
for atom2, [sp2, rv2, s2, f2] in enumerate(zip(sv.atom2sp, sv.atom2coord,
atom2s,atom2s[1:])):
if atom2>atom1:
continue # skip
oo2f = funct(me, sp1, rv1, sp2, rv2, **kw)
cp_block_pack_u(oo2f, s1, f1, s2, f2, res)
return res, norbs