from __future__ import print_function, division
from pynao.m_coulomb_am import coulomb_am
import numpy as np
[docs]def comp_coulomb_den(sv, ao_log=None, funct=coulomb_am, dtype=np.float64, **kvargs):
"""
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
Output parameters:
res: 2D np.array, dtype
matrix elements (real-space overlap) for the whole system
"""
from pynao.m_ao_matelem import ao_matelem_c
aome = ao_matelem_c(sv.ao_log.rr, sv.ao_log.pp)
me = aome.init_one_set(sv.ao_log) if ao_log is None else aome.init_one_set(ao_log)
atom2s = np.zeros((sv.natm+1), dtype=np.int32)
for atom, sp in enumerate(sv.atom2sp):
atom2s[atom+1] = atom2s[atom] + me.ao1.sp2norbs[sp]
norbs = atom2s[-1]
# dim triangular matrix: n*(n+1)/2
res = np.zeros((norbs, norbs), dtype=dtype)
for atom1, [sp1, rv1, s1, f1] in enumerate(zip(sv.atom2sp, sv.atom2coord, atom2s,
atom2s[1:])):
for atom2, [sp2, rv2, s2, f2] in enumerate(zip(sv.atom2sp, sv.atom2coord,
atom2s, atom2s[1:])):
oo2f = funct(me,sp1,rv1,sp2,rv2,**kvargs)
res[s1:f1,s2:f2] = oo2f[:,:]
return res