diff --git a/examples/dft2/24-reparameterize_xc_functional.py b/examples/dft2/24-reparameterize_xc_functional.py deleted file mode 100644 index e5a90faf8..000000000 --- a/examples/dft2/24-reparameterize_xc_functional.py +++ /dev/null @@ -1,71 +0,0 @@ -#!/usr/bin/env python - -''' -Running DFT calculations with reparameterized functionals in Libxc. - -See also -* Example 24-custom_xc_functional.py to customize XC functionals using the - functionals provided by Libxc library. -''' - -from pyscf import gto -from pyscf import dft, dft2 -import numpy as np - -# Patch the libxc module with the new implementation -dft.libxc = dft2.libxc -dft.numint.libxc = dft2.libxc -dft.numint.LibXCMixin.libxc = dft2.libxc - -XC_ID_B97_2 = 410 - -param_b97_1 = np.array([0.789518, 0.573805, 0.660975, 0.0, 0.0, - 0.0820011, 2.71681, -2.87103, 0.0, 0.0, - 0.955689, 0.788552, -5.47869, 0.0, 0.0, - 0.21 -]) - -param_b97_3 = np.array([0.7334648, 0.292527, 3.338789, -10.51158, 10.60907, - 0.5623649, -1.32298, 6.359191, -7.464002, 1.827082, - 1.13383, -2.811967, 7.431302, -1.969342, -11.74423, - 2.692880E-01 -]) - -mol = gto.M( - atom = ''' - O 0. 0. 0. - H 0. -0.757 0.587 - H 0. 0.757 0.587 ''', - basis = 'ccpvdz') - -# Run normal B97-1 -print('Normal B97-1') -mf = dft.RKS(mol, 'B97-1') -e_b971 = mf.kernel() - -# Run normal B97-2 -print('\nNormal B97-2') -mf.xc = 'B97-2' -e_b972 = mf.kernel() - -# Run normal B97-3 -print('\nNormal B97-3') -mf.xc = 'B97-3' -e_b973 = mf.kernel() - -# Construct XC named `myfunctional` based on B97-2, but set its parameter to be B97-1 -print('\nReparameterized B97-2: will be the same as B97-1') -dft2.libxc.register_custom_functional_('myfunctional', 'B97-2', - ext_params={XC_ID_B97_2: param_b97_1}) -mf.xc = 'myfunctional' -e = mf.kernel() -print('difference:', e - e_b971) - -# Update the parameter of `myfunctional` to be the same as B97-3 and rerun -print('\nReparameterized B97-2: will be the same as B97-3') -dft2.libxc.update_custom_functional_('myfunctional', - ext_params={XC_ID_B97_2: param_b97_3}) -e = mf.kernel() -print('difference:', e - e_b973) -print() - diff --git a/examples/mcpdft/03-metaGGA_functionals.py b/examples/mcpdft/03-metaGGA_functionals.py deleted file mode 100644 index 574ef0c94..000000000 --- a/examples/mcpdft/03-metaGGA_functionals.py +++ /dev/null @@ -1,25 +0,0 @@ -# patch dft2.libxc (must be done before loading mpcdft module) -from pyscf import dft -from pyscf import dft2 -dft.libxc = dft2.libxc -#!/usr/bin/env/python -from pyscf import gto, scf, mcpdft - -mol = gto.M ( - atom = 'O 0 0 0; O 0 0 1.2', - basis = 'ccpvdz', - spin = 2) - -mf = scf.RHF (mol).run () - -# MC23: meta-hybrid on-top functional [PNAS, 122, 1, 2025, e2419413121; https://doi.org/10.1073/pnas.2419413121] - -# State-Specific -mc = mcpdft.CASCI(mf, 'MC23', 6, 8) -mc.kernel() - -# State-average -nroots=2 -mc = mcpdft.CASCI(mf, 'MC23', 6, 8) -mc.fcisolver.nroots=nroots -mc.kernel()[0] diff --git a/pyscf/dft2/__init__.py b/pyscf/dft2/__init__.py deleted file mode 100644 index e4d222da9..000000000 --- a/pyscf/dft2/__init__.py +++ /dev/null @@ -1,21 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2022 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# - -# NOTE 08/04/2025 -# This module should function with PySCF v2.10.0 specifically -# The next release will include this module, which should then be removed -# from PySCF-Forge -from pyscf.dft2 import libxc diff --git a/pyscf/dft2/libxc.py b/pyscf/dft2/libxc.py deleted file mode 100644 index 4b0e40606..000000000 --- a/pyscf/dft2/libxc.py +++ /dev/null @@ -1,1333 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- -# Copyright 2014-2021 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Authors: Qiming Sun -# Susi Lehtola -# Dayou Zhang -# This file is adapted from `dft/libxc.py` of the PySCF core module -# commit 25eaa9572977b903de24d5c11ad345cecd744728 -# -# Implemented functions register_custom_functional_, update_custom_functional_, -# and unregister_custom_functional_. To make this new functionality possible, -# and to improve the memory management of the LibXC data structure, interface -# functions which originally took arguments of type c_int were modified to take -# arguments of type xc_func_type from LibXC. The LibXC data structure are cached -# and managed in the python layer to avoid construction and destruction at every -# API call. - -''' -XC functional, the interface to libxc -(https://libxc.gitlab.io) -''' - -import sys -import warnings -import ctypes -import math -import numpy -from functools import lru_cache, cached_property -from pyscf import lib -from pyscf.dft.xc.utils import remove_dup, format_xc_code -from pyscf.dft import xc_deriv -from pyscf.dft.libxc import XC_CODES, XC, PROBLEMATIC_XC, XC_KEYS, XC_ALIAS, _NAME_WITH_DASH -from pyscf import __config__ - -_itrf = lib.load_library('libxc_itrf2') -_itrf.LIBXC_is_lda.restype = ctypes.c_int -_itrf.LIBXC_is_gga.restype = ctypes.c_int -_itrf.LIBXC_is_meta_gga.restype = ctypes.c_int -_itrf.LIBXC_needs_laplacian.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_needs_laplacian.restype = ctypes.c_int -_itrf.LIBXC_is_hybrid.restype = ctypes.c_int -_itrf.LIBXC_is_nlc.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_is_nlc.restype = ctypes.c_int -_itrf.LIBXC_is_cam_rsh.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_is_cam_rsh.restype = ctypes.c_int -_itrf.LIBXC_xc_type.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_xc_type.restype = ctypes.c_int -_itrf.LIBXC_max_deriv_order.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_max_deriv_order.restype = ctypes.c_int -_itrf.LIBXC_hybrid_coeff.argtypes = (ctypes.c_void_p, ) -_itrf.LIBXC_hybrid_coeff.restype = ctypes.c_double -_itrf.LIBXC_nlc_coeff.argtypes = (ctypes.c_void_p, ctypes.POINTER(ctypes.c_double)) -_itrf.LIBXC_rsh_coeff.argtypes = (ctypes.c_void_p, ctypes.POINTER(ctypes.c_double)) - -_itrf.LIBXC_xc_reference.argtypes = (ctypes.c_int, ctypes.c_void_p, ctypes.POINTER(ctypes.c_char_p)) - -_itrf.LIBXC_xc_func_init.restype = ctypes.c_void_p -_itrf.LIBXC_xc_func_end.argtypes = (ctypes.c_int, ctypes.c_void_p) -_itrf.LIBXC_xc_func_type_size.restype = ctypes.c_size_t -_itrf.LIBXC_eval_xc.argtypes = (ctypes.c_int, ctypes.c_void_p, - ctypes.POINTER(ctypes.c_double), - ctypes.c_int, ctypes.c_int, - ctypes.c_int, ctypes.c_int, - ctypes.c_int, - ctypes.c_void_p, - ctypes.c_void_p) -_itrf.LIBXC_xc_func_set_params.argtypes = (ctypes.c_int, ctypes.c_void_p, - ctypes.POINTER(ctypes.c_double), - ctypes.c_double) - -_itrf.xc_version_string.restype = ctypes.c_char_p -_itrf.xc_reference.restype = ctypes.c_char_p -_itrf.xc_reference_doi.restype = ctypes.c_char_p -_itrf.xc_number_of_functionals.restype = ctypes.c_int -_itrf.xc_available_functional_numbers.argtypes = (numpy.ctypeslib.ndpointer(dtype=numpy.intc, - ndim=1, flags=("W", "C", "A")), ) -_itrf.xc_functional_get_name.argtypes = (ctypes.c_int, ) -_itrf.xc_functional_get_name.restype = ctypes.c_char_p -_itrf.xc_functional_get_number.argtypes = (ctypes.c_char_p, ) -_itrf.xc_functional_get_number.restype = ctypes.c_int -_itrf.xc_func_get_info.argtypes = (ctypes.c_void_p, ) -_itrf.xc_func_get_info.restype = ctypes.c_void_p -_itrf.xc_func_info_get_n_ext_params.argtypes = (ctypes.c_void_p, ) -_itrf.xc_func_info_get_n_ext_params.restype = ctypes.c_int -_itrf.xc_func_set_ext_params.argtypes = (ctypes.c_void_p, ctypes.POINTER(ctypes.c_double)) - -_XC_FUNC_TYPE_SIZE = _itrf.LIBXC_xc_func_type_size() - -def libxc_version(): - '''Returns the version of libxc''' - return _itrf.xc_version_string().decode("UTF-8") - -def libxc_reference(): - '''Returns the reference to libxc''' - return _itrf.xc_reference().decode("UTF-8") - -def libxc_reference_doi(): - '''Returns the reference to libxc''' - return _itrf.xc_reference_doi().decode("UTF-8") - -__version__ = libxc_version() -__reference__ = libxc_reference() -__reference_doi__ = libxc_reference_doi() - -def available_libxc_functionals(): - # Number of functionals is - nfunc = _itrf.xc_number_of_functionals() - # Get functional numbers - func_ids = numpy.zeros(nfunc, dtype=numpy.intc) - _itrf.xc_available_functional_numbers(func_ids) - # Returned array - return {_itrf.xc_functional_get_name(x).decode("UTF-8").upper() : x for x in func_ids} - -def xc_reference(xc_code): - '''Returns the reference to the individual XC functional''' - xc = _get_xc(xc_code) - return _xc_reference(xc.nfunc, xc.xc_arr) - -def _xc_reference(nfunc, xc_arr): - refs = [] - c_refs = (8 * nfunc * ctypes.c_char_p)() - _itrf.LIBXC_xc_reference(nfunc, xc_arr, c_refs) - for ref in c_refs: - if ref: - refs.append(ref.decode("UTF-8")) - return refs - -def xc_type(xc_code): - '''Returns the type of a functional as a str - ''' - xc = _get_xc(xc_code) - return xc.xc_type - -def _xc_type(nfunc, xc_arr): - if nfunc == 0: - return 'HF' - else: - t = _itrf.LIBXC_xc_type(nfunc, xc_arr) - return ('LDA', 'GGA', 'MGGA', 'UNKNOWN')[t] - -def is_lda(xc_code): - return xc_type(xc_code) == 'LDA' - -def is_hybrid_xc(xc_code): - if xc_code is None: - return False - elif isinstance(xc_code, str): - if 'HF' in xc_code: - return True - xc = _get_xc(xc_code) - return _is_hybrid_xc(xc.xc_objs, xc.hyb, xc.facs, xc_code) - elif numpy.issubdtype(type(xc_code), numpy.integer): - xc = _get_xc(xc_code) - return _is_hybrid_xc(xc.xc_objs, xc.hyb, xc.facs, xc_code) - else: - return any((is_hybrid_xc(x) for x in xc_code)) - -def _is_hybrid_xc(xc_objs, hyb, facs, xc_code): - if _hybrid_coeff(xc_objs, hyb, facs) != 0: - return True - if _rsh_coeff(xc_objs, hyb, facs, xc_code) != (0, 0, 0): - return True - return False - -def is_meta_gga(xc_code): - '''Returns True if a functional is a meta GGA - ''' - return xc_type(xc_code) == 'MGGA' - -def is_gga(xc_code): - '''Returns True if a functional is a GGA - ''' - return xc_type(xc_code) == 'GGA' - -def is_nlc(xc_code): - # identify nlc by xc_code itself if enable_nlc is None - if isinstance(xc_code, str) or \ - numpy.issubdtype(type(xc_code), numpy.integer): - xc = _get_xc(xc_code) - return xc.is_nlc - else: - return any((is_nlc(x) for x in xc_code)) - -def _is_nlc(nfunc, xc_arr): - return _itrf.LIBXC_is_nlc(nfunc, xc_arr) - -def needs_laplacian(xc_code): - xc = _get_xc(xc_code) - return xc.needs_laplacian - -def _needs_laplacian(nfunc, xc_arr): - return _itrf.LIBXC_needs_laplacian(nfunc, xc_arr) != 0 - -def max_deriv_order(xc_code): - xc = _get_xc(xc_code) - return xc.max_deriv_order - -def _max_deriv_order(nfunc, xc_arr): - if nfunc == 0: - return 3 - else: - return _itrf.LIBXC_max_deriv_order(nfunc, xc_arr) - -def test_deriv_order(xc_code, deriv, raise_error=False): - support = deriv <= max_deriv_order(xc_code) - if not support and raise_error: - from pyscf.dft import xcfun - msg = ('libxc library does not support derivative order %d for %s' % - (deriv, xc_code)) - try: - if xcfun.test_deriv_order(xc_code, deriv, raise_error=False): - msg += (''' - This functional derivative is supported in the xcfun library. - The following code can be used to change the libxc library to xcfun library: - - from pyscf.dft import xcfun - mf._numint.libxc = xcfun -''') - raise NotImplementedError(msg) - except KeyError as e: - sys.stderr.write('\n'+msg+'\n') - sys.stderr.write('%s not found in xcfun library\n\n' % xc_code) - raise e - return support - -def hybrid_coeff(xc_code, spin=0): - '''Support recursively defining hybrid functional - ''' - xc = _get_xc(xc_code, spin=spin) - return xc.hybrid_coeff - -def _hybrid_coeff(xc_objs, hyb, facs): - hybs = (fac * _itrf.LIBXC_hybrid_coeff(xc) for fac, xc in zip(facs, xc_objs)) - return hyb[0] + sum(hybs) - -def nlc_coeff(xc_code): - '''Get NLC coefficients - ''' - xc = _get_xc(xc_code) - return xc.nlc_coeff - -def _nlc_coeff(xc_objs, hyb, facs): - nlc_pars = [] - nlc_tmp = (ctypes.c_double*2)() - for fac, xc in zip(facs, xc_objs): - if _itrf.LIBXC_is_nlc(ctypes.c_int(1), xc): - _itrf.LIBXC_nlc_coeff(xc, nlc_tmp) - nlc_pars.append((tuple(nlc_tmp), fac)) - return tuple(nlc_pars) - -def rsh_coeff(xc_code): - '''Range-separated parameter and HF exchange components: omega, alpha, beta - - Exc_RSH = c_LR * LR_HFX + c_SR * SR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec - = alpha * HFX + beta * SR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec - = alpha * LR_HFX + hyb * SR_HFX + (1-c_SR) * Ex_SR + (1-c_LR) * Ex_LR + Ec - - SR_HFX = < pi | (1-erf(-omega r_{12}))/r_{12} | iq > - LR_HFX = < pi | erf(-omega r_{12})/r_{12} | iq > - alpha = c_LR - beta = c_SR - c_LR = hyb - alpha - ''' - if xc_code is None: - return 0, 0, 0 - - xc = _get_xc(xc_code) - return xc.rsh_coeff - -def _rsh_coeff(xc_objs, hyb, facs, xc_code): - check_omega = True - if isinstance(xc_code, str) and ',' in xc_code: - # Parse only X part for the RSH coefficients. This is to handle - # exceptions for C functionals such as M11. - xc_code = format_xc_code(xc_code) - xc_code = xc_code.split(',')[0] + ',' - if 'SR_HF' in xc_code or 'LR_HF' in xc_code or 'RSH(' in xc_code: - check_omega = False - - hyb, alpha, omega = hyb - if omega == 0: - # SR and LR Coulomb share the same coefficients - # Note: this change breaks compatibility with pyscf-2.7 - assert hyb == alpha - beta = 0. - else: - beta = hyb - alpha - rsh_pars = [omega, alpha, beta] - rsh_tmp = (ctypes.c_double*3)() - for fac, xc in zip(facs, xc_objs): - _itrf.LIBXC_rsh_coeff(xc, rsh_tmp) - if rsh_pars[0] == 0: - rsh_pars[0] = rsh_tmp[0] - elif check_omega: - # Check functional is actually a CAM functional - if rsh_tmp[0] != 0 and not _itrf.LIBXC_is_cam_rsh(xc): - raise KeyError('Libxc functional %i employs a range separation ' - 'kernel that is not supported in PySCF' % xc.xc_code) - # Check omega - if (rsh_tmp[0] != 0 and rsh_pars[0] != rsh_tmp[0]): - raise ValueError('Different values of omega found for RSH functionals') - rsh_pars[1] += rsh_tmp[1] * fac - rsh_pars[2] += rsh_tmp[2] * fac - return tuple(rsh_pars) - -def parse_xc_name(xc_name='LDA,VWN'): - '''Convert the XC functional name to libxc library internal ID. - ''' - fn_facs = parse_xc(xc_name)[1] - return fn_facs[0][0], fn_facs[1][0] - -def parse_xc(description): - r'''Rules to input functional description: - - * The given functional description must be a one-line string. - * The functional description is case-insensitive. - * The functional description string has two parts, separated by ",". The - first part describes the exchange functional, the second is the correlation - functional. - - - If "," was not in string, the entire string is considered as a - compound XC functional (including both X and C functionals, such as b3lyp). - - To input only X functional (without C functional), leave the second - part blank. E.g. description='slater,' means pure LDA functional. - - To neglect X functional (just apply C functional), leave the first - part blank. E.g. description=',vwn' means pure VWN functional. - - If compound XC functional is specified, no matter whether it is in the - X part (the string in front of comma) or the C part (the string behind - comma), both X and C functionals of the compound XC functional will be - used. - - * The functional name can be placed in arbitrary order. Two name needs to - be separated by operators "+" or "-". Blank spaces are ignored. - NOTE the parser only reads operators "+" "-" "*". / is not in support. - * A functional name can have at most one factor. If the factor is not - given, it is set to 1. Compound functional can be scaled as a unit. For - example '0.5*b3lyp' is equivalent to - 'HF*0.1 + .04*LDA + .36*B88, .405*LYP + .095*VWN' - * String "HF" stands for exact exchange (HF K matrix). Putting "HF" in - correlation functional part is the same to putting "HF" in exchange - part. - * String "RSH" means range-separated operator. Its format is - RSH(omega, alpha, beta). Another way to input RSH is to use keywords - SR_HF and LR_HF: "SR_HF(0.1) * alpha_plus_beta" and "LR_HF(0.1) * - alpha" where the number in parenthesis is the value of omega. - * Be careful with the libxc convention on GGA functional, in which the LDA - contribution has been included. - - Args: - description : str - A string to describe the linear combination of different XC functionals. - The X and C functional are separated by comma like '.8*LDA+.2*B86,VWN'. - If "HF" was appeared in the string, it stands for the exact exchange. - - Returns: - decoded XC description, with the data structure - (hybrid, alpha, omega), ((libxc-Id, fac), (libxc-Id, fac), ...) - ''' # noqa: E501 - - hyb = [0, 0, 0] # hybrid, alpha, omega (== SR_HF, LR_HF, omega) - if description is None: - return tuple(hyb), () - elif numpy.issubdtype(type(description), numpy.integer): - return tuple(hyb), ((description, 1.),) - elif not isinstance(description, str): #isinstance(description, (tuple,list)): - return parse_xc('%s,%s' % tuple(description)) - - description = description.upper() - if '-D3' in description or '-D4' in description: - from pyscf.scf.dispersion import parse_dft - description, _, _ = parse_dft(description) - description = description.upper() - - if (description in ('B3P86', 'B3LYP', 'X3LYP') and - not getattr(parse_xc, 'b3lyp5_warned', False) and - not hasattr(__config__, 'B3LYP_WITH_VWN5')): - parse_xc.b3lyp5_warned = True - warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, ' - 'corresponding to the original definition by Stephens et al. (issue 1480) ' - 'and the same as the B3LYP functional in Gaussian. ' - 'To restore the VWN5 definition, you can put the setting ' - '"B3LYP_WITH_VWN5 = True" in pyscf_conf.py') - - def assign_omega(omega, hyb_or_sr, lr=0): - if hyb[2] == omega or omega == 0: - hyb[0] += hyb_or_sr - hyb[1] += lr - elif hyb[2] == 0: - hyb[0] += hyb_or_sr - hyb[1] += lr - hyb[2] = omega - else: - raise ValueError('Different values of omega found for RSH functionals') - fn_facs = [] - - def parse_token(token, ftype, search_xc_alias=False): - if token: - if token[0] == '-': - sign = -1 - token = token[1:] - else: - sign = 1 - if '*' in token: - fac, key = token.split('*') - if fac[0].isalpha(): - fac, key = key, fac - fac = sign * float(fac.replace('E_', 'E-')) - else: - fac, key = sign, token - - if key[:3] == 'RSH': - # RSH(alpha; beta; omega): Range-separated-hybrid functional - # See also utils.format_xc_code - alpha, beta, omega = [float(x) for x in key[4:-1].split(';')] - assign_omega(omega, fac*(alpha+beta), fac*alpha) - elif key == 'HF': - hyb[0] += fac - hyb[1] += fac # also add to LR_HF - elif 'SR_HF' in key: - if '(' in key: - omega = float(key.split('(')[1].split(')')[0]) - assign_omega(omega, fac, 0) - else: # Assuming this omega the same to the existing omega - hyb[0] += fac - elif 'LR_HF' in key: - if '(' in key: - omega = float(key.split('(')[1].split(')')[0]) - assign_omega(omega, 0, fac) - else: - hyb[1] += fac # == alpha - elif key.isdigit(): - fn_facs.append((int(key), fac)) - else: - if search_xc_alias and key in XC_ALIAS: - x_id = XC_ALIAS[key] - elif key in XC_CODES: - x_id = XC_CODES[key] - else: - possible_xc_for = fpossible_dic[ftype] - possible_xc = XC_KEYS.intersection(possible_xc_for(key)) - if possible_xc: - if len(possible_xc) > 1: - sys.stderr.write('Possible xc_code %s matches %s. ' - % (list(possible_xc), key)) - for x_id in possible_xc: # Prefer X functional - if '_X_' in x_id: - break - else: - x_id = possible_xc.pop() - sys.stderr.write('XC parser takes %s\n' % x_id) - sys.stderr.write('You can add prefix to %s for a ' - 'specific functional (e.g. GGA_X_%s, ' - 'HYB_MGGA_X_%s)\n' - % (key, key, key)) - else: - x_id = possible_xc.pop() - x_id = XC_CODES[x_id] - else: - # Some libxc functionals may not be listed in the - # XC_CODES table. Query libxc directly - x_id = _itrf.xc_functional_get_number(ctypes.c_char_p(key.encode())) - if x_id == -1: - raise KeyError(f"LibXCFunctional: name '{key}' not found.") - if isinstance(x_id, str): - hyb1, fn_facs1 = parse_xc(x_id) - # Recursively scale the composed functional, to support e.g. '0.5*b3lyp' - if hyb1[0] != 0 or hyb1[1] != 0: - assign_omega(hyb1[2], hyb1[0]*fac, hyb1[1]*fac) - fn_facs.extend([(xid, c*fac) for xid, c in fn_facs1]) - elif x_id is None: - raise NotImplementedError('%s functional %s' % (ftype, key)) - else: - fn_facs.append((x_id, fac)) - def possible_x_for(key): - return {key, - 'LDA_X_'+key, 'GGA_X_'+key, 'MGGA_X_'+key, - 'HYB_GGA_X_'+key, 'HYB_MGGA_X_'+key} - def possible_xc_for(key): - return {key, 'LDA_XC_'+key, 'GGA_XC_'+key, 'MGGA_XC_'+key, - 'HYB_LDA_XC_'+key, 'HYB_GGA_XC_'+key, 'HYB_MGGA_XC_'+key} - def possible_k_for(key): - return {key, - 'LDA_K_'+key, 'GGA_K_'+key,} - def possible_x_k_for(key): - return possible_x_for(key).union(possible_k_for(key)) - def possible_c_for(key): - return {key, - 'LDA_C_'+key, 'GGA_C_'+key, 'MGGA_C_'+key} - fpossible_dic = {'X': possible_x_for, - 'C': possible_c_for, - 'compound XC': possible_xc_for, - 'K': possible_k_for, - 'X or K': possible_x_k_for} - - description = format_xc_code(description) - - if '-' in description: # To handle e.g. M06-L - for key in _NAME_WITH_DASH: - if key in description: - description = description.replace(key, _NAME_WITH_DASH[key]) - - if ',' in description: - x_code, c_code = description.split(',') - for token in x_code.replace('-', '+-').replace(';+', ';').split('+'): - parse_token(token, 'X or K') - for token in c_code.replace('-', '+-').replace(';+', ';').split('+'): - parse_token(token, 'C') - else: - for token in description.replace('-', '+-').replace(';+', ';').split('+'): - # dftd3 cannot be used in a custom xc description - assert '-d3' not in token - parse_token(token, 'compound XC', search_xc_alias=True) - return tuple(hyb), tuple(remove_dup(fn_facs)) - -def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None): - r'''Interface to call libxc library to evaluate XC functional, potential - and functional derivatives. - - * The given functional xc_code must be a one-line string. - * The functional xc_code is case-insensitive. - * The functional xc_code string has two parts, separated by ",". The - first part describes the exchange functional, the second part sets the - correlation functional. - - - If "," not appeared in string, the entire string is treated as the - name of a compound functional (containing both the exchange and - the correlation functional) which was declared in the functional - aliases list. The full list of functional aliases can be obtained by - calling the function pyscf.dft.xcfun.XC_ALIAS.keys() . - - If the string was not found in the aliased functional list, it is - treated as X functional. - - - To input only X functional (without C functional), leave the second - part blank. E.g. description='slater,' means a functional with LDA - contribution only. - - - To neglect the contribution of X functional (just apply C functional), - leave blank in the first part, e.g. description=',vwn' means a - functional with VWN only. - - - If compound XC functional is specified, no matter whether it is in the - X part (the string in front of comma) or the C part (the string behind - comma), both X and C functionals of the compound XC functional will be - used. - - * The functional name can be placed in arbitrary order. Two names need to - be separated by operators "+" or "-". Blank spaces are ignored. - NOTE the parser only reads operators "+" "-" "*". / is not supported. - - * A functional name can have at most one factor. If the factor is not - given, it is set to 1. Compound functional can be scaled as a unit. For - example '0.5*b3lyp' is equivalent to - 'HF*0.1 + .04*LDA + .36*B88, .405*LYP + .095*VWN' - - * String "HF" stands for exact exchange (HF K matrix). "HF" can be put in - the correlation functional part (after comma). Putting "HF" in the - correlation part is the same to putting "HF" in the exchange part. - - * String "RSH" means range-separated operator. Its format is - RSH(omega, alpha, beta). Another way to input RSH is to use keywords - SR_HF and LR_HF: "SR_HF(0.1) * alpha_plus_beta" and "LR_HF(0.1) * - alpha" where the number in parenthesis is the value of omega. - - * Be careful with the libxc convention of GGA functional, in which the LDA - contribution is included. - - Args: - xc_code : str - A string to describe the linear combination of different XC functionals. - The X and C functional are separated by comma like '.8*LDA+.2*B86,VWN'. - If "HF" (exact exchange) is appeared in the string, the HF part will - be skipped. If an empty string "" is given, the returns exc, vxc,... - will be vectors of zeros. - rho : ndarray - Shape of ((*,N)) for electron density (and derivatives) if spin = 0; - Shape of ((*,N),(*,N)) for alpha/beta electron density (and derivatives) if spin > 0; - where N is number of grids. - rho (*,N) are ordered as (den,grad_x,grad_y,grad_z,laplacian,tau) - where grad_x = d/dx den, laplacian = \nabla^2 den, tau = 1/2(\nabla f)^2 - In spin unrestricted case, - rho is ((den_u,grad_xu,grad_yu,grad_zu,laplacian_u,tau_u) - (den_d,grad_xd,grad_yd,grad_zd,laplacian_d,tau_d)) - - Kwargs: - spin : int - spin polarized if spin > 0 - relativity : int - No effects. - verbose : int or object of :class:`Logger` - No effects. - - Returns: - ex, vxc, fxc, kxc - - where - - * vxc = (vrho, vsigma, vlapl, vtau) for restricted case - - * vxc for unrestricted case - | vrho[:,2] = (u, d) - | vsigma[:,3] = (uu, ud, dd) - | vlapl[:,2] = (u, d) - | vtau[:,2] = (u, d) - - * fxc for restricted case: - (v2rho2, v2rhosigma, v2sigma2, v2lapl2, v2tau2, v2rholapl, v2rhotau, v2lapltau, v2sigmalapl, v2sigmatau) - - * fxc for unrestricted case: - | v2rho2[:,3] = (u_u, u_d, d_d) - | v2rhosigma[:,6] = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd) - | v2sigma2[:,6] = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd) - | v2lapl2[:,3] - | v2tau2[:,3] = (u_u, u_d, d_d) - | v2rholapl[:,4] - | v2rhotau[:,4] = (u_u, u_d, d_u, d_d) - | v2lapltau[:,4] - | v2sigmalapl[:,6] - | v2sigmatau[:,6] = (uu_u, uu_d, ud_u, ud_d, dd_u, dd_d) - - * kxc for restricted case: - (v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, - v3rho2lapl, v3rho2tau, - v3rhosigmalapl, v3rhosigmatau, - v3rholapl2, v3rholapltau, v3rhotau2, - v3sigma2lapl, v3sigma2tau, - v3sigmalapl2, v3sigmalapltau, v3sigmatau2, - v3lapl3, v3lapl2tau, v3lapltau2, v3tau3) - - * kxc for unrestricted case: - | v3rho3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d) - | v3rho2sigma[:,9] = (u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd) - | v3rhosigma2[:,12] = (u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd) - | v3sigma3[:,10] = (uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd) - | v3rho2lapl[:,6] - | v3rho2tau[:,6] = (u_u_u, u_u_d, u_d_u, u_d_d, d_d_u, d_d_d) - | v3rhosigmalapl[:,12] - | v3rhosigmatau[:,12] = (u_uu_u, u_uu_d, u_ud_u, u_ud_d, u_dd_u, u_dd_d, - d_uu_u, d_uu_d, d_ud_u, d_ud_d, d_dd_u, d_dd_d) - | v3rholapl2[:,6] - | v3rholapltau[:,8] - | v3rhotau2[:,6] = (u_u_u, u_u_d, u_d_d, d_u_u, d_u_d, d_d_d) - | v3sigma2lapl[:,12] - | v3sigma2tau[:,12] = (uu_uu_u, uu_uu_d, uu_ud_u, uu_ud_d, uu_dd_u, uu_dd_d, - ud_ud_u, ud_ud_d, ud_dd_u, ud_dd_d, dd_dd_u, dd_dd_d) - | v3sigmalapl2[:,9] - | v3sigmalapltau[:,12] - | v3sigmatau2[:,9] = (uu_u_u, uu_u_d, uu_d_d, ud_u_u, ud_u_d, ud_d_d, dd_u_u, dd_u_d, dd_d_d) - | v3lapl3[:,4] - | v3lapl2tau[:,6] - | v3lapltau2[:,6] - | v3tau3[:,4] = (u_u_u, u_u_d, u_d_d, d_d_d) - - see also libxc_itrf.c - ''' # noqa: E501 - outbuf = _eval_xc(xc_code, rho, spin, deriv, omega) - exc = outbuf[0] - vxc = fxc = kxc = None - xctype = xc_type(xc_code) - if xctype == 'LDA' and spin == 0: - if deriv > 0: - vxc = [outbuf[1]] - if deriv > 1: - fxc = [outbuf[2]] - if deriv > 2: - kxc = [outbuf[3]] - elif xctype == 'GGA' and spin == 0: - if deriv > 0: - vxc = [outbuf[1], outbuf[2]] - if deriv > 1: - fxc = [outbuf[3], outbuf[4], outbuf[5]] - if deriv > 2: - kxc = [outbuf[6], outbuf[7], outbuf[8], outbuf[9]] - elif xctype == 'LDA' and spin == 1: - if deriv > 0: - vxc = [outbuf[1:3].T] - if deriv > 1: - fxc = [outbuf[3:6].T] - if deriv > 2: - kxc = [outbuf[6:10].T] - elif xctype == 'GGA' and spin == 1: - if deriv > 0: - vxc = [outbuf[1:3].T, outbuf[3:6].T] - if deriv > 1: - fxc = [outbuf[6:9].T, outbuf[9:15].T, outbuf[15:21].T] - if deriv > 2: - kxc = [outbuf[21:25].T, outbuf[25:34].T, outbuf[34:46].T, outbuf[46:56].T] - elif xctype == 'MGGA' and spin == 0: - if deriv > 0: - vxc = [outbuf[1], outbuf[2], None, outbuf[3]] - if deriv > 1: - fxc = [ - # v2rho2, v2rhosigma, v2sigma2, - outbuf[4], outbuf[5], outbuf[6], - # v2lapl2, v2tau2, - None, outbuf[9], - # v2rholapl, v2rhotau, - None, outbuf[7], - # v2lapltau, v2sigmalapl, v2sigmatau, - None, None, outbuf[8]] - if deriv > 2: - # v3lapltau2 might not be strictly 0 - # outbuf[18] = 0 - kxc = [ - # v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, - outbuf[10], outbuf[11], outbuf[12], outbuf[13], - # v3rho2lapl, v3rho2tau, - None, outbuf[14], - # v3rhosigmalapl, v3rhosigmatau, - None, outbuf[15], - # v3rholapl2, v3rholapltau, v3rhotau2, - None, None, outbuf[16], - # v3sigma2lapl, v3sigma2tau, - None, outbuf[17], - # v3sigmalapl2, v3sigmalapltau, v3sigmatau2, - None, None, outbuf[18], - # v3lapl3, v3lapl2tau, v3lapltau2, v3tau3) - None, None, None, outbuf[19]] - elif xctype == 'MGGA' and spin == 1: - if deriv > 0: - vxc = [outbuf[1:3].T, outbuf[3:6].T, None, outbuf[6:8].T] - if deriv > 1: - # v2lapltau might not be strictly 0 - # outbuf[39:43] = 0 - fxc = [ - # v2rho2, v2rhosigma, v2sigma2, - outbuf[8:11].T, outbuf[11:17].T, outbuf[17:23].T, - # v2lapl2, v2tau2, - None, outbuf[33:36].T, - # v2rholapl, v2rhotau, - None, outbuf[23:27].T, - # v2lapltau, v2sigmalapl, v2sigmatau, - None, None, outbuf[27:33].T] - if deriv > 2: - # v3lapltau2 might not be strictly 0 - # outbuf[204:216] = 0 - kxc = [ - # v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, - outbuf[36:40].T, outbuf[40:49].T, outbuf[49:61].T, outbuf[61:71].T, - # v3rho2lapl, v3rho2tau, - None, outbuf[71:77].T, - # v3rhosigmalapl, v3rhosigmatau, - None, outbuf[77:89].T, - # v3rholapl2, v3rholapltau, v3rhotau2, - None, None, outbuf[89:95].T, - # v3sigma2lapl, v3sigma2tau, - None, outbuf[95:107].T, - # v3sigmalapl2, v3sigmalapltau, v3sigmatau2, - None, None, outbuf[107:116].T, - # v3lapl3, v3lapl2tau, v3lapltau2, v3tau3) - None, None, None, outbuf[116:120].T] - return exc, vxc, fxc, kxc - -_GGA_SORT = { - (1, 2): numpy.array([ - 6, 7, 9, 10, 11, 8, 12, 13, 14, 15, 16, 17, 18, 19, 20, - ]), - (1, 3): numpy.array([ - 21, 22, 25, 26, 27, 23, 28, 29, 30, 34, 35, 36, 37, 38, 39, 24, 31, 32, - 33, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, - ]), - (1, 4): numpy.array([ - 56, 57, 61, 62, 63, 58, 64, 65, 66, 73, 74, 75, 76, 77, 78, 59, 67, 68, - 69, 79, 80, 81, 82, 83, 84, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 60, - 70, 71, 72, 85, 86, 87, 88, 89, 90, 101, 102, 103, 104, 105, 106, 107, - 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, - 122, 123, 124, 125, - ]) -} -_MGGA_SORT = { - (0, 2): numpy.array([ - 0, # v2rho2 - 1, # v2rhosigma - 3, # v2rhotau - 2, # v2sigma2 - 4, # v2sigmatau - 5, # v2tau2 - ]) + 4, - (0, 3): numpy.array([ - 0, # v3rho3 - 1, # v3rho2sigma - 4, # v3rho2tau - 2, # v3rhosigma2 - 5, # v3rhosigmatau - 6, # v3rhotau2 - 3, # v3sigma3 - 7, # v3sigma2tau - 8, # v3sigmatau2 - 9, # v3tau3 - ]) + 10, - (0, 4): numpy.array([ - 0, # v4rho4 - 1, # v4rho3sigma - 5, # v4rho3tau - 2, # v4rho2sigma2 - 6, # v4rho2sigmatau - 7, # v4rho2tau2 - 3, # v4rhosigma3 - 8, # v4rhosigma2tau - 9, # v4rhosigmatau2 - 10, # v4rhotau3 - 4, # v4sigma4 - 11, # v4sigma3tau - 12, # v4sigma2tau2 - 13, # v4sigmatau3 - 14, # v4tau4 - ]) + 20, - (1, 2): numpy.array([ - 8, 9, 11, 12, 13, 23, 24, 10, 14, 15, 16, 25, 26, 17, 18, 19, 27, 28, - 20, 21, 29, 30, 22, 31, 32, 33, 34, 35, - ]), - (1, 3): numpy.array([ - 36, 37, 40, 41, 42, 71, 72, 38, 43, 44, 45, 73, 74, 49, 50, 51, 77, 78, - 52, 53, 79, 80, 54, 81, 82, 89, 90, 91, 39, 46, 47, 48, 75, 76, 55, 56, - 57, 83, 84, 58, 59, 85, 86, 60, 87, 88, 92, 93, 94, 61, 62, 63, 95, 96, - 64, 65, 97, 98, 66, 99, 100, 107, 108, 109, 67, 68, 101, 102, 69, 103, - 104, 110, 111, 112, 70, 105, 106, 113, 114, 115, 116, 117, 118, 119, - ]), - (1, 4): numpy.array([ - 120, 121, 125, 126, 127, 190, 191, 122, 128, 129, 130, 192, 193, 137, - 138, 139, 198, 199, 140, 141, 200, 201, 142, 202, 203, 216, 217, 218, - 123, 131, 132, 133, 194, 195, 143, 144, 145, 204, 205, 146, 147, 206, - 207, 148, 208, 209, 219, 220, 221, 155, 156, 157, 225, 226, 158, 159, - 227, 228, 160, 229, 230, 249, 250, 251, 161, 162, 231, 232, 163, 233, - 234, 252, 253, 254, 164, 235, 236, 255, 256, 257, 267, 268, 269, 270, - 124, 134, 135, 136, 196, 197, 149, 150, 151, 210, 211, 152, 153, 212, - 213, 154, 214, 215, 222, 223, 224, 165, 166, 167, 237, 238, 168, 169, - 239, 240, 170, 241, 242, 258, 259, 260, 171, 172, 243, 244, 173, 245, - 246, 261, 262, 263, 174, 247, 248, 264, 265, 266, 271, 272, 273, 274, - 175, 176, 177, 275, 276, 178, 179, 277, 278, 180, 279, 280, 295, 296, - 297, 181, 182, 281, 282, 183, 283, 284, 298, 299, 300, 184, 285, 286, - 301, 302, 303, 313, 314, 315, 316, 185, 186, 287, 288, 187, 289, 290, - 304, 305, 306, 188, 291, 292, 307, 308, 309, 317, 318, 319, 320, 189, - 293, 294, 310, 311, 312, 321, 322, 323, 324, 325, 326, 327, 328, 329, - ]) -} - -def eval_xc1(xc_code, rho, spin=0, deriv=1, omega=None): - '''Similar to eval_xc. - Returns an array with the order of derivatives following xcfun convention. - ''' - out = _eval_xc(xc_code, rho, spin, deriv=deriv, omega=omega) - xctype = xc_type(xc_code) - idx = _libxc_to_xcfun_indices(xctype, spin, deriv) - return out[idx] - -def _libxc_to_xcfun_indices(xctype, spin=0, deriv=1): - if deriv <= 1: - return slice(None) - elif xctype == 'LDA' or xctype == 'HF': - return slice(None) - elif xctype == 'GGA': - if spin == 0: - return slice(None) - else: - idx = [numpy.arange(6)] # up to deriv=1 - for i in range(2, deriv+1): - idx.append(_GGA_SORT[(spin, i)]) - else: # MGGA - if spin == 0: - idx = [numpy.arange(4)] # up to deriv=1 - else: - idx = [numpy.arange(8)] # up to deriv=1 - for i in range(2, deriv+1): - idx.append(_MGGA_SORT[(spin, i)]) - return numpy.hstack(idx) - -def _eval_xc(xc_code, rho, spin=0, deriv=1, omega=None): - xc = _get_xc(xc_code, spin, omega) - assert deriv <= xc.max_deriv_order - xctype = xc.xc_type - assert xctype in ('HF', 'LDA', 'GGA', 'MGGA') - - rho = numpy.asarray(rho, order='C', dtype=numpy.double) - if xctype == 'MGGA' and rho.shape[-2] == 6: - rho = numpy.asarray(rho[...,[0,1,2,3,5],:], order='C') - - if xc.needs_laplacian: - raise NotImplementedError('laplacian in meta-GGA method') - - nvar, xlen = xc_deriv._XC_NVAR[xctype, spin] - ngrids = rho.shape[-1] - rho = rho.reshape(spin+1,nvar,ngrids) - outlen = lib.comb(xlen+deriv, deriv) - out = numpy.zeros((outlen,ngrids)) - n = xc.nfunc - if n > 0: - _itrf.LIBXC_eval_xc(n, - xc.xc_arr, - (ctypes.c_double*n)(*xc.facs), - spin, deriv, - nvar, ngrids, - outlen, - rho.ctypes, - out.ctypes) - return out - -def eval_xc_eff(xc_code, rho, deriv=1, omega=None): - r'''Returns the derivative tensor against the density parameters - - [density_a, (nabla_x)_a, (nabla_y)_a, (nabla_z)_a, tau_a] - - or spin-polarized density parameters - - [[density_a, (nabla_x)_a, (nabla_y)_a, (nabla_z)_a, tau_a], - [density_b, (nabla_x)_b, (nabla_y)_b, (nabla_z)_b, tau_b]]. - - It differs from the eval_xc method in the derivatives of non-local part. - The eval_xc method returns the XC functional derivatives to sigma - (|\nabla \rho|^2) - - Args: - rho: 2-dimensional or 3-dimensional array - Total density or (spin-up, spin-down) densities (and their - derivatives if GGA or MGGA functionals) on grids - - Kwargs: - deriv: int - derivative orders - omega: float - define the exponent in the attenuated Coulomb for RSH functional - ''' - xctype = xc_type(xc_code) - rho = numpy.asarray(rho, order='C', dtype=numpy.double) - if xctype == 'MGGA' and rho.shape[-2] == 6: - rho = numpy.asarray(rho[...,[0,1,2,3,5],:], order='C') - - spin_polarized = rho.ndim >= 2 and rho.shape[0] == 2 - if spin_polarized: - spin = 1 - else: - spin = 0 - out = eval_xc1(xc_code, rho, spin, deriv, omega) - return xc_deriv.transform_xc(rho, out, xctype, spin, deriv) - -def define_xc_(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)): - '''Define XC functional. See also :func:`eval_xc` for the rules of input description. - - Args: - ni : an instance of :class:`NumInt` - - description : str - A string to describe the linear combination of different XC functionals. - The X and C functional are separated by comma like '.8*LDA+.2*B86,VWN'. - If "HF" was appeared in the string, it stands for the exact exchange. - - Kwargs: - xctype : str - 'LDA' or 'GGA' or 'MGGA' - hyb : float - hybrid functional coefficient - rsh : a list of three floats - coefficients (omega, alpha, beta) for range-separated hybrid functional. - omega is the exponent factor in attenuated Coulomb operator e^{-omega r_{12}}/r_{12} - alpha is the coefficient for long-range part, hybrid coefficient - can be obtained by alpha + beta - - Examples: - - >>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz') - >>> mf = dft.RKS(mol) - >>> define_xc_(mf._numint, '.2*HF + .08*LDA + .72*B88, .81*LYP + .19*VWN') - >>> mf.kernel() - -76.3783361189611 - >>> define_xc_(mf._numint, 'LDA*.08 + .72*B88 + .2*HF, .81*LYP + .19*VWN') - >>> mf.kernel() - -76.3783361189611 - >>> def eval_xc(xc_code, rho, *args, **kwargs): - ... exc = 0.01 * rho**2 - ... vrho = 0.01 * 3 * rho**2 - ... vxc = (vrho, None, None, None) - ... fxc = None # 2nd order functional derivative - ... kxc = None # 3rd order functional derivative - ... return exc, vxc, fxc, kxc - >>> define_xc_(mf._numint, eval_xc, xctype='LDA') - >>> mf.kernel() - 48.8525211046668 - ''' - if isinstance(description, str): - def _eval_xc(xc_code, rho, *args, **kwargs): - return eval_xc(description, rho, *args, **kwargs) - ni.eval_xc = _eval_xc - ni.hybrid_coeff = lambda *args, **kwargs: hybrid_coeff(description) - ni.rsh_coeff = lambda *args: rsh_coeff(description) - ni._xc_type = lambda *args: xc_type(description) - - elif callable(description): - ni.eval_xc = _eval_xc = description - ni.hybrid_coeff = lambda *args, **kwargs: hyb - ni.rsh_coeff = lambda *args, **kwargs: rsh - ni._xc_type = lambda *args: xctype - - else: - raise ValueError('Unknown description %s' % description) - - def _eval_xc1(xc_code, rho, spin=0, deriv=1, omega=None): - libxc_out = _eval_xc(xc_code, rho, spin, deriv=deriv, omega=omega) - nvar, xlen = xc_deriv._XC_NVAR[xctype, spin] - outlen = lib.comb(xlen+deriv, deriv) - exc, vxc, fxc, kxc = libxc_out[:4] - out = [exc] - if deriv > 0: - assert vxc is not None - out.extend([x for x in vxc if x is not None]) - if deriv > 1: - assert fxc is not None - if xctype == 'GGA': - assert len(fxc) == 3, 'fxc for GGA should be arranged as (v2rho2, v2rhosigma, v2sigma2)' - elif xctype == 'MGGA': - if len(fxc) == 10: - fxc = [fxc[i] for i in [0, 1, 2, 6, 4, 9]] - else: - assert len(fxc) == 6, ( - 'fxc for MGGA should be arranged as\n' - '(v2rho2, v2rhosigma, v2sigma2, v2tau2, v2rhotau, v2sigmatau)\nor\n' - '(v2rho2, v2rhosigma, v2sigma2, v2lapl2, v2tau2, ' - 'v2rholapl, v2rhotau, v2lapltau, v2sigmalapl, v2sigmatau)') - assert all(x is not None for x in fxc) - out.extend(fxc) - if deriv > 2: - assert kxc is not None - out.extend([x for x in kxc if x is not None]) - if spin == 1: - # Returns of eval_xc are structured as [grid_id,deriv_component] - # for each term in libxc_out. Change the shape to [deriv_comp, grid_id] - out = [x.T for x in out] - out = numpy.vstack(out)[:outlen] - assert len(out) == outlen - idx = _libxc_to_xcfun_indices(xctype, spin, deriv) - return out[idx] - ni.eval_xc1 = _eval_xc1 - return ni - -def define_xc(ni, description, xctype='LDA', hyb=0, rsh=(0,0,0)): - return define_xc_(ni.copy(), description, xctype, hyb, rsh) -define_xc.__doc__ = define_xc_.__doc__ - -class XCFunctionalCache: - ''' - Data structure for caching LibXC C data structures - - Attributes for XCFunctionalCache: - xc_code : str - A string to describe the XC functional. See `parse_xc` for details. - spin : int - spin polarized if spin > 0 - xc_arr : ctypes array containing all `xc_func_type` struct of each - functional component - xc_objs : list of ctypes pointers of the `xc_func_type` struct of each - functional component. The objects point to the same memory address - in `xc_arr` - nfunc : int - number of LibXC functional components - hyb : hybrid coefficients as produced by `parse_xc()` - fn_ids : tuple of int - LibXC ID of each functional component - facs : tuple of float - factor of each functional component - - ''' - def __init__(self, xc_code, spin=0, omega=None, density_threshold=None): - if isinstance(xc_code, str) and '__VV10' in xc_code: - raise RuntimeError('Deprecated notation for NLC functional.') - self.xc_code = xc_code - self.spin = spin - self.hyb, fn_facs = parse_xc(xc_code) - if fn_facs: - self.fn_ids, self.facs = zip(*fn_facs) - else: - self.fn_ids = tuple() - self.facs = tuple() - self.nfunc = len(self.facs) - if omega is not None: - self.hyb = self.hyb[:2] + (float(omega),) - if self.hyb[2] != 0: - # Current implementation does not support different omegas for - # different RSH functionals if there are multiple RSHs - omega = [self.hyb[2]] * len(self.facs) - else: - omega = [0] * len(self.facs) - fn_ids_ctypes = (ctypes.c_int * self.nfunc)(*self.fn_ids) - problem_xc = set(self.fn_ids).intersection(PROBLEMATIC_XC) - if problem_xc: - problem_xc = [PROBLEMATIC_XC[k] for k in problem_xc] - warnings.warn('Libxc functionals %s may have discrepancy to xcfun ' - 'library.\n' % problem_xc) - - self.xc_arr = _itrf.LIBXC_xc_func_init(self.nfunc, fn_ids_ctypes, spin) - - if self.xc_arr is None: - del self.xc_arr - raise ValueError(f"{xc_code} is not a valid functional") - if density_threshold or any(omega): - self._set_omega_density_threshold_(omega, density_threshold) - self.xc_objs = [ctypes.cast(self.xc_arr + _XC_FUNC_TYPE_SIZE * i, ctypes.c_void_p) - for i in range(self.nfunc)] - - @cached_property - def xc_type(self): - return _xc_type(self.nfunc, self.xc_arr) - - @cached_property - def is_nlc(self): - return _is_nlc(self.nfunc, self.xc_arr) - - @cached_property - def needs_laplacian(self): - return _needs_laplacian(self.nfunc, self.xc_arr) - - @cached_property - def max_deriv_order(self): - return _max_deriv_order(self.nfunc, self.xc_arr) - - @cached_property - def hybrid_coeff(self): - return _hybrid_coeff(self.xc_objs, self.hyb, self.facs) - - @cached_property - def nlc_coeff(self): - return _nlc_coeff(self.xc_objs, self.hyb, self.facs) - - @cached_property - def rsh_coeff(self): - return _rsh_coeff(self.xc_objs, self.hyb, self.facs, self.xc_code) - - def customize_(self, ext_params=None, omega=None, hyb=None, facs=None, - density_threshold=None, callback=None, clear=True): - ''' - Customize the functional in place. See `register_custom_functional_` for - a detailed description of the arguments. - - Additional argument: - clear: bool - If True, clear all cached properties. Set to False may have a slight performace - gain, but one must manually clear cached properties that might change prior to - any productive calculations, preferably in the callback. - ''' - obj_by_id = self.obj_by_id() - if density_threshold or (omega and any(omega)): - self._set_omega_density_threshold_(omega, density_threshold) - if hyb is not None: - self.hyb = hyb - if facs is not None: - self.facs = facs - if ext_params is not None: - for xid, param in ext_params.items(): - param = numpy.asarray(param, dtype=numpy.double) - func = obj_by_id[xid] - info = _itrf.xc_func_get_info(func) - n = _itrf.xc_func_info_get_n_ext_params(info) - assert param.size == n, \ - f"""Unexpected size of external parameters for functional {xid}. -Expected {n} but {param.size} provided.""" - _itrf.xc_func_set_ext_params(func, param.ctypes.data_as(ctypes.POINTER(ctypes.c_double))) - if callable(callback): - callback(self, obj_by_id, self.spin) - - if clear: - # clear all cached properties - self._clear_cached_property('xc_type') - self._clear_cached_property('is_nlc') - self._clear_cached_property('needs_laplacian') - self._clear_cached_property('max_deriv_order') - self._clear_cached_property('hybrid_coeff') - self._clear_cached_property('nlc_coeff') - self._clear_cached_property('rsh_coeff') - - def _clear_cached_property(self, key): - try: - delattr(self, key) - except AttributeError: - pass - - def _set_omega_density_threshold_(self, omega, density_threshold): - if density_threshold is None: - density_threshold = 0. - _itrf.LIBXC_xc_func_set_params(self.nfunc, self.xc_arr, - (ctypes.c_double * self.nfunc)(*omega), - density_threshold) - - def __del__(self): - try: - _itrf.LIBXC_xc_func_end(self.nfunc, self.xc_arr) - except AttributeError: - # skip destructor call when XC failed to initialize - pass - - def obj_by_id(self): - ''' - Returns a dict where LibXC functional IDs are used as the key, - and the corresponding ctypes void pointer of the `xc_func_type` - structs are used as the value. This dict is useful for users who - wish to perform low-level operations on the underlying LibXC data - structure. - ''' - return {xid: func for xid, func in zip(self.fn_ids, self.xc_objs)} - -_CUSTOM_FUNC_R = {} -_CUSTOM_FUNC_U = {} - -@lru_cache(100) -def _get_system_xc(xc_code, spin, omega=None): - return XCFunctionalCache(xc_code, spin, omega) - -def _get_xc(xc_code, spin=0, omega=None): - if omega is not None: - return _get_system_xc(xc_code, spin, omega) - try: - if spin > 0: - return _CUSTOM_FUNC_U[xc_code] - else: - return _CUSTOM_FUNC_R[xc_code] - except KeyError: - return _get_system_xc(xc_code, spin) - -def register_custom_functional_(new_xc_code, based_on_xc_code, ext_params=None, omega=None, - hyb=None, facs=None, density_threshold=None, callback=None, - restricted=True, unrestricted=True): - ''' - Register a custom functional with name `new_xc_code` based on an existing LibXC - functional with xc_code `based_on_xc_code`. When this function is called, the - functional with xc_code `based_on_xc_code` will be constructed first. Then users can - perform customization to the functional by setting the respective keyword arguments. - Common customization include setting functional external parameters, adjust mixing - factors of each functional components, and adjusting hybrid coefficients. Users can - also use a callback function to perform more advanced modification to the LibXC data - structure using ctypes or modify the parameters in XCFunctionalCache. Once registered, - users may use the custom functional by passing `new_xc_code` to any PySCF module that - uses LibXC. If `new_xc_code` has previously been registered, the older functional will - be replaced with the new definition. In this case, however, users are encouraged to use - `update_custom_functional_` instead to avoid reconstructing the LibXC C data structure. - - Args: - new_xc_code : str - The string identifier of the customized functional - based_on_xc_code : str - The string identifier of the parent functional - ext_params : dict, with LibXC functional integer ID as key, and an array-like - object containing the functional parameters as value. - If not None, set the external parameters of the functional componet from - the dict using the xc_func_set_ext_params API. - omega : float - If not None, set the omega value in self.hyb. If hyb argument is set, - this argument is ignored. - hyb : tuple of float - If not None, set the hybrid coefficients self.hyb - facs : tuple of float - If not None, set the factors of each functional component self.facs - callback : callable - Callback function to perform more advanced customization. - The function takes 3 parameters: - xc: a XCFunctionalCache instance of the functional - obj_by_id: a dict constructed through xc.obj_by_id, used to access - underlying LibXC C data structure - spin: 0 for restricted and 1 for unrestricted - The callback function will make inplace modifications to xc object or - the underlying LibXC C data structure. Depending on the `restricted` - and `unrestricted` arguments, the callback may be called more than once. - restricted : bool - If True, create the spin-restricted version of the custom functional - unrestricted : bool - If True, create the spin-unrestricted version of the custom functional - ''' - def build(spin): - xc = XCFunctionalCache(based_on_xc_code, spin) - xc.xc_code = new_xc_code - xc.customize_(ext_params, omega, hyb, facs, density_threshold, callback, False) - return xc - - try: - parse_xc(new_xc_code) - warnings.warn('Attempting to register a custom functional with xc_code ' - f'"{new_xc_code}", which is a valid built-in xc_code. ' - 'Consider using a different xc_code to avoid confusion.') - except Exception: - pass - - if restricted: - _CUSTOM_FUNC_R[new_xc_code] = build(0) - if unrestricted: - _CUSTOM_FUNC_U[new_xc_code] = build(1) - -def update_custom_functional_(xc_code, ext_params=None, omega=None, - hyb=None, facs=None, density_threshold=None, callback=None, - restricted=True, unrestricted=True, clear=True): - ''' - Update the custom functional with name `xc_code` that was previously registered - through `register_custom_functional_`. See `register_custom_functional_` for - a detailed description of the keyword arguments. - - Additional keyword argument: - clear: bool - If True, clear all cached properties. Clearing cached properties are - required if one or more cached properties (such as hybrid_coeff) are - changed after the update. Set to False allows users to manually handle - cache; one can clear only the properties that are expected to change - for a slight performace gain. - ''' - def update(xc): - xc.customize_(ext_params, omega, hyb, facs, density_threshold, callback, clear) - - if restricted: - update(_CUSTOM_FUNC_R[xc_code]) - if unrestricted: - update(_CUSTOM_FUNC_U[xc_code]) - -def unregister_custom_functional_(xc_code): - ''' - Unregister a custom functional with name `xc_code` that was previously registered - through `register_custom_functional_`. - ''' - try: - del _CUSTOM_FUNC_R[xc_code] - except KeyError: - pass - - try: - del _CUSTOM_FUNC_U[xc_code] - except KeyError: - pass - diff --git a/pyscf/dft2/test/dm_h4.npy b/pyscf/dft2/test/dm_h4.npy deleted file mode 100644 index 3581a8c8e..000000000 Binary files a/pyscf/dft2/test/dm_h4.npy and /dev/null differ diff --git a/pyscf/dft2/test/test_grad_metagga_mcpdft.py b/pyscf/dft2/test/test_grad_metagga_mcpdft.py deleted file mode 100644 index 6406f3c64..000000000 --- a/pyscf/dft2/test/test_grad_metagga_mcpdft.py +++ /dev/null @@ -1,146 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2025 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Author: Matthew Hennefarth - -import unittest - -from pyscf import scf, gto, df, dft -from pyscf.data.nist import BOHR -from pyscf import mcpdft - -def diatomic( - atom1, - atom2, - r, - fnal, - basis, - ncas, - nelecas, - nstates, - charge=None, - spin=None, - symmetry=False, - cas_irrep=None, - density_fit=False, - grids_level=9, -): - """Used for checking diatomic systems to see if the Lagrange Multipliers are working properly.""" - global mols - xyz = "{:s} 0.0 0.0 0.0; {:s} {:.3f} 0.0 0.0".format(atom1, atom2, r) - mol = gto.M( - atom=xyz, - basis=basis, - charge=charge, - spin=spin, - symmetry=symmetry, - verbose=0, - output="/dev/null", - ) - mols.append(mol) - mf = scf.RHF(mol) - if density_fit: - mf = mf.density_fit(auxbasis=df.aug_etb(mol)) - - mc = mcpdft.CASSCF(mf.run(), fnal, ncas, nelecas, grids_level=grids_level) - if spin is None: - spin = mol.nelectron % 2 - - ss = spin * (spin + 2) * 0.25 - mc.fix_spin_(ss=ss, shift=2) - - if nstates > 1: - mc = mc.state_average( - [ - 1.0 / float(nstates), - ] - * nstates, - ) - - mc.conv_tol = 1e-12 - mc.conv_grad_tol = 1e-6 - mo = None - if symmetry and (cas_irrep is not None): - mo = mc.sort_mo_by_irrep(cas_irrep) - - mc_grad = mc.run(mo).nuc_grad_method() - mc_grad.conv_rtol = 1e-12 - return mc_grad - - -def setUpModule(): - global mols, original_grids - from importlib import reload - from pyscf import dft2 - dft.libxc = dft2.libxc - reload (mcpdft) - reload (mcpdft.otfnal) - mols = [] - original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False - - -def tearDownModule(): - global mols, diatomic, original_grids - [m.stdout.close() for m in mols] - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids - del mols, diatomic, original_grids - - -class KnownValues(unittest.TestCase): - - def test_grad_lih_ssmc2322_sto3g(self): - mc = diatomic("Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 1, grids_level=1) - de = mc.kernel()[1, 0] / BOHR - - # Numerical from this software - # PySCF commit: f2c2d3f963916fb64ae77241f1b44f24fa484d96 - # PySCF-forge commit: e82ba940654cd0b91f799e889136a316fda34b10 - DE_REF = -1.0641645070 - - self.assertAlmostEqual(de, DE_REF, 5) - - def test_grad_lih_sa2mc2322_sto3g(self): - mc = diatomic("Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 2, grids_level=1) - - # Numerical from this software - # PySCF commit: f2c2d3f963916fb64ae77241f1b44f24fa484d96 - # PySCF-forge commit: e82ba940654cd0b91f799e889136a316fda34b10 - DE_REF = [-1.0510225010, -0.8963063432] - - for state in range(2): - with self.subTest(state=state): - de = mc.kernel(state=state)[1, 0] / BOHR - self.assertAlmostEqual(de, DE_REF[state], 5) - - def test_grad_lih_sa2mc2322_sto3g_df(self): - mc = diatomic( - "Li", "H", 0.8, "MC23", "STO-3G", 2, 2, 2, grids_level=1, density_fit=df - ) - - # Numerical from this software - # PySCF commit: f2c2d3f963916fb64ae77241f1b44f24fa484d96 - # PySCF-forge commit: ee6ac742fbc79d170bc4b63ef2b2c4b49478c53a - DE_REF = [-1.0510303416, -0.8963992331] - - for state in range(2): - with self.subTest(state=state): - de = mc.kernel(state=state)[1, 0] / BOHR - self.assertAlmostEqual(de, DE_REF[state], 5) - - -if __name__ == "__main__": - print("Full Tests for MC-PDFT gradients with meta-GGA functionals") - unittest.main() diff --git a/pyscf/dft2/test/test_libxc.py b/pyscf/dft2/test/test_libxc.py deleted file mode 100644 index 42bef7ec6..000000000 --- a/pyscf/dft2/test/test_libxc.py +++ /dev/null @@ -1,501 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2019 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -# This file is adapted from `dft/test/test_libxc.py` of the PySCF core module -# commit 25eaa9572977b903de24d5c11ad345cecd744728 - -import os -import ctypes -import unittest -import numpy -from pyscf import gto, scf -from pyscf import dft2 as dft -from pyscf import dft as dft1 -from pyscf import lib - -def setUpModule(): - global mol, mf, ao, rho, dm - mol = gto.Mole() - mol.verbose = 0 - mol.output = None - mol.atom = 'h 0 0 0; h 1 .5 0; h 0 4 1; h 1 0 .2' - mol.basis = 'aug-ccpvdz' - mol.build() - #dm = scf.RHF(mol).run(conv_tol=1e-14).make_rdm1() - dm = numpy.load(os.path.realpath(os.path.join(__file__, '..', 'dm_h4.npy'))) - with lib.temporary_env(dft1.radi, ATOM_SPECIFIC_TREUTLER_GRIDS=False): - mf = dft1.RKS(mol) - mf.grids.atom_grid = {"H": (50, 110)} - mf.prune = None - mf.grids.build(with_non0tab=False) - nao = mol.nao_nr() - ao = dft1.numint.eval_ao(mol, mf.grids.coords, deriv=1) - rho = dft1.numint.eval_rho(mol, ao, dm, xctype='GGA') - -def tearDownModule(): - global mol, mf, ao, rho - del mol, mf, ao, rho - -class KnownValues(unittest.TestCase): - @classmethod - def setUpClass(cls): - cls.original_grids = dft1.radi.ATOM_SPECIFIC_TREUTLER_GRIDS - dft1.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False - - @classmethod - def tearDownClass(cls): - dft1.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = cls.original_grids - - def test_parse_xc(self): - hyb, fn_facs = dft.libxc.parse_xc('.5*HF+.5*B3LYP5,VWN*.5') - self.assertAlmostEqual(hyb[0], .6, 12) - self.assertEqual([x[0] for x in fn_facs], [1,106,131,7]) - self.assertTrue(numpy.allclose([x[1] for x in fn_facs], - (0.04, 0.36, 0.405, 0.595))) - hyb, fn_facs = dft.libxc.parse_xc('HF,') - self.assertEqual(hyb[0], 1) - self.assertEqual(fn_facs, ()) - - hyb, fn_facs = dft.libxc.parse_xc('B88 - SLATER') - self.assertEqual(fn_facs, ((106, 1), (1, -1))) - hyb, fn_facs = dft.libxc.parse_xc('B88 -SLATER*.5') - self.assertEqual(fn_facs, ((106, 1), (1, -0.5))) - - hyb = dft.libxc.hybrid_coeff('0.5*B3LYP\n+0.25*B3LYP') - self.assertAlmostEqual(hyb, .15, 12) - - hyb, fn_facs = dft.libxc.parse_xc('0.6*CAM_B3LYP+0.4*B3P86V5') - self.assertTrue(numpy.allclose(hyb, (.08, .08, 0))) - self.assertTrue(numpy.allclose(fn_facs, - ((433, 0.6), (1, 0.032), (106, 0.288), (132, 0.324), (7, 0.076)))) - rsh = dft.libxc.rsh_coeff('0.6*CAM_B3LYP+0.4*B3P86V5') - rsh1 = dft.libxc.rsh_coeff('CAM_B3LYP') - hyb = dft.libxc.hybrid_coeff('B3P86V5') - self.assertTrue(numpy.allclose(rsh, (0.33, 0.6*rsh1[1]+.4*hyb, 0.6*rsh1[2]))) - - hyb, fn_facs = dft.libxc.parse_xc('0.4*B3P86V5+0.6*CAM_B3LYP') - self.assertTrue(numpy.allclose(hyb, (.08, .08, 0))) - self.assertTrue(numpy.allclose(fn_facs, - ((1, 0.032), (106, 0.288), (132, 0.324), (7, 0.076), (433, 0.6)))) - rsh = dft.libxc.rsh_coeff('0.4*B3P86V5+0.6*CAM_B3LYP') - rsh1 = dft.libxc.rsh_coeff('CAM_B3LYP') - hyb = dft.libxc.hybrid_coeff('B3P86V5') - self.assertTrue(numpy.allclose(rsh, (0.33, 0.6*rsh1[1]+.4*hyb, 0.6*rsh1[2]))) - - hyb, fn_facs = dft.libxc.parse_xc('0.5*SR-HF(0.3) + .8*HF + .22*LR_HF') - self.assertEqual(hyb, (1.3, 1.02, 0.3)) - - hyb, fn_facs = dft.libxc.parse_xc('0.5*SR-HF + .22*LR_HF(0.3) + .8*HF') - self.assertEqual(hyb, (1.3, 1.02, 0.3)) - - hyb, fn_facs = dft.libxc.parse_xc('0.5*SR-HF + .8*HF + .22*LR_HF(0.3)') - self.assertEqual(hyb, (1.3, 1.02, 0.3)) - - hyb, fn_facs = dft.libxc.parse_xc('0.5*RSH(2.04;0.56;0.3) + 0.5*BP86') - self.assertEqual(hyb, (1.3, 1.02, 0.3)) - self.assertEqual(fn_facs, ((106, 0.5), (132, 0.5))) - - hyb, fn_facs = dft.libxc.parse_xc('0.5*RSH(.3, 2.04, 0.56) + 0.5*BP86') - self.assertEqual(hyb, (1.3, 1.02, 0.3)) - self.assertEqual(fn_facs, ((106, 0.5), (132, 0.5))) - - rsh = dft.libxc.rsh_coeff('0.5*HSE06+.001*HF') - self.assertTrue(numpy.allclose(rsh, (0.11, .001, 0.125))) - - rsh = dft.libxc.rsh_coeff('0.5*wb97+0.001*HF') - self.assertTrue(numpy.allclose(rsh, (0.4, 0.501, -.5))) - - self.assertRaises(ValueError, dft.libxc.parse_xc, 'SR_HF(0.3) + LR_HF(.5)') - self.assertRaises(ValueError, dft.libxc.parse_xc, 'LR-HF(0.3) + SR-HF(.5)') - - hyb = dft.libxc.hybrid_coeff('0.5*B3LYP+0.2*HF') - self.assertAlmostEqual(hyb, .3, 12) - - hyb = dft.libxc.hybrid_coeff('M05') - self.assertAlmostEqual(hyb, 0.28, 9) - - hyb, fn_facs = dft.libxc.parse_xc('APBE,') - self.assertEqual(fn_facs, ((184, 1),)) - - hyb, fn_facs = dft.libxc.parse_xc('LDA0') - self.assertEqual(fn_facs, ((177, 1),)) - - #hyb, fn_facs = dft.libxc.parse_xc('TF,') - #self.assertEqual(fn_facs, [(50, 1)]) - - hyb, fn_facs = dft.libxc.parse_xc("9.999e-5*HF,") - self.assertEqual(hyb, (9.999e-5, 9.999e-5, 0)) - - ref = ((1, 1), (7, 1)) - self.assertEqual(dft.libxc.parse_xc_name('LDA,VWN'), (1,7)) - self.assertEqual(dft.libxc.parse_xc(('LDA','VWN'))[1], ref) - self.assertEqual(dft.libxc.parse_xc((1, 7))[1], ref) - self.assertEqual(dft.libxc.parse_xc('1, 7')[1], ref) - self.assertEqual(dft.libxc.parse_xc(7)[1], ((7,1),)) - - self.assertEqual(dft.libxc.parse_xc('M11-L')[1], ((226,1),(75,1))) - self.assertEqual(dft.libxc.parse_xc('M11L' )[1], ((226,1),(75,1))) - self.assertEqual(dft.libxc.parse_xc('M11-L,M11L' )[1], ((226,1),(75,1))) - self.assertEqual(dft.libxc.parse_xc('M11_L,M11-L')[1], ((226,1),(75,1))) - self.assertEqual(dft.libxc.parse_xc('M11L,M11_L' )[1], ((226,1),(75,1))) - - self.assertEqual(dft.libxc.parse_xc('Xpbe,')[1], ((123,1),)) - self.assertEqual(dft.libxc.parse_xc('pbe,' )[1], ((101,1),)) - self.assertEqual(dft.libxc.parse_xc('gga_x_pbe_gaussian' )[1], ((321,1),)) - - - hyb, fn_facs = dft.libxc.parse_xc('PBE*.4+LDA') - self.assertEqual(fn_facs, ((101, 0.4), (130, 0.4), (1, 1))) - self.assertRaises(KeyError, dft.libxc.parse_xc, 'PBE+VWN') - - self.assertTrue (dft.libxc.is_meta_gga('m05')) - self.assertFalse(dft.libxc.is_meta_gga('pbe0')) - self.assertFalse(dft.libxc.is_meta_gga('tf,')) - self.assertFalse(dft.libxc.is_meta_gga('vv10')) - self.assertTrue (dft.libxc.is_gga('PBE0')) - self.assertFalse(dft.libxc.is_gga('m05')) - self.assertFalse(dft.libxc.is_gga('tf,')) - self.assertTrue (dft.libxc.is_lda('tf,')) - self.assertFalse(dft.libxc.is_lda('vv10')) - self.assertTrue (dft.libxc.is_hybrid_xc('m05')) - self.assertTrue (dft.libxc.is_hybrid_xc('pbe0,')) - self.assertTrue (dft.libxc.is_hybrid_xc('m05,')) - self.assertFalse(dft.libxc.is_hybrid_xc('vv10')) - self.assertTrue (dft.libxc.is_hybrid_xc((402,'vv10'))) - self.assertTrue (dft.libxc.is_hybrid_xc(('402','vv10'))) - self.assertTrue (dft.libxc.is_nlc('b97mv')) - self.assertTrue (dft.libxc.is_nlc('lc-vv10')) - self.assertTrue (dft.libxc.is_nlc('scanl-vv10')) - self.assertTrue (dft.libxc.is_nlc('b97mv+pbe')) - self.assertTrue (dft.libxc.is_nlc((402, 'b97mv'))) - self.assertTrue (dft.libxc.is_nlc(('402', 'b97mv'))) - - def test_libxc_cam_beta(self): - rsh_tmp = dft.libxc.rsh_coeff(1) - beta = rsh_tmp[2] - self.assertEqual(beta, 0) - - rsh_tmp = dft.libxc.rsh_coeff(433) - rsh_tmp = dft.libxc.rsh_coeff(1) - beta = rsh_tmp[2] - self.assertEqual(beta, 0) - - dft.libxc.is_hybrid_xc(1) - rsh_tmp = dft.libxc.rsh_coeff(1) - beta = rsh_tmp[2] - self.assertEqual(beta, 0) - - def test_nlc_coeff(self): - self.assertEqual(dft.libxc.nlc_coeff('0.5*vv10'), (((5.9, 0.0093), .5),)) - self.assertEqual(dft.libxc.nlc_coeff('pbe+vv10'), (((5.9, 0.0093), 1),)) - - def test_lda(self): - e,v,f,k = dft.libxc.eval_xc('lda,', rho[0], deriv=3) - self.assertAlmostEqual(numpy.dot(rho[0], e) , -789.1150849798871 , 8) - self.assertAlmostEqual(numpy.dot(rho[0], v[0]), -1052.1534466398498, 8) - self.assertAlmostEqual(numpy.dot(rho[0], f[0]), -1762.3340626646932, 8) - self.assertAlmostEqual(numpy.dot(rho[0], k[0]), 1202284274.6255436 , 3) - - e,v,f,k = dft.libxc.eval_xc('lda,', [rho[0]*.5]*2, spin=1, deriv=3) - self.assertAlmostEqual(numpy.dot(rho[0], e) , -789.1150849798871 , 8) - self.assertAlmostEqual(numpy.dot(rho[0], v[0].T[0]), -1052.1534466398498, 8) - self.assertAlmostEqual(numpy.dot(rho[0], v[0].T[1]), -1052.1534466398498, 8) - self.assertAlmostEqual(numpy.dot(rho[0], f[0].T[0]), -1762.3340626646932*2, 8) - self.assertAlmostEqual(numpy.dot(rho[0], f[0].T[2]), -1762.3340626646932*2, 8) - self.assertAlmostEqual(numpy.dot(rho[0], k[0].T[0]), 1202284274.6255436*4, 2) - self.assertAlmostEqual(numpy.dot(rho[0], k[0].T[3]), 1202284274.6255436*4, 2) - - def test_lyp(self): - e,v,f = dft.libxc.eval_xc(',LYP', rho, deriv=2)[:3] - self.assertAlmostEqual(numpy.dot(rho[0],e), -62.114576182676615, 8) - self.assertAlmostEqual(numpy.dot(rho[0],v[0]),-81.771670866308455, 8) - self.assertAlmostEqual(numpy.dot(rho[0],v[1]), 27.485383255125743, 8) - self.assertAlmostEqual(numpy.dot(rho[0],f[0]), 186.823806251777, 2) - self.assertAlmostEqual(numpy.dot(rho[0],f[1]), -3391.2428894571085, 6) - self.assertAlmostEqual(numpy.dot(rho[0],f[2]), 0, 8) - self.assertAlmostEqual(abs(f[2]).sum(), 0, 3) - - #def test_tpss(self): - # #FIXME: raised numerical error - # rho_a = numpy.array([[3.67808547e-08,-2.02358682e-08, 2.16729780e-07, 2.27036045e-07,-1.47795869e-07,-1.45668997e-09]]).T - # e, v = dft.libxc.eval_xc('tpss,', rho_a, spin=0, deriv=1)[:2] - # rho_b = numpy.array([[4.53272893e-06, 4.18968775e-06,-2.83034672e-06, 2.61832978e-06, 5.63360737e-06, 8.97541777e-07]]).T - # e, v = dft.libxc.eval_xc('tpss,', (rho_a, rho_b), spin=1, deriv=1)[:2] - - #TDOO: enable this test when https://gitlab.com/libxc/libxc/-/issues/561 is solved - @unittest.skip('hse03 and hse06 fxc have large numerical errors in Libxc') - def test_hse06(self): - ni = dft.numint.NumInt() - rho = numpy.array([.235, 1.5e-9, 2e-9, 1e-9])[:,None] - xc = 'hse06' - fxc1 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2] - rho = numpy.array([rho*.5]*2) - fxc2 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2] - fxc2 = (fxc2[0,:,0] + fxc2[0,:,1])/2 - self.assertAlmostEqual(abs(fxc2 - fxc1).max(), 0, 12) - - rho = numpy.array([.235, 0, 0, 0])[:,None] - xc = 'hse06' - fxc1 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2] - rho = numpy.array([rho*.5]*2) - fxc2 = ni.eval_xc_eff(xc, rho, deriv=2, xctype='GGA')[2] - fxc2 = (fxc2[0,:,0] + fxc2[0,:,1])/2 - self.assertAlmostEqual(abs(fxc2 - fxc1).max(), 0, 12) - - def test_define_xc_gga(self): - def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None): - # A fictitious XC functional to demonstrate the usage - rho0, dx, dy, dz = rho[:4] - gamma = (dx**2 + dy**2 + dz**2) - exc = .01 * rho0**2 + .02 * (gamma+.001)**.5 - vrho = .01 * 2 * rho0 - vgamma = .02 * .5 * (gamma+.001)**(-.5) - vlapl = None - vtau = None - vxc = (vrho, vgamma, vlapl, vtau) - fxc = None # 2nd order functional derivative - kxc = None # 3rd order functional derivative - return exc, vxc, fxc, kxc - - mf = dft1.RKS(mol) - ni = dft.libxc.define_xc(mf._numint, eval_xc, 'GGA', hyb=0.2) - numpy.random.seed(1) - rho = numpy.random.random((4,10)) - exc, vxc = ni.eval_xc(None, rho, 0, deriv=1)[:2] - self.assertAlmostEqual(lib.fp(exc), 0.0012441814416833327, 9) - self.assertAlmostEqual(lib.fp(vxc[0]), 0.0065565189784811129, 9) - self.assertAlmostEqual(lib.fp(vxc[1]), 0.0049270110162854116, 9) - - with self.assertRaises(AssertionError): - ni.eval_xc_eff(None, rho, deriv=2) - - # Ensure the xc_code in the input has no impact to the result. - n, exc, vxc = ni.nr_rks(mol, mf.grids, 'm06x', dm) - self.assertAlmostEqual(n, 4, 5) - self.assertAlmostEqual(lib.fp(exc), 0.01197588220700074, 6) - self.assertAlmostEqual(lib.fp(vxc), -0.043974912389152986, 6) - - ni = dft.libxc.define_xc(ni, eval_xc, 'MGGA') - with self.assertRaises(AssertionError): - ni.eval_xc_eff(None, rho, deriv=1) - - mf = mf.define_xc_('0.5*B3LYP+0.5*B3LYP') - exc0, vxc0 = mf._numint.eval_xc(None, rho, 0, deriv=1)[:2] - exc1, vxc1 = dft.libxc.eval_xc('0.5*B3LYP+0.5*B3LYP', rho, 0, deriv=1)[:2] - self.assertAlmostEqual(abs(exc0-exc1).max(), 0, 9) - self.assertAlmostEqual(abs(vxc0[0]-vxc1[0]).max(), 0, 9) - self.assertAlmostEqual(abs(vxc0[1]-vxc1[1]).max(), 0, 9) - - self.assertRaises(ValueError, dft.libxc.define_xc, mf._numint, 0.1) - - ni = dft.libxc.define_xc(mf._numint, 'PBE') - exc1, vxc1 = dft.libxc.eval_xc('PBE', rho, 0, deriv=1)[:2] - exc0, vxc0 = ni.eval_xc(None, rho, 0, deriv=1)[:2] - self.assertAlmostEqual(abs(exc0-exc1).max(), 0, 9) - self.assertAlmostEqual(abs(vxc0[0]-vxc1[0]).max(), 0, 9) - self.assertAlmostEqual(abs(vxc0[1]-vxc1[1]).max(), 0, 9) - - def test_define_xc_mgga(self): - def eval_mgga_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None): - # A fictitious XC functional to demonstrate the usage - rho0, dx, dy, dz, tau = rho - gamma = (dx**2 + dy**2 + dz**2) - exc = .01 * rho0**2 + .02 * (gamma+.001)**.5 + .01 * tau**2 - vrho = .01 * 2 * rho0 - vgamma = .02 * .5 * (gamma+.001)**(-.5) - vtau = .02 * tau - vxc = (vrho, vgamma, vtau) - zeros = numpy.zeros_like(rho0) - fxc = (zeros, zeros, zeros) - kxc = None # 3rd order functional derivative - return exc, vxc, fxc, kxc - - mf = dft1.RKS(mol) - ni = dft.libxc.define_xc(mf._numint, eval_mgga_xc, 'MGGA') - numpy.random.seed(1) - rho = numpy.random.random((5,10)) - exc, vxc = ni.eval_xc(None, rho, 0, deriv=1)[:2] - self.assertAlmostEqual(lib.fp(exc), 0.016353861820705, 9) - self.assertAlmostEqual(lib.fp(vxc[0]), 0.0065565189784811129, 9) - self.assertAlmostEqual(lib.fp(vxc[1]), 0.0049270110162854116, 9) - - with self.assertRaises(AssertionError): - ni.eval_xc_eff(None, rho, deriv=2) - - def test_m05x(self): - rho =(numpy.array([1., 1., 0., 0., 0., 0.165 ]).reshape(-1,1), - numpy.array([.8, 1., 0., 0., 0., 0.1050]).reshape(-1,1)) - #test_ref = numpy.array([-1.57876583, -2.12127045,-2.11264351,-0.00315462, - # 0.00000000, -0.00444560, 3.45640232, 4.4349756]) # libxc-4.3.4 - test_ref = numpy.array([-1.57730394, -2.12127045,-2.11297165,-0.00315462, - 0.00000000, -0.00446935, 3.45640232, 4.42563831]) # libxc-5.1.2 - exc, vxc, fxc, kxc = dft.libxc.eval_xc('1.38888888889*m05,', rho, 1, deriv=1) - self.assertAlmostEqual(float(exc[0])*1.8, test_ref[0], 5) - self.assertAlmostEqual(abs(vxc[0]-test_ref[1:3]).max(), 0, 6) - self.assertAlmostEqual(abs(vxc[1]-test_ref[3:6]).max(), 0, 6) - self.assertAlmostEqual(abs(vxc[3]-test_ref[6:8]).max(), 0, 5) - - exc, vxc, fxc, kxc = dft.libxc.eval_xc('1.38888888889*m05,', rho[0], 0, deriv=1) - self.assertAlmostEqual(float(exc[0]), -0.5746231988116002, 5) - self.assertAlmostEqual(float(vxc[0][0]), -0.8806121005703862, 6) - self.assertAlmostEqual(float(vxc[1][0]), -0.0032300155406846756, 7) - self.assertAlmostEqual(float(vxc[3][0]), 0.4474953100487698, 5) - - def test_camb3lyp(self): - rho = numpy.array([1., 1., 0.1, 0.1]).reshape(-1,1) - exc, vxc, fxc, kxc = dft.libxc.eval_xc('camb3lyp', rho, 0, deriv=1) - self.assertAlmostEqual(float(exc[0]), -0.5752559666317147, 7) - self.assertAlmostEqual(float(vxc[0][0]), -0.7709812578936763, 7) - self.assertAlmostEqual(float(vxc[1][0]), -0.0029862221286189846, 7) - - self.assertEqual(dft.libxc.rsh_coeff('camb3lyp'), (0.33, 0.65, -0.46)) - - rho = numpy.array([1., 1., 0.1, 0.1]).reshape(-1,1) - exc, vxc, fxc, kxc = dft.libxc.eval_xc('RSH(0.5,0.65,-0.46) + 0.46*ITYH + .35*B88,', rho, 0, deriv=1) - self.assertAlmostEqual(float(exc[0]), -0.48916154057161476, 9) - self.assertAlmostEqual(float(vxc[0][0]), -0.6761177630311709, 9) - self.assertAlmostEqual(float(vxc[1][0]), -0.002949151742087167, 9) - - def test_ityh(self): - rho = numpy.array([1., 1., 0.1, 0.1]).reshape(-1,1) - exc, vxc, fxc, kxc = dft.libxc.eval_xc('ityh,', rho, 0, deriv=1) - self.assertAlmostEqual(float(exc[0]), -0.6359945579326314, 7) - self.assertAlmostEqual(float(vxc[0][0]), -0.8712041561251518, 7) - self.assertAlmostEqual(float(vxc[1][0]), -0.003911167644579979, 7) - self.assertEqual(dft.libxc.rsh_coeff('ityh,'), (0.2, 0.0, 0.0)) - - def test_lcwpbe(self): - rho = numpy.array([1., 1., 0.1, 0.1]).reshape(-1,1) - exc, vxc, fxc, kxc = dft.libxc.eval_xc('LC_wPBE', rho, 0, deriv=1) - exc1, vxc1, fxc1, kxc1 = dft.libxc.eval_xc('RSH(0.4,1.0,-1.0)+wpbeh,pbe', rho, 0, deriv=1) - self.assertAlmostEqual(float(exc[0]), float(exc1[0]), 7) - self.assertAlmostEqual(float(vxc[0][0]), float(vxc1[0][0]), 7) - self.assertAlmostEqual(float(vxc[1][0]), float(vxc1[1][0]), 7) - - def test_deriv_order(self): - self.assertTrue(dft.libxc.test_deriv_order('lda', 3, raise_error=False)) - self.assertTrue(dft.libxc.test_deriv_order('m05', 2, raise_error=False)) - self.assertTrue(dft.libxc.test_deriv_order('camb3lyp', 3, True)) - self.assertTrue(dft.libxc.test_deriv_order('pbe0', 3, True)) - self.assertRaises(KeyError, dft.libxc.test_deriv_order, 'OL2', 3, True) - - def test_xc_type(self): - self.assertEqual(dft.libxc.xc_type(416), 'GGA') - self.assertEqual(dft.libxc.xc_type('hf'), 'HF') - self.assertEqual(dft.libxc.xc_type(',vwn'), 'LDA') - self.assertEqual(dft.libxc.xc_type('lda+b3lyp'), 'GGA') - self.assertEqual(dft.libxc.xc_type('wb97x_v'), 'GGA') - self.assertEqual(dft.libxc.xc_type('wb97m_v'), 'MGGA') - self.assertEqual(dft.libxc.xc_type('bp86'), 'GGA') - - def test_m06(self): - rho = numpy.array([0.11939647, -0.18865577, -0.11633254, 0.01779666, -0.55475521, 0.07092032, - 0.11850155, -0.19297934, -0.11581427, 0.01373251, -0.57534216, 0.06596468,]) - rho = rho.reshape(2,6,1) - exc, vxc, fxc, kxc = dft.libxc.eval_xc(',m06', rho[0], 0, deriv=3) - exc_ref = numpy.array([-0.02809518]) - vxc_ref = numpy.array([-0.0447466, -0.002186, 0.0154814]) - fxc_ref = numpy.array([-0.66060719, 0.35848677, -1.37754938, - -1.63288299, 0.21999345, 1.20222259]) - kxc_ref = numpy.array([-8.91475141e+00, 2.57071018e+01, -4.37413046e+01, 9.92064877e+01, - -1.47859948e+01, 1.51506522e+01, -4.83492461e+00, -1.75713898e+01, - -2.82029718e+01, 6.66847006e+01]) - self.assertAlmostEqual(exc[0] - exc_ref[0], 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([vxc[i] for i in [0,1,3]])-vxc_ref).max(), 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([fxc[i] for i in [0,1,2,4,6,9]])-fxc_ref).max(), 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([kxc[i] for i in [0,1,2,3,5,7,10,12,15,19]])-kxc_ref).max(), 0, 7) - - exc, vxc, fxc, kxc = dft.libxc.eval_xc(',m06', rho, 1, deriv=3) - exc_ref = numpy.array([-0.03730162]) - vxc_ref = numpy.array([-0.06168318,-0.0610502,-0.00957516,0,-0.01988745, 0.03161512, 0.04391615]) - fxc_ref = numpy.array([ 0.62207817,-1.02032218, 0.59318092, - -0.00267104, 0, 0.2768053, 0.26785479, 0, 0.02509479, - -1.34211352, 0,-0.0399852, 0, 0,-1.28860549, - -1.69697623,-0.03422618,-2.1835459, - -0.22523468, 0.26116507, 0.26545614,-0.23488143, - 1.52817397,-0.13694284, 0, 0,-0.13797703, 1.73484678]) - kxc_ref = numpy.array([ 1.49056834e+01, 4.93235407e+00, 5.13177936e+00, 6.83179846e+00, - -1.97455841e+01, 0,-8.12321586e-01,-1.51226035e+00, 0,-1.66826458e+00,-8.99591711e-01, 0, -9.19638526e+00, - 1.90010442e+01, 0, 1.34796896e-01, 0, 0,-7.73620959e-01,-6.10018848e-01, 0, 2.77414513e-01, 0, 0, 6.16594005e+00, - 2.67162133e+01, 0,-9.09592410e-01, 0, 0,-9.28025786e-01, 0, 0, 0, 3.76462242e+01, - 2.74381367e-01,-2.08192677e+00, 5.25100873e-01, 6.10851890e-01,-2.02662181e+00,-3.07707014e+00, - -4.45125589e-01, 6.81329714e-01, 0, 0,-4.69145558e-01, -2.42540301e+00,-2.45259140e+00,-5.57504449e-01, 0, 0, 6.02633158e-01, 2.58871134e+00, - 2.87554064e+00, 2.85563442e-01,-1.94628304e-01, -1.88482952e-01, 2.90513272e-01, 3.07132889e+00, - -1.08483575e+00, 1.57821239e+00, 0, 0, 1.59013101e+00, 1.61019571e+00, 0, 0, 0, 0, 1.62235586e+00,-3.60275312e+00, - -3.90972925e+01,-1.10619931e-02,-1.12015763e-02, 0, 0, 0,-1.11455329e-02,-1.12861703e-02, -4.81639494e+01, - 7.24911656e+01, 4.35367189e-02, 4.40860771e-02, 1.00010975e+02]) - self.assertAlmostEqual(exc[0] - exc_ref[0], 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([vxc[i] for i in [0,1,3]])-vxc_ref).max(), 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([fxc[i] for i in [0,1,2,4,6,9]])-fxc_ref).max(), 0, 7) - self.assertAlmostEqual(abs(numpy.hstack([kxc[i] for i in [0,1,2,3,5,7,10,12,15,19]])-kxc_ref).max(), 0, 6) - - def test_dft_parser(self): - from pyscf.dft.dft_parser import parse_dft - self.assertEqual(parse_dft('wb97m-d3bj'), ('wb97m-v', False, 'd3bj')) - self.assertEqual(parse_dft('b3lyp-d3zerom'), ('b3lyp', '', 'd3zerom')) - self.assertEqual(parse_dft('wb97x-d3bj'), ('wb97x-v', False, 'd3bj')) - self.assertEqual(parse_dft('wb97x-d3zero2b'), ('wb97x', '', 'd3zero2b')) - - def test_set_param(self): - XC_ID_B97_2 = 410 - - param_b97_1 = numpy.array([0.789518, 0.573805, 0.660975, 0.0, 0.0, - 0.0820011, 2.71681, -2.87103, 0.0, 0.0, - 0.955689, 0.788552, -5.47869, 0.0, 0.0, - 0.21 - ]) - - param_b97_3 = numpy.array([0.7334648, 0.292527, 3.338789, -10.51158, 10.60907, - 0.5623649, -1.32298, 6.359191, -7.464002, 1.827082, - 1.13383, -2.811967, 7.431302, -1.969342, -11.74423, - 2.692880E-01 - ]) - - e1, v1, f1 = dft.libxc.eval_xc('B97-1', rho, 0, deriv=2)[:3] - - e3, v3, f3 = dft.libxc.eval_xc('B97-3', rho, 0, deriv=2)[:3] - - dft.libxc.register_custom_functional_('myfunc1', 'B97-2', - ext_params={XC_ID_B97_2: param_b97_1}) - self.assertAlmostEqual(dft.libxc.hybrid_coeff('myfunc1'), 0.21, 10) - e, v, f = dft.libxc.eval_xc('myfunc1', rho, 0, deriv=2)[:3] - self.assertTrue(numpy.allclose(e, e1, 1e-7)) - self.assertTrue(numpy.allclose(v, v1, 1e-7)) - self.assertTrue(numpy.allclose(f, f1, 1e-7)) - - dft.libxc.register_custom_functional_('myfunc2', 'B97-2', - ext_params={XC_ID_B97_2: param_b97_3}) - self.assertAlmostEqual(dft.libxc.hybrid_coeff('myfunc2'), 2.692880E-01, 10) - e, v, f = dft.libxc.eval_xc('myfunc2', rho, 0, deriv=2)[:3] - self.assertTrue(numpy.allclose(e, e3, 1e-7)) - self.assertTrue(numpy.allclose(v, v3, 1e-7)) - self.assertTrue(numpy.allclose(f, f3, 1e-7)) - - dft.libxc.update_custom_functional_('myfunc1', - ext_params={XC_ID_B97_2: param_b97_3}) - self.assertAlmostEqual(dft.libxc.hybrid_coeff('myfunc1'), 2.692880E-01, 10) - e, v, f = dft.libxc.eval_xc('myfunc1', rho, 0, deriv=2)[:3] - self.assertTrue(numpy.allclose(e, e3, 1e-7)) - self.assertTrue(numpy.allclose(v, v3, 1e-7)) - self.assertTrue(numpy.allclose(f, f3, 1e-7)) - - dft.libxc.unregister_custom_functional_('myfunc1') - dft.libxc.unregister_custom_functional_('myfunc2') - self.assertRaises(KeyError, dft.libxc.xc_type, 'myfunc1') - self.assertRaises(KeyError, dft.libxc.xc_type, 'myfunc2') - -if __name__ == "__main__": - print("Test libxc") - unittest.main() diff --git a/pyscf/dft2/test/test_lpdft.py b/pyscf/dft2/test/test_lpdft.py deleted file mode 100644 index b9b4544e8..000000000 --- a/pyscf/dft2/test/test_lpdft.py +++ /dev/null @@ -1,149 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2023 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Author: Matthew Hennefarth - -import tempfile, h5py -import numpy as np -from pyscf import gto, scf, dft, fci, lib -from pyscf import mcpdft -import unittest - - -def get_lih (r, n_states=2, functional='ftLDA,VWN3', basis='sto3g'): - mol = gto.M (atom='Li 0 0 0\nH {} 0 0'.format (r), basis=basis, - output='/dev/null', verbose=0) - mf = scf.RHF (mol).run () - if n_states == 2: - mc = mcpdft.CASSCF (mf, functional, 2, 2, grids_level=1) - - else: - mc = mcpdft.CASSCF(mf, functional, 5, 2, grids_level=1) - - mc.fix_spin_(ss=0) - weights = [1.0/float(n_states), ] * n_states - - mc = mc.multi_state(weights, "lin") - mc = mc.run() - return mc - -def get_water(functional='tpbe', basis='6-31g'): - mol = gto.M(atom=''' - O 0. 0.000 0.1174 - H 0. 0.757 -0.4696 - H 0. -0.757 -0.4696 - ''',symmetry=True, basis=basis, output='/dev/null', verbose=0) - - mf = scf.RHF(mol).run() - - weights = [0.5, 0.5] - solver1 = fci.direct_spin1_symm.FCI(mol) - solver1.wfnsym = 'A1' - solver1.spin = 0 - solver2 = fci.direct_spin1_symm.FCI(mol) - solver2.wfnsym = 'A2' - solver2.spin = 2 - - mc = mcpdft.CASSCF(mf, functional, 4, 4, grids_level=1) - mc.chkfile = tempfile.NamedTemporaryFile().name - # mc.chk_ci = True - mc = mc.multi_state_mix([solver1, solver2], weights, "lin") - mc.run() - return mc - -def get_water_triplet(functional='tPBE', basis="6-31G"): - mol = gto.M(atom=''' - O 0. 0.000 0.1174 - H 0. 0.757 -0.4696 - H 0. -0.757 -0.4696 - ''', symmetry=True, basis=basis, output='/dev/null', verbose=0) - - mf = scf.RHF(mol).run() - - weights = np.ones(3) / 3 - solver1 = fci.direct_spin1_symm.FCI(mol) - solver1.spin = 2 - solver1 = fci.addons.fix_spin(solver1, shift=.2, ss=2) - solver1.nroots = 1 - solver2 = fci.direct_spin0_symm.FCI(mol) - solver2.spin = 0 - solver2.nroots = 2 - - mc = mcpdft.CASSCF(mf, functional, 4, 4, grids_level=1) - mc.chkfile = tempfile.NamedTemporaryFile().name - # mc.chk_ci = True - mc = mc.multi_state_mix([solver1, solver2], weights, "lin") - mc.run() - return mc - - -def setUpModule(): - global lih, lih_4, lih_tpbe, lih_tpbe0, lih_mc23, water, t_water, original_grids - - from importlib import reload - from pyscf import dft2 - dft.libxc = dft2.libxc - reload (mcpdft) - reload (mcpdft.otfnal) - - original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False - lih = get_lih(1.5) - lih_4 = get_lih(1.5, n_states=4, basis="6-31G") - lih_tpbe = get_lih(1.5, functional="tPBE") - lih_tpbe0 = get_lih(1.5, functional="tPBE0") - lih_mc23 = get_lih(1.5, functional="MC23") - water = get_water() - t_water = get_water_triplet() - -def tearDownModule(): - global lih, lih_4, lih_tpbe0, lih_tpbe, t_water, water, original_grids, lih_mc23 - dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = original_grids - lih.mol.stdout.close() - lih_4.mol.stdout.close() - lih_tpbe0.mol.stdout.close() - lih_tpbe.mol.stdout.close() - water.mol.stdout.close() - t_water.mol.stdout.close() - del lih, lih_4, lih_tpbe0, lih_tpbe, t_water, water, original_grids, lih_mc23 - -class KnownValues(unittest.TestCase): - - def assertListAlmostEqual(self, first_list, second_list, expected): - self.assertTrue(len(first_list) == len(second_list)) - for first, second in zip(first_list, second_list): - self.assertAlmostEqual(first, second, expected) - - def test_lih_mc23_adiabat(self): - e_mcscf_mc23_avg = np.dot(lih_mc23.e_mcscf, lih_mc23.weights) - hcoup = abs(lih_mc23.lpdft_ham[1,0]) - hdiag = lih_mc23.get_lpdft_diag() - - # Reference values from - # - PySCF commit 9a0bb6ddded7049bdacdaf4cfe422f7ce826c2c7 - # - PySCF-forge commit eb0ad96f632994d2d1846009ecce047193682526 - E_MCSCF_AVG_EXPECTED = -7.78902182 - E_MC23_EXPECTED = [-7.94539408, -7.80094952] - HCOUP_EXPECTED = 0.01285147 - HDIAG_EXPECTED = [-7.94424147, -7.80210214] - - self.assertAlmostEqual(e_mcscf_mc23_avg, E_MCSCF_AVG_EXPECTED, 7) - self.assertAlmostEqual(hcoup, HCOUP_EXPECTED, 7) - self.assertAlmostEqual(lib.fp(hdiag), lib.fp(HDIAG_EXPECTED), 7) - self.assertAlmostEqual(lib.fp(lih_mc23.e_states), lib.fp(E_MC23_EXPECTED), 7) - -if __name__ == "__main__": - print("Full Tests for Linearized-PDFT") - unittest.main() diff --git a/pyscf/dft2/test/test_mgga.py b/pyscf/dft2/test/test_mgga.py deleted file mode 100644 index 26cbad5a1..000000000 --- a/pyscf/dft2/test/test_mgga.py +++ /dev/null @@ -1,232 +0,0 @@ -#!/usr/bin/env python -# Copyright 2014-2024 The PySCF Developers. All Rights Reserved. -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Author: Bhavnesh Jangid - -import numpy as np -from pyscf import gto, scf, dft, fci -from pyscf import mcpdft -import unittest - -''' -In this unit-test, test the MCPDFT energies calculated for the LiH -molecule at the state-specific and state-average (2-states) using -1. Meta-GGA functional (tM06L) -2. Hybrid-meta-GGA functional tM06L0 -3. MC23 Functional - -Test the MCPDFT energies calculated for the triplet water molecule at the -4. Meta-GGA functional (M06L) -5. MC23 Functional - -Note: The reference values are generated from -OpenMolcas v24.10, tag 682-gf74be507d - -The OpenMolcas results were obtained with this grid settings -&SEWARD -Grid Input -RQuad=TA -NR=100 -LMAX=41 -NOPrun -NOSCreening -''' - -# To be consistent with OpenMolcas Grid Settings. Grids_att is defined as below -# Source: pyscf/mcpdft/test/test_diatomic_energies.py - -om_ta_alpha = [0.8, 0.9, # H, He - 1.8, 1.4, # Li, Be - 1.3, 1.1, 0.9, 0.9, 0.9, 0.9, # B - Ne - 1.4, 1.3, # Na, Mg - 1.3, 1.2, 1.1, 1.0, 1.0, 1.0, # Al - Ar - 1.5, 1.4, # K, Ca - 1.3, 1.2, 1.2, 1.2, 1.2, 1.2, 1.2, 1.1, 1.1, 1.1, # Sc - Zn - 1.1, 1.0, 0.9, 0.9, 0.9, 0.9] # Ga - Kr - -def om_treutler_ahlrichs(n, chg, *args, **kwargs): - ''' - "Treutler-Ahlrichs" as implemented in OpenMolcas - ''' - r = np.empty(n) - dr = np.empty(n) - alpha = om_ta_alpha[chg-1] - step = 2.0 / (n+1) # = numpy.pi / (n+1) - ln2 = alpha / np.log(2) - for i in range(n): - x = (i+1)*step - 1 # = numpy.cos((i+1)*step) - r [i] = -ln2*(1+x)**.6 * np.log((1-x)/2) - dr[i] = (step #* numpy.sin((i+1)*step) - * ln2*(1+x)**.6 *(-.6/(1+x)*np.log((1-x)/2)+1/(1-x))) - return r[::-1], dr[::-1] - -my_grids = {'atom_grid': (99,590), - 'radi_method': om_treutler_ahlrichs, - 'prune': False, - 'radii_adjust': None} - -def get_lih (r, stateaverage=False, functional='tM06L', basis='sto3g'): - mol = gto.M (atom='Li 0 0 0\nH {} 0 0'.format (r), basis=basis, - output='/dev/null', verbose=0) - mf = scf.RHF (mol).run () - if functional == 'tM06L0': - tM06L0 = 't' + mcpdft.hyb('M06L',0.25, hyb_type='average') - mc = mcpdft.CASSCF(mf, tM06L0, 5, 2, grids_attr=my_grids) - else: - mc = mcpdft.CASSCF(mf, functional, 5, 2, grids_attr=my_grids) - - if stateaverage: - mc = mc.state_average_([0.5, 0.5]) - - mc.fix_spin_(ss=0) - mc = mc.run() - return mc - -def get_water_triplet(functional='tM06L', basis='6-31G'): - mol = gto.M(atom=''' - O 0. 0.000 0.1174 - H 0. 0.757 -0.4696 - H 0. -0.757 -0.4696 - ''',basis=basis, spin=2,output='/dev/null', verbose=0) - - mf = scf.RHF(mol).run() - - mc = mcpdft.CASSCF(mf, functional, 2, 2, grids_attr=my_grids) - solver1 = fci.direct_spin1.FCI(mol) - solver1 = fci.addons.fix_spin(solver1, ss=2) - mc.fcisolver = solver1 - mc = mc.run() - return mc - -def setUpModule(): - global get_lih, lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2,lih_tm06l0 - global get_water_triplet, water_tm06l, water_tmc23 - global lih_tmc23_2, lih_tmc23_sa2_2, water_tmc23_2 - - from importlib import reload - from pyscf import dft2 - dft.libxc = dft2.libxc - reload (mcpdft) - reload (mcpdft.otfnal) - - # register otfnal tMC23_2 which is identical to MC23 - mc232_preset = mcpdft.otfnal.OT_PRESET['MC23'] - mcpdft.otfnal.register_otfnal('MC23_2', mc232_preset) - - lih_tm06l = get_lih(1.5, functional='tM06L') - lih_tmc23 = get_lih(1.5, functional='MC23') - lih_tmc23_2 = get_lih(1.5, functional='tMC23_2') - lih_tm06l_sa2 = get_lih(1.5, stateaverage=True, functional='tM06L') - lih_tmc23_sa2 = get_lih(1.5, stateaverage=True, functional='MC23') - lih_tmc23_sa2_2 = get_lih(1.5, stateaverage=True, functional='tmc23_2') - lih_tm06l0 = get_lih(1.5, functional='tM06L0') - water_tm06l = get_water_triplet() - water_tmc23 = get_water_triplet(functional='MC23') - water_tmc23_2 = get_water_triplet(functional='TMc23_2') - -def tearDownModule(): - global lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 - global lih_tm06l0, water_tm06l, water_tmc23 - global lih_tmc23_2, lih_tmc23_sa2_2, water_tmc23_2 - - lih_tm06l.mol.stdout.close() - lih_tmc23.mol.stdout.close() - lih_tmc23_2.mol.stdout.close() - lih_tm06l_sa2.mol.stdout.close() - lih_tmc23_sa2.mol.stdout.close() - lih_tmc23_sa2_2.mol.stdout.close() - lih_tm06l0.mol.stdout.close() - water_tm06l.mol.stdout.close() - water_tmc23.mol.stdout.close() - water_tmc23_2.mol.stdout.close() - - mcpdft.otfnal.unregister_otfnal('tMC23_2') - - del lih_tm06l, lih_tmc23, lih_tm06l_sa2, lih_tmc23_sa2 - del lih_tm06l0, water_tm06l, water_tmc23 - del lih_tmc23_2, lih_tmc23_sa2_2, water_tmc23_2 - -class KnownValues(unittest.TestCase): - - def assertListAlmostEqual(self, first_list, second_list, expected): - self.assertTrue(len(first_list) == len(second_list)) - for first, second in zip(first_list, second_list): - self.assertAlmostEqual(first, second, expected) - - def test_tmc23(self): - e_mcscf = lih_tmc23.e_mcscf - epdft = lih_tmc23.e_tot - - sa_e_mcscf = lih_tmc23_sa2.e_mcscf - sa_epdft = lih_tmc23_sa2.e_states - - # The CAS and MCPDFT reference values are generated using - # OpenMolcas v24.10, tag 682-gf74be507d - E_CASSCF_EXPECTED = -7.88214917 - E_MCPDFT_EXPECTED = -7.95098727 - SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] - SA_E_MCPDFT_EXPECTED = [-7.95093826, -7.80604012] - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 6) - self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 6) - - def test_tmc23_2(self): - e_mcscf = lih_tmc23_2.e_mcscf - epdft = lih_tmc23_2.e_tot - - sa_e_mcscf = lih_tmc23_sa2_2.e_mcscf - sa_epdft = lih_tmc23_sa2_2.e_states - - # The CAS and MCPDFT reference values are generated using - # OpenMolcas v24.10, tag 682-gf74be507d - E_CASSCF_EXPECTED = -7.88214917 - E_MCPDFT_EXPECTED = -7.95098727 - SA_E_CASSCF_EXPECTED = [-7.88205449, -7.74391704] - SA_E_MCPDFT_EXPECTED = [-7.95093826, -7.80604012] - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - self.assertListAlmostEqual(sa_e_mcscf, SA_E_CASSCF_EXPECTED, 6) - self.assertListAlmostEqual(sa_epdft, SA_E_MCPDFT_EXPECTED, 6) - - def test_water_triplet_tmc23(self): - e_mcscf = water_tmc23.e_mcscf - epdft = water_tmc23.e_tot - - # The CAS and MCPDFT reference values are generated using - # OpenMolcas v24.10, tag 682-gf74be507d - E_CASSCF_EXPECTED = -75.72365496 - E_MCPDFT_EXPECTED = -76.02630019 - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - - def test_water_triplet_tmc23_2(self): - e_mcscf = water_tmc23_2.e_mcscf - epdft = water_tmc23_2.e_tot - - # The CAS and MCPDFT reference values are generated using - # OpenMolcas v24.10, tag 682-gf74be507d - E_CASSCF_EXPECTED = -75.72365496 - E_MCPDFT_EXPECTED = -76.02630019 - - self.assertAlmostEqual(e_mcscf, E_CASSCF_EXPECTED, 6) - self.assertAlmostEqual(epdft, E_MCPDFT_EXPECTED, 6) - -if __name__ == "__main__": - print("Full Tests for MGGAs, Hybrid-MGGAs, and MC23") - unittest.main() diff --git a/pyscf/lib/CMakeLists.txt b/pyscf/lib/CMakeLists.txt index c53aa5088..d20477b6c 100644 --- a/pyscf/lib/CMakeLists.txt +++ b/pyscf/lib/CMakeLists.txt @@ -171,15 +171,6 @@ else() message(STATUS "OCCRI C extension build disabled by user") endif() -# Build the LibXC interface v2 library -set (LIBXC_ITRF2_SOURCE_FILES "dft/libxc_itrf2.c") -add_library (clib_libxc_itrf2 SHARED ${LIBXC_ITRF2_SOURCE_FILES}) -target_link_libraries (clib_libxc_itrf2 xc ${OPENMP_C_PROPERTIES}) -set_target_properties (clib_libxc_itrf2 PROPERTIES - CLEAN_DIRECT_OUTPUT 1 - LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR} - OUTPUT_NAME "xc_itrf2") - # Build the DSRG library set (DSRG_SOURCE_FILES "dsrg/dsrg_helper.c") add_library (clib_dsrg SHARED ${DSRG_SOURCE_FILES}) diff --git a/pyscf/lib/dft/libxc_itrf2.c b/pyscf/lib/dft/libxc_itrf2.c deleted file mode 100644 index 02ebfbf91..000000000 --- a/pyscf/lib/dft/libxc_itrf2.c +++ /dev/null @@ -1,1151 +0,0 @@ -/* Copyright 2014-2018 The PySCF Developers. All Rights Reserved. - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. - - * - * Authors: Qiming Sun - * Susi Lehtola - * Xing Zhang - * Dayou Zhang - * - * libxc from - * https://libxc.gitlab.io - * - * This file is adapted from `lib/dft/libxc_itrf.c` of the PySCF core module - * commit 25eaa9572977b903de24d5c11ad345cecd744728 - */ - -#include -#include -#include -#include -#include -#include "config.h" -#define MAX(X,Y) ((X) > (Y) ? (X) : (Y)) -#define MIN(X,Y) ((X) < (Y) ? (X) : (Y)) -#define MAX_THREADS 256 - -// TODO: register python signal -#define raise_error return - -/* Extracted from libxc:functionals.c since this function is not exposed - * currently. See issue #2756. - */ -static int -xc_func_find_ext_params_name(const xc_func_type *p, const char *name) { - int ii; - assert(p != NULL && p->info->ext_params.n > 0); - for(ii=0; iiinfo->ext_params.n; ii++){ - if(strcmp(p->info->ext_params.names[ii], name) == 0) { - return ii; - } - } - return -1; -} - -/* Extracted from comments of libxc:gga.c - - sigma_st = grad rho_s . grad rho_t - zk = energy density per unit particle - - vrho_s = d n*zk / d rho_s - vsigma_st = d n*zk / d sigma_st - - v2rho2_st = d^2 n*zk / d rho_s d rho_t - v2rhosigma_svx = d^2 n*zk / d rho_s d sigma_tv - v2sigma2_stvx = d^2 n*zk / d sigma_st d sigma_vx - - v3rho3_stv = d^3 n*zk / d rho_s d rho_t d rho_v - v3rho2sigma_stvx = d^3 n*zk / d rho_s d rho_t d sigma_vx - v3rhosigma2_svxyz = d^3 n*zk / d rho_s d sigma_vx d sigma_yz - v3sigma3_stvxyz = d^3 n*zk / d sigma_st d sigma_vx d sigma_yz - - if nspin == 2 - rho(2) = (u, d) - sigma(3) = (uu, ud, dd) - - * vxc(N*5): - vrho(2) = (u, d) - vsigma(3) = (uu, ud, dd) - - * fxc(N*45): - v2rho2(3) = (u_u, u_d, d_d) - v2rhosigma(6) = (u_uu, u_ud, u_dd, d_uu, d_ud, d_dd) - v2sigma2(6) = (uu_uu, uu_ud, uu_dd, ud_ud, ud_dd, dd_dd) - v2lapl2(3) - vtau2(3) - v2rholapl(4) - v2rhotau(4) - v2lapltau(4) - v2sigmalapl(6) - v2sigmatau(6) - - * kxc(N*35): - v3rho3(4) = (u_u_u, u_u_d, u_d_d, d_d_d) - v3rho2sigma(9) = (u_u_uu, u_u_ud, u_u_dd, u_d_uu, u_d_ud, u_d_dd, d_d_uu, d_d_ud, d_d_dd) - v3rhosigma2(12) = (u_uu_uu, u_uu_ud, u_uu_dd, u_ud_ud, u_ud_dd, u_dd_dd, d_uu_uu, d_uu_ud, d_uu_dd, d_ud_ud, d_ud_dd, d_dd_dd) - v3sigma(10) = (uu_uu_uu, uu_uu_ud, uu_uu_dd, uu_ud_ud, uu_ud_dd, uu_dd_dd, ud_ud_ud, ud_ud_dd, ud_dd_dd, dd_dd_dd) - - */ - -#define LDA_NVAR 1 -#define GGA_NVAR 4 -#define MGGA_NVAR 5 - -/* - * rho_u/rho_d = (den,grad_x,grad_y,grad_z,laplacian,tau) - * In spin restricted case (spin == 1), rho_u is assumed to be the - * spin-free quantities, rho_d is not used. - */ -static void _eval_rho(double *restrict rho, double *restrict rho_u, const int spin, - const int nvar, const int np, int ld_rho_u) -{ - int i; - double *sigma, *tau; - double *gxu, *gyu, *gzu, *gxd, *gyd, *gzd; - double *tau_u, *tau_d; - double *rho_d = rho_u + ld_rho_u * nvar; - - switch (nvar) { - case LDA_NVAR: - if (spin == 1) { - for (i = 0; i < np; i++) { - rho[i*2+0] = rho_u[i]; - rho[i*2+1] = rho_d[i]; - } - } else { - for (i = 0; i < np; i++) { - rho[i] = rho_u[i]; - } - } - break; - case GGA_NVAR: - if (spin == 1) { - sigma = rho + np * 2; - gxu = rho_u + ld_rho_u; - gyu = rho_u + ld_rho_u * 2; - gzu = rho_u + ld_rho_u * 3; - gxd = rho_d + ld_rho_u; - gyd = rho_d + ld_rho_u * 2; - gzd = rho_d + ld_rho_u * 3; - for (i = 0; i < np; i++) { - rho[i*2+0] = rho_u[i]; - rho[i*2+1] = rho_d[i]; - sigma[i*3+0] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i]; - sigma[i*3+1] = gxu[i]*gxd[i] + gyu[i]*gyd[i] + gzu[i]*gzd[i]; - sigma[i*3+2] = gxd[i]*gxd[i] + gyd[i]*gyd[i] + gzd[i]*gzd[i]; - } - } else { - sigma = rho + np; - gxu = rho_u + ld_rho_u; - gyu = rho_u + ld_rho_u * 2; - gzu = rho_u + ld_rho_u * 3; - for (i = 0; i < np; i++) { - rho[i] = rho_u[i]; - sigma[i] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i]; - } - } - break; - case MGGA_NVAR: - if (spin == 1) { - sigma = rho + np * 2; - tau = sigma + np * 3; - gxu = rho_u + ld_rho_u; - gyu = rho_u + ld_rho_u * 2; - gzu = rho_u + ld_rho_u * 3; - gxd = rho_d + ld_rho_u; - gyd = rho_d + ld_rho_u * 2; - gzd = rho_d + ld_rho_u * 3; - tau_u = rho_u + ld_rho_u * 4; - tau_d = rho_d + ld_rho_u * 4; - for (i = 0; i < np; i++) { - rho[i*2+0] = rho_u[i]; - rho[i*2+1] = rho_d[i]; - tau[i*2+0] = tau_u[i]; - tau[i*2+1] = tau_d[i]; - } - for (i = 0; i < np; i++) { - sigma[i*3+0] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i]; - sigma[i*3+1] = gxu[i]*gxd[i] + gyu[i]*gyd[i] + gzu[i]*gzd[i]; - sigma[i*3+2] = gxd[i]*gxd[i] + gyd[i]*gyd[i] + gzd[i]*gzd[i]; - } - } else { - sigma = rho + np; - tau = sigma + np; - gxu = rho_u + ld_rho_u; - gyu = rho_u + ld_rho_u * 2; - gzu = rho_u + ld_rho_u * 3; - tau_u = rho_u + ld_rho_u * 4; - for (i = 0; i < np; i++) { - rho[i] = rho_u[i]; - sigma[i] = gxu[i]*gxu[i] + gyu[i]*gyu[i] + gzu[i]*gzu[i]; - tau[i] = tau_u[i]; - } - } - break; - } -} -static void _eval_xc(const xc_func_type *func_x, const int spin, const int deriv, const int np, - double *restrict rho, double *restrict exc, const int offset, const int blksize) -{ - double *sigma, *tau; - double *lapl = rho; - double *vrho = NULL; - double *vsigma = NULL; - double *vlapl = NULL; - double *vtau = NULL; - double *v2rho2 = NULL; - double *v2rhosigma = NULL; - double *v2sigma2 = NULL; - double *v2lapl2 = NULL; - double *v2tau2 = NULL; - double *v2rholapl = NULL; - double *v2rhotau = NULL; - double *v2sigmalapl = NULL; - double *v2sigmatau = NULL; - double *v2lapltau = NULL; - double *v3rho3 = NULL; - double *v3rho2sigma = NULL; - double *v3rhosigma2 = NULL; - double *v3sigma3 = NULL; - double *v3rho2lapl = NULL; - double *v3rho2tau = NULL; - double *v3rhosigmalapl = NULL; - double *v3rhosigmatau = NULL; - double *v3rholapl2 = NULL; - double *v3rholapltau = NULL; - double *v3rhotau2 = NULL; - double *v3sigma2lapl = NULL; - double *v3sigma2tau = NULL; - double *v3sigmalapl2 = NULL; - double *v3sigmalapltau = NULL; - double *v3sigmatau2 = NULL; - double *v3lapl3 = NULL; - double *v3lapl2tau = NULL; - double *v3lapltau2 = NULL; - double *v3tau3 = NULL; - double *v4rho4 = NULL; - double *v4rho3sigma = NULL; - double *v4rho3lapl = NULL; - double *v4rho3tau = NULL; - double *v4rho2sigma2 = NULL; - double *v4rho2sigmalapl = NULL; - double *v4rho2sigmatau = NULL; - double *v4rho2lapl2 = NULL; - double *v4rho2lapltau = NULL; - double *v4rho2tau2 = NULL; - double *v4rhosigma3 = NULL; - double *v4rhosigma2lapl = NULL; - double *v4rhosigma2tau = NULL; - double *v4rhosigmalapl2 = NULL; - double *v4rhosigmalapltau= NULL; - double *v4rhosigmatau2 = NULL; - double *v4rholapl3 = NULL; - double *v4rholapl2tau = NULL; - double *v4rholapltau2 = NULL; - double *v4rhotau3 = NULL; - double *v4sigma4 = NULL; - double *v4sigma3lapl = NULL; - double *v4sigma3tau = NULL; - double *v4sigma2lapl2 = NULL; - double *v4sigma2lapltau = NULL; - double *v4sigma2tau2 = NULL; - double *v4sigmalapl3 = NULL; - double *v4sigmalapl2tau = NULL; - double *v4sigmalapltau2 = NULL; - double *v4sigmatau3 = NULL; - double *v4lapl4 = NULL; - double *v4lapl3tau = NULL; - double *v4lapl2tau2 = NULL; - double *v4lapltau3 = NULL; - double *v4tau4 = NULL; - - switch (func_x->info->family) { - case XC_FAMILY_LDA: -#ifdef XC_FAMILY_HYB_LDA - case XC_FAMILY_HYB_LDA: -#endif - // ex is the energy density - // NOTE libxc library added ex/ec into vrho/vcrho - // vrho = rho d ex/d rho + ex, see work_lda.c:L73 - if (spin == 1) { - if (deriv > 0) { - vrho = exc + np; - } - if (deriv > 1) { - v2rho2 = vrho + np * 2; - } - if (deriv > 2) { - v3rho3 = v2rho2 + np * 3; - } - if (deriv > 3) { - v4rho4 = v3rho3 + np * 4; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset * 2; - } - if (deriv > 1) { - v2rho2 += offset * 3; - } - if (deriv > 2) { - v3rho3 += offset * 4; - } - if (deriv > 3) { - v4rho4 += offset * 5; - } - } else { - if (deriv > 0) { - vrho = exc + np; - } - if (deriv > 1) { - v2rho2 = vrho + np; - } - if (deriv > 2) { - v3rho3 = v2rho2 + np; - } - if (deriv > 3) { - v4rho4 = v3rho3 + np; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset; - } - if (deriv > 1) { - v2rho2 += offset; - } - if (deriv > 2) { - v3rho3 += offset; - } - if (deriv > 3) { - v4rho4 += offset; - } - } - xc_lda(func_x, blksize, rho, exc, vrho, v2rho2, v3rho3, v4rho4); - break; - case XC_FAMILY_GGA: -#ifdef XC_FAMILY_HYB_GGA - case XC_FAMILY_HYB_GGA: -#endif - if (spin == 1) { - sigma = rho + blksize * 2; - if (deriv > 0) { - vrho = exc + np; - vsigma = vrho + np * 2; - } - if (deriv > 1) { - v2rho2 = vsigma + np * 3; - v2rhosigma = v2rho2 + np * 3; - v2sigma2 = v2rhosigma + np * 6; // np*6 - } - if (deriv > 2) { - v3rho3 = v2sigma2 + np * 6; - v3rho2sigma = v3rho3 + np * 4; - v3rhosigma2 = v3rho2sigma + np * 9; - v3sigma3 = v3rhosigma2 + np * 12; // np*10 - } - if (deriv > 3) { - v4rho4 = v3sigma3 + np * 10 ; - v4rho3sigma = v4rho4 + np * 5 ; - v4rho2sigma2 = v4rho3sigma + np * 4*3 ; - v4rhosigma3 = v4rho2sigma2 + np * 3*6 ; - v4sigma4 = v4rhosigma3 + np * 2*10; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset * 2; - vsigma += offset * 3; - } - if (deriv > 1) { - v2rho2 += offset * 3; - v2rhosigma += offset * 6; - v2sigma2 += offset * 6; - } - if (deriv > 2) { - v3rho3 += offset * 4; - v3rho2sigma += offset * 9; - v3rhosigma2 += offset * 12; - v3sigma3 += offset * 10; - } - if (deriv > 3) { - v4rho4 += offset * 5; - v4rho3sigma += offset * 4*3; - v4rho2sigma2 += offset * 3*6; - v4rhosigma3 += offset * 2*10; - v4sigma4 += offset * 15; - } - } else { - sigma = rho + blksize; - if (deriv > 0) { - vrho = exc + np; - vsigma = vrho + np; - } - if (deriv > 1) { - v2rho2 = vsigma + np; - v2rhosigma = v2rho2 + np; - v2sigma2 = v2rhosigma + np; - } - if (deriv > 2) { - v3rho3 = v2sigma2 + np; - v3rho2sigma = v3rho3 + np; - v3rhosigma2 = v3rho2sigma + np; - v3sigma3 = v3rhosigma2 + np; - } - if (deriv > 3) { - v4rho4 = v3sigma3 + np; - v4rho3sigma = v4rho4 + np; - v4rho2sigma2 = v4rho3sigma + np; - v4rhosigma3 = v4rho2sigma2 + np; - v4sigma4 = v4rhosigma3 + np; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset; - vsigma += offset; - } - if (deriv > 1) { - v2rho2 += offset; - v2rhosigma += offset; - v2sigma2 += offset; - } - if (deriv > 2) { - v3rho3 += offset; - v3rho2sigma += offset; - v3rhosigma2 += offset; - v3sigma3 += offset; - } - if (deriv > 3) { - v4rho4 += offset; - v4rho3sigma += offset; - v4rho2sigma2 += offset; - v4rhosigma3 += offset; - v4sigma4 += offset; - } - } - xc_gga(func_x, blksize, rho, sigma, - exc, vrho, vsigma, - v2rho2, v2rhosigma, v2sigma2, - v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3, - v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4); - break; - case XC_FAMILY_MGGA: -#ifdef XC_FAMILY_HYB_MGGA - case XC_FAMILY_HYB_MGGA: -#endif - if (spin == 1) { - sigma = rho + blksize * 2; - tau = sigma + blksize * 3; - if (deriv > 0) { - vrho = exc + np; - vsigma = vrho + np * 2; - vtau = vsigma + np * 3; - } - if (deriv > 1) { - v2rho2 = vtau + np * 2; - v2rhosigma = v2rho2 + np * 3; - v2sigma2 = v2rhosigma + np * 6; - v2rhotau = v2sigma2 + np * 6; - v2sigmatau = v2rhotau + np * 4; - v2tau2 = v2sigmatau + np * 6; - } - if (deriv > 2) { - v3rho3 = v2tau2 + np * 3 ; - v3rho2sigma = v3rho3 + np * 4 ; - v3rhosigma2 = v3rho2sigma + np * 9 ; - v3sigma3 = v3rhosigma2 + np * 12; - v3rho2tau = v3sigma3 + np * 10; - v3rhosigmatau = v3rho2tau + np * 6 ; - v3rhotau2 = v3rhosigmatau + np * 12; - v3sigma2tau = v3rhotau2 + np * 6 ; - v3sigmatau2 = v3sigma2tau + np * 12; - v3tau3 = v3sigmatau2 + np * 9 ; - } - if (deriv > 3) { - v4rho4 = v3tau3 + np * 4 ; - v4rho3sigma = v4rho4 + np * 5 ; - v4rho2sigma2 = v4rho3sigma + np * 4*3 ; - v4rhosigma3 = v4rho2sigma2 + np * 3*6 ; - v4sigma4 = v4rhosigma3 + np * 2*10 ; - v4rho3tau = v4sigma4 + np * 15 ; - v4rho2sigmatau = v4rho3tau + np * 4*2 ; - v4rho2tau2 = v4rho2sigmatau + np * 3*3*2; - v4rhosigma2tau = v4rho2tau2 + np * 3*3 ; - v4rhosigmatau2 = v4rhosigma2tau + np * 2*6*2; - v4rhotau3 = v4rhosigmatau2 + np * 2*3*3; - v4sigma3tau = v4rhotau3 + np * 2*4 ; - v4sigma2tau2 = v4sigma3tau + np * 10*2 ; - v4sigmatau3 = v4sigma2tau2 + np * 6*3 ; - v4tau4 = v4sigmatau3 + np * 3*4 ; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset * 2; - vsigma += offset * 3; - vtau += offset * 2; - } - if (deriv > 1) { - v2rho2 += offset * 3; - v2rhosigma += offset * 6; - v2sigma2 += offset * 6; - v2rhotau += offset * 4; - v2sigmatau += offset * 6; - v2tau2 += offset * 3; - } - if (deriv > 2) { - v3rho3 += offset * 4 ; - v3rho2sigma += offset * 9 ; - v3rhosigma2 += offset * 12; - v3sigma3 += offset * 10; - v3rho2tau += offset * 6 ; - v3rhosigmatau += offset * 12; - v3rhotau2 += offset * 6 ; - v3sigma2tau += offset * 12; - v3sigmatau2 += offset * 9 ; - v3tau3 += offset * 4 ; - } - if (deriv > 3) { - v4rho4 += offset * 5 ; - v4rho3sigma += offset * 4*3 ; - v4rho2sigma2 += offset * 3*6 ; - v4rhosigma3 += offset * 2*10 ; - v4sigma4 += offset * 15 ; - v4rho3tau += offset * 4*2 ; - v4rho2sigmatau += offset * 3*3*2; - v4rho2tau2 += offset * 3*3 ; - v4rhosigma2tau += offset * 2*6*2; - v4rhosigmatau2 += offset * 2*3*3; - v4rhotau3 += offset * 2*4 ; - v4sigma3tau += offset * 10*2 ; - v4sigma2tau2 += offset * 6*3 ; - v4sigmatau3 += offset * 3*4 ; - v4tau4 += offset * 5 ; - } - } else { - sigma = rho + blksize; - tau = sigma + blksize; - if (deriv > 0) { - vrho = exc + np; - vsigma = vrho + np; - vtau = vsigma + np; - } - if (deriv > 1) { - v2rho2 = vtau + np; - v2rhosigma = v2rho2 + np; - v2sigma2 = v2rhosigma + np; - v2rhotau = v2sigma2 + np; - v2sigmatau = v2rhotau + np; - v2tau2 = v2sigmatau + np; - } - if (deriv > 2) { - v3rho3 = v2tau2 + np; - v3rho2sigma = v3rho3 + np; - v3rhosigma2 = v3rho2sigma + np; - v3sigma3 = v3rhosigma2 + np; - v3rho2tau = v3sigma3 + np; - v3rhosigmatau = v3rho2tau + np; - v3rhotau2 = v3rhosigmatau + np; - v3sigma2tau = v3rhotau2 + np; - v3sigmatau2 = v3sigma2tau + np; - v3tau3 = v3sigmatau2 + np; - } - if (deriv > 3) { - v4rho4 = v3tau3 + np; - v4rho3sigma = v4rho4 + np; - v4rho2sigma2 = v4rho3sigma + np; - v4rhosigma3 = v4rho2sigma2 + np; - v4sigma4 = v4rhosigma3 + np; - v4rho3tau = v4sigma4 + np; - v4rho2sigmatau = v4rho3tau + np; - v4rho2tau2 = v4rho2sigmatau + np; - v4rhosigma2tau = v4rho2tau2 + np; - v4rhosigmatau2 = v4rhosigma2tau + np; - v4rhotau3 = v4rhosigmatau2 + np; - v4sigma3tau = v4rhotau3 + np; - v4sigma2tau2 = v4sigma3tau + np; - v4sigmatau3 = v4sigma2tau2 + np; - v4tau4 = v4sigmatau3 + np; - } - - // set offset - exc += offset; - if (deriv > 0) { - vrho += offset; - vsigma += offset; - vtau += offset; - } - if (deriv > 1) { - v2rho2 += offset; - v2rhosigma += offset; - v2sigma2 += offset; - v2rhotau += offset; - v2sigmatau += offset; - v2tau2 += offset; - } - if (deriv > 2) { - v3rho3 += offset; - v3rho2sigma += offset; - v3rhosigma2 += offset; - v3sigma3 += offset; - v3rho2tau += offset; - v3rhosigmatau += offset; - v3rhotau2 += offset; - v3sigma2tau += offset; - v3sigmatau2 += offset; - v3tau3 += offset; - } - if (deriv > 3) { - v4rho4 += offset; - v4rho3sigma += offset; - v4rho2sigma2 += offset; - v4rhosigma3 += offset; - v4sigma4 += offset; - v4rho3tau += offset; - v4rho2sigmatau += offset; - v4rho2tau2 += offset; - v4rhosigma2tau += offset; - v4rhosigmatau2 += offset; - v4rhotau3 += offset; - v4sigma3tau += offset; - v4sigma2tau2 += offset; - v4sigmatau3 += offset; - v4tau4 += offset; - } - } - xc_mgga(func_x, blksize, rho, sigma, lapl, tau, - exc, vrho, vsigma, vlapl, vtau, - v2rho2, v2rhosigma, v2rholapl, v2rhotau, v2sigma2, - v2sigmalapl, v2sigmatau, v2lapl2, v2lapltau, v2tau2, - v3rho3, v3rho2sigma, v3rho2lapl, v3rho2tau, v3rhosigma2, - v3rhosigmalapl, v3rhosigmatau, v3rholapl2, v3rholapltau, - v3rhotau2, v3sigma3, v3sigma2lapl, v3sigma2tau, - v3sigmalapl2, v3sigmalapltau, v3sigmatau2, v3lapl3, - v3lapl2tau, v3lapltau2, v3tau3, - v4rho4, v4rho3sigma, v4rho3lapl, v4rho3tau, v4rho2sigma2, - v4rho2sigmalapl, v4rho2sigmatau, v4rho2lapl2, v4rho2lapltau, - v4rho2tau2, v4rhosigma3, v4rhosigma2lapl, v4rhosigma2tau, - v4rhosigmalapl2, v4rhosigmalapltau, v4rhosigmatau2, - v4rholapl3, v4rholapl2tau, v4rholapltau2, v4rhotau3, - v4sigma4, v4sigma3lapl, v4sigma3tau, v4sigma2lapl2, - v4sigma2lapltau, v4sigma2tau2, v4sigmalapl3, v4sigmalapl2tau, - v4sigmalapltau2, v4sigmatau3, v4lapl4, v4lapl3tau, - v4lapl2tau2, v4lapltau3, v4tau4); - break; - default: - fprintf(stderr, "functional %d '%s' is not implemented\n", - func_x->info->number, func_x->info->name); - raise_error; - } -} - -int LIBXC_is_lda(const xc_func_type *func) -{ - int lda; - switch(func->info->family) - { - case XC_FAMILY_LDA: -#ifdef XC_FAMILY_HYB_LDA - case XC_FAMILY_HYB_LDA: -#endif - lda = 1; - break; - default: - lda = 0; - } - - return lda; -} - -int LIBXC_is_gga(const xc_func_type *func) -{ - int gga; - switch(func->info->family) - { - case XC_FAMILY_GGA: -#ifdef XC_FAMILY_HYB_GGA - case XC_FAMILY_HYB_GGA: -#endif - gga = 1; - break; - default: - gga = 0; - } - - return gga; -} - -int LIBXC_is_meta_gga(const xc_func_type *func) -{ - int mgga; - switch(func->info->family) - { - case XC_FAMILY_MGGA: -#ifdef XC_FAMILY_HYB_MGGA - case XC_FAMILY_HYB_MGGA: -#endif - mgga = 1; - break; - default: - mgga = 0; - } - - return mgga; -} - -int LIBXC_needs_laplacian(const int nfunc, const xc_func_type *func) -{ - int i; - for (i = 0; i < nfunc; i++) { - if (func[i].info->flags & XC_FLAGS_NEEDS_LAPLACIAN) { - return 1; - } - } - return 0; -} - -int LIBXC_is_hybrid(const xc_func_type *func) -{ - int hyb; - -#if XC_MAJOR_VERSION <= 7 - switch(func->info->family) - { -#ifdef XC_FAMILY_HYB_LDA - case XC_FAMILY_HYB_LDA: -#endif - case XC_FAMILY_HYB_GGA: - case XC_FAMILY_HYB_MGGA: - hyb = 1; - break; - default: - hyb = 0; - } -#else - hyb = (xc_hyb_type(func) == XC_HYB_HYBRID); -#endif - - return hyb; -} - -double LIBXC_hybrid_coeff(const xc_func_type *func) -{ - double factor; - -#if XC_MAJOR_VERSION <= 7 - switch(func->info->family) - { -#ifdef XC_FAMILY_HYB_LDA - case XC_FAMILY_HYB_LDA: -#endif - case XC_FAMILY_HYB_GGA: - case XC_FAMILY_HYB_MGGA: - factor = xc_hyb_exx_coef(func); - break; - default: - factor = 0; - } - -#else - if(xc_hyb_type(func) == XC_HYB_HYBRID) - factor = xc_hyb_exx_coef(func); - else - factor = 0.0; -#endif - - return factor; -} - -void LIBXC_nlc_coeff(const xc_func_type *func, double *nlc_pars) { - - XC(nlc_coef)(func, &nlc_pars[0], &nlc_pars[1]); -} - -void LIBXC_rsh_coeff(const xc_func_type *func, double *rsh_pars) { - - rsh_pars[0] = 0.0; - rsh_pars[1] = 0.0; - rsh_pars[2] = 0.0; - -#if XC_MAJOR_VERSION <= 7 - XC(hyb_cam_coef)(func, &rsh_pars[0], &rsh_pars[1], &rsh_pars[2]); -#else - switch(xc_hyb_type(func)) { - case(XC_HYB_HYBRID): - case(XC_HYB_CAM): - XC(hyb_cam_coef)(func, &rsh_pars[0], &rsh_pars[1], &rsh_pars[2]); - } -#endif -} - -int LIBXC_is_cam_rsh(const xc_func_type *func) { -#if XC_MAJOR_VERSION <= 7 - return func->info->flags & XC_FLAGS_HYB_CAM; -#else - return (xc_hyb_type(func) == XC_HYB_CAM); -#endif -} - -/* - * XC_FAMILY_LDA 1 - * XC_FAMILY_GGA 2 - * XC_FAMILY_MGGA 4 - * XC_FAMILY_LCA 8 - * XC_FAMILY_OEP 16 - * XC_FAMILY_HYB_GGA 32 - * XC_FAMILY_HYB_MGGA 64 - * XC_FAMILY_HYB_LDA 128 - */ -int LIBXC_xc_type(const int nfunc, const xc_func_type *func) -{ - int i; - int type = -1; - for (i = 0; i < nfunc; i++) { - int t; - switch(func[i].info->family) - { - case XC_FAMILY_LDA: -#ifdef XC_FAMILY_HYB_LDA - case XC_FAMILY_HYB_LDA: -#endif - t = 0; - break; - case XC_FAMILY_GGA: -#ifdef XC_FAMILY_HYB_GGA - case XC_FAMILY_HYB_GGA: -#endif - t = 1; - break; - case XC_FAMILY_MGGA: -#ifdef XC_FAMILY_HYB_MGGA - case XC_FAMILY_HYB_MGGA: -#endif - t = 2; - break; - default: - return -1; - } - type = MAX(t, type); - } - return type; -} - -//static int xc_output_length(int nvar, int deriv) -//{ -// int i; -// int len = 1; -// for (i = 1; i <= nvar; i++) { -// len *= deriv + i; -// len /= i; -// } -// return len; -//} -// offsets = [xc_output_length(nvar, i) for i in range(deriv+1) -// for nvar in [1,2,3,5,7]] -static const int xc_nvar1_offsets[] = {0, 1, 2, 3, 4, 5}; -static const int xc_nvar2_offsets[] = {0, 1, 3, 6, 10, 15}; -static const int xc_nvar3_offsets[] = {0, 1, 4, 10, 20, 35}; -static const int xc_nvar5_offsets[] = {0, 1, 6, 21, 56, 126}; -static const int xc_nvar7_offsets[] = {0, 1, 8, 36, 120, 330}; - -static void axpy(double *restrict dst, double *restrict src, const double fac, - const int np, const int nsrc) -{ - int i, j; - for (j = 0; j < nsrc; j++) { - #pragma omp parallel for schedule(static) - for (i = 0; i < np; i++) { - dst[j*np+i] += fac * src[i*nsrc+j]; - } - } -} -static const int vseg1[] = {2, 3, 2}; -static const int fseg1[] = {3, 6, 6, 4, 6, 3}; -static const int kseg1[] = {4, 9, 12, 10, 6, 12, 6, 12, 9, 4}; -static const int lseg1[] = {5, 12, 18, 20, 15, 8, 18, 9, 24, 18, 8, 20, 18, 12, 5}; -static const int *seg1[] = {NULL, vseg1, fseg1, kseg1, lseg1}; - -static void merge_xc(double *restrict dst, double *restrict ebuf, const double fac, - const int spin, const int deriv, const int nvar, const int np, - const int outlen, const int type) -{ - int order, nsrc, i; - for (i = 0; i < np; i++) { - dst[i] += fac * ebuf[i]; - } - - const int *offsets0, *offsets1; - double *pout, *pin; - switch (type) { - case XC_FAMILY_GGA: -#ifdef XC_FAMILY_HYB_GGA - case XC_FAMILY_HYB_GGA: -#endif - offsets0 = xc_nvar2_offsets; - break; - case XC_FAMILY_MGGA: -#ifdef XC_FAMILY_HYB_MGGA - case XC_FAMILY_HYB_MGGA: -#endif - offsets0 = xc_nvar3_offsets; - break; - default: //case XC_FAMILY_LDA: - offsets0 = xc_nvar1_offsets; - } - - if (spin == 0) { - switch (nvar) { - case LDA_NVAR: - offsets1 = xc_nvar1_offsets; - break; - case GGA_NVAR: - offsets1 = xc_nvar2_offsets; - break; - default: // MGGA_NVAR - offsets1 = xc_nvar3_offsets; - break; - } - - for (order = 1; order <= deriv; order++) { - pout = dst + offsets1[order] * np; - pin = ebuf + offsets0[order] * np; - nsrc = offsets0[order+1] - offsets0[order]; - #pragma omp parallel for schedule(static) - for (i = 0; i < np * nsrc; i++) { - pout[i] += fac * pin[i]; - } - } - return; - } - - switch (nvar) { - case LDA_NVAR: - offsets1 = xc_nvar2_offsets; - break; - case GGA_NVAR: - offsets1 = xc_nvar5_offsets; - break; - default: // MGGA_NVAR - offsets1 = xc_nvar7_offsets; - break; - } - - int terms; - const int *pseg1; - pin = ebuf + np; - for (order = 1; order <= deriv; order++) { - pseg1 = seg1[order]; - pout = dst + offsets1[order] * np; - terms = offsets0[order+1] - offsets0[order]; - for (i = 0; i < terms; i++) { - nsrc = pseg1[i]; - axpy(pout, pin, fac, np, nsrc); - pin += np * nsrc; - pout += np * nsrc; - } - } -} - -// omega is the range separation parameter mu in xcfun -void LIBXC_eval_xc(const int nfn, xc_func_type *fn_obj, const double *fac, const int spin, - const int deriv, const int nvar, const int np, const int outlen, - double *restrict rho_u, double *restrict output) -{ - assert(deriv <= 4); - double *ebuf = malloc(sizeof(double) * np * outlen); - - double *rhobufs[MAX_THREADS]; - int offsets[MAX_THREADS+1]; -#pragma omp parallel -{ - const int iblk = omp_get_thread_num(); - const int nblk = omp_get_num_threads(); - assert(nblk <= MAX_THREADS); - - int blksize = np / nblk; - int ioff = iblk * blksize; - const int np_mod = np % nblk; - if (iblk < np_mod) { - blksize += 1; - } - if (np_mod > 0) { - ioff += MIN(iblk, np_mod); - } - offsets[iblk] = ioff; - if (iblk == nblk-1) { - offsets[nblk] = np; - assert(ioff + blksize == np); - } - - double *rho_priv = malloc(sizeof(double) * blksize * 7); - rhobufs[iblk] = rho_priv; - _eval_rho(rho_priv, rho_u+ioff, spin, nvar, blksize, np); -} - - int i; - xc_func_type *func = fn_obj; - for (i = 0; i < nfn; i++) { -#pragma omp parallel -{ - const int iblk = omp_get_thread_num(); - const int offset = offsets[iblk]; - const int blksize = offsets[iblk+1] - offset; - _eval_xc(func, spin, deriv, np, rhobufs[iblk], ebuf, offset, blksize); -} - - merge_xc(output, ebuf, fac[i], - spin, deriv, nvar, np, outlen, func->info->family); - ++func; - } - free(ebuf); -#pragma omp parallel -{ - const int iblk = omp_get_thread_num(); - free(rhobufs[iblk]); -} -} - -int LIBXC_max_deriv_order(const int nfunc, const xc_func_type *func) -{ - int ord = 4; - int i, o; - const int DERIV_FLAGS_TABLE[] = { - XC_FLAGS_HAVE_EXC, /* order 0 */ - XC_FLAGS_HAVE_VXC, /* order 1 */ - XC_FLAGS_HAVE_FXC, /* order 2 */ - XC_FLAGS_HAVE_KXC, /* order 3 */ - XC_FLAGS_HAVE_LXC /* order 4 */ - }; - - for (i = 0; i < nfunc; i++) { - /* find the minimum order of all functionals */ - const int flag = func[i].info->flags; - for (o = ord; o > 0; o--) { - if (flag & DERIV_FLAGS_TABLE[o]) { - ord = o; - break; - } - } - if (o == -1) return -1; - } - - return ord; -} - -void LIBXC_xc_reference(const int nfn, const xc_func_type *func, const char **refs) -{ - int i, f; - int ref_count = 0; - const char **current_ref = refs; - for (f = 0; f < nfn; f++) { - for (i = 0; ref_count < XC_MAX_REFERENCES * nfn - 1; i++) { - if (func[f].info->refs[i] == NULL || func[f].info->refs[i]->ref == NULL) { - break; - } - *current_ref++ = func[f].info->refs[i]->ref; - ref_count++; - } - } - *current_ref = NULL; -} - -int LIBXC_is_nlc(const int nfunc, const xc_func_type *func) -{ - int i; - for (i = 0; i < nfunc; i++) { - if (func[i].info->flags & XC_FLAGS_VV10) { - return 1; - } - } - return 0; -} - -void LIBXC_xc_func_end(const int nfn, xc_func_type *func) -{ - int i; - for (i = 0; i < nfn; i++) { - xc_func_end(func + i); - } - free(func); -} - -void LIBXC_xc_func_set_params(const int nfn, xc_func_type *fn_obj, const double *omega, - const double dens_threshold) -{ - int i, j; - xc_func_type *func = fn_obj; - for (i = 0; i < nfn; i++) { - if (dens_threshold > 0) { - xc_func_set_dens_threshold(func, dens_threshold); - } - - // set the range-separated parameter - if (omega[i] != 0) { - // skip if func is not a RSH functional - if ( xc_func_find_ext_params_name(func, "_omega") >= 0 ) { - xc_func_set_ext_params_name(func, "_omega", omega[i]); - } - // Recursively set the sub-functionals if they are RSH - // functionals - for (j = 0; j < func->n_func_aux; j++) { - if ( xc_func_find_ext_params_name(func->func_aux[j], "_omega") >= 0 ) { - xc_func_set_ext_params_name(func->func_aux[j], "_omega", omega[i]); - } - } - } - - // alpha and beta are hardcoded in many functionals in the libxc - // code, e.g. the coefficients of B88 (=1-alpha) and - // ITYH (=-beta) in cam-b3lyp. Overwriting func->cam_alpha and - // func->cam_beta does not update the coefficients accordingly. - //func->cam_alpha = alpha; - //func->cam_beta = beta; - // However, the parameters can be set with the libxc function - //void xc_func_set_ext_params_name(xc_func_type *p, const char *name, double par); - // since libxc 5.1.0 -#if defined XC_SET_RELATIVITY - xc_lda_x_set_params(func, relativity); -#endif - ++func; - } -} - -xc_func_type *LIBXC_xc_func_init(const int nfn, const int *xc_ids, const int spin) -{ - int i; - const int s = spin > 0 ? XC_POLARIZED : XC_UNPOLARIZED; - xc_func_type *func; - func = malloc(nfn * sizeof(xc_func_type)); - for (i = 0; i < nfn; i++) { - const int xc_id = xc_ids[i]; - if (xc_func_init(func + i, xc_id, s) != 0){ - fprintf(stderr, "XC functional %d not found\n", xc_id); - LIBXC_xc_func_end(i, func); - raise_error NULL; - } - } - return func; -} - -size_t LIBXC_xc_func_type_size() -{ - return sizeof(xc_func_type); -} diff --git a/pyscf/mcdcft/dcfnal.py b/pyscf/mcdcft/dcfnal.py index 98c754f5a..f2e941f92 100644 --- a/pyscf/mcdcft/dcfnal.py +++ b/pyscf/mcdcft/dcfnal.py @@ -31,7 +31,7 @@ from pyscf.lib import logger from pyscf import dft, lib -from pyscf.dft2 import libxc +from pyscf.dft import libxc import numpy as np import copy