<commit_before>#!/usr/bin/env python
#
# Author: <NAME> <<EMAIL>>
#

'''
Gaussian and planewaves mixed density fitting
Ref:
'''

import time
import copy
import tempfile
import ctypes
import numpy
import h5py
import scipy.linalg
from pyscf import lib
from pyscf import gto
from pyscf.lib import logger
import pyscf.df
import pyscf.df.mdf
from pyscf.pbc import gto as pgto
from pyscf.pbc.df import incore
from pyscf.pbc.df import outcore
from pyscf.pbc import tools
from pyscf.pbc.df import ft_ao
from pyscf.df.mdf import _uncontract_basis
from pyscf.pbc.df import mdf_jk
from pyscf.pbc.df.mdf_jk import zdotNN, zdotCN, zdotNC
from pyscf.pbc.df import mdf_ao2mo
from pyscf.pbc.df import pwdf

KPT_DIFF_TOL = 1e-6

#
# Split the Coulomb potential to two parts.  Computing short range part in
# real space, long range part in reciprocal space.
#

#
# Be very careful with the smooth function exponents.
# The smooth function can be considered as an effective nuclear (distribution)
# potential.  This potential leads to attraction effects somewhere in the
# space.  The attraction potential can interact with diffuse AO basis, which
# may artificially produce a dipole
#
# A relatively steep function is local in space thus maybe have smaller region
# of attraction potential.  An atom-specific eta based on atomic radius might
# be needed
#
def estimate_eta(cell, cutoff=1e-12):
    '''The exponent of the smooth gaussian model density, requiring that at
    boundary, density ~ 4pi rmax^2 exp(-eta*rmax^2) ~ 1e-12
    '''
    rmax = max(lib.norm(cell.lattice_vectors(), axis=0))
    eta = max(-numpy.log(cutoff/(4*numpy.pi*rmax**2))/rmax**2, .1)
    return eta

def make_modrho_basis(cell, auxbasis=None, drop_eta=1.):
    auxcell = copy.copy(cell)
    if auxbasis is None:
        _basis = _uncontract_basis(cell, auxbasis)
    elif isinstance(auxbasis, str):
        uniq_atoms = set([a[0] for a in cell._atom])
        _basis = auxcell.format_basis(dict([(a, auxbasis) for a in uniq_atoms]))
    else:
        _basis = auxcell.format_basis(auxbasis)
    auxcell._basis = _basis
    auxcell._atm, auxcell._bas, auxcell._env = \
            auxcell.make_env(cell._atom, auxcell._basis, cell._env[:gto.PTR_ENV_START])

# Note libcint library will multiply the norm of the integration over spheric
# part sqrt(4pi) to the basis.
    half_sph_norm = numpy.sqrt(.25/numpy.pi)
    steep_shls = []
    ndrop = 0
    for ib in range(len(auxcell._bas)):
        l = auxcell.bas_angular(ib)
        np = auxcell.bas_nprim(ib)
        nc = auxcell.bas_nctr(ib)
        es = auxcell.bas_exp(ib)
        ptr = auxcell._bas[ib,gto.PTR_COEFF]
        cs = auxcell._env[ptr:ptr+np*nc].reshape(nc,np).T

#        if l > 0:
#            continue
        if numpy.any(es < drop_eta):
            es = es[es>=drop_eta]
            cs = cs[es>=drop_eta]
            np, ndrop = len(es), ndrop+np-len(es)
            pe = auxcell._bas[ib,gto.PTR_EXP]
            auxcell._bas[ib,gto.NPRIM_OF] = np
            auxcell._env[pe:pe+np] = es

        if np > 0:
# int1 is the multipole value. l*2+2 is due to the radial part integral
# \int (r^l e^{-ar^2} * Y_{lm}) (r^l Y_{lm}) r^2 dr d\Omega
            int1 = gto.mole._gaussian_int(l*2+2, es)
            s = numpy.einsum('pi,p->i', cs, int1)
# The auxiliary basis normalization factor is not a must for density expansion.<commit_msg>Bug fix for pbc.df.mdf auxbasis generation<commit_after>#!/usr/bin/env python
#
# Author: <NAME> <<EMAIL>>
#

'''
Gaussian and planewaves mixed density fitting
Ref:
'''

import time
import copy
import tempfile
import ctypes
import numpy
import h5py
import scipy.linalg
from pyscf import lib
from pyscf import gto
from pyscf.lib import logger
import pyscf.df
import pyscf.df.mdf
from pyscf.pbc import gto as pgto
from pyscf.pbc.df import incore
from pyscf.pbc.df import outcore
from pyscf.pbc import tools
from pyscf.pbc.df import ft_ao
from pyscf.df.mdf import _uncontract_basis
from pyscf.pbc.df import mdf_jk
from pyscf.pbc.df.mdf_jk import zdotNN, zdotCN, zdotNC
from pyscf.pbc.df import mdf_ao2mo
from pyscf.pbc.df import pwdf

KPT_DIFF_TOL = 1e-6

#
# Split the Coulomb potential to two parts.  Computing short range part in
# real space, long range part in reciprocal space.
#

#
# Be very careful with the smooth function exponents.
# The smooth function can be considered as an effective nuclear (distribution)
# potential.  This potential leads to attraction effects somewhere in the
# space.  The attraction potential can interact with diffuse AO basis, which
# may artificially produce a dipole
#
# A relatively steep function is local in space thus maybe have smaller region
# of attraction potential.  An atom-specific eta based on atomic radius might
# be needed
#
def estimate_eta(cell, cutoff=1e-12):
    '''The exponent of the smooth gaussian model density, requiring that at
    boundary, density ~ 4pi rmax^2 exp(-eta*rmax^2) ~ 1e-12
    '''
    rmax = max(lib.norm(cell.lattice_vectors(), axis=0))
    eta = max(-numpy.log(cutoff/(4*numpy.pi*rmax**2))/rmax**2, .1)
    return eta

def make_modrho_basis(cell, auxbasis=None, drop_eta=1.):
    auxcell = copy.copy(cell)
    if auxbasis is None:
        _basis = _uncontract_basis(cell, auxbasis)
    elif isinstance(auxbasis, str):
        uniq_atoms = set([a[0] for a in cell._atom])
        _basis = auxcell.format_basis(dict([(a, auxbasis) for a in uniq_atoms]))
    else:
        _basis = auxcell.format_basis(auxbasis)
    auxcell._basis = _basis
    auxcell._atm, auxcell._bas, auxcell._env = \
            auxcell.make_env(cell._atom, auxcell._basis, cell._env[:gto.PTR_ENV_START])

# Note libcint library will multiply the norm of the integration over spheric
# part sqrt(4pi) to the basis.
    half_sph_norm = numpy.sqrt(.25/numpy.pi)
    steep_shls = []
    ndrop = 0
    for ib in range(len(auxcell._bas)):
        l = auxcell.bas_angular(ib)
        np = auxcell.bas_nprim(ib)
        nc = auxcell.bas_nctr(ib)
        es = auxcell.bas_exp(ib)
        ptr = auxcell._bas[ib,gto.PTR_COEFF]
        cs = auxcell._env[ptr:ptr+np*nc].reshape(nc,np).T

#        if l > 0:
#            continue
        if numpy.any(es < drop_eta):
            cs = cs[es>=drop_eta]
            es = es[es>=drop_eta]
            np, ndrop = len(es), ndrop+np-len(es)
            pe = auxcell._bas[ib,gto.PTR_EXP]
            auxcell._bas[ib,gto.NPRIM_OF] = np
            auxcell._env[pe:pe+np] = es

        if np > 0:
# int1 is the multipole value. l*2+2 is due to the radial part integral
# \int (r^l e^{-ar^2} * Y_{lm}) (r^l Y_{lm}) r^2 dr d\Omega
            int1 = gto.mole._gaussian_int(l*2+2, es)
            s = numpy.einsum('pi,p->i', cs, int1)
# The auxiliary basis normalization factor is not a must for density expansion.
{'fraction_non_alphanumeric': 0.10336146897675107, 'fraction_numerical': 0.010365763364430623, 'mean_word_length': 3.184634448574969, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4758181', 'n_tokens_mistral': 2456, 'n_tokens_neox': 2355, 'n_words': 704}