diff --git a/src/fqe/bitstring.py b/src/fqe/bitstring.py index 97ee938..543e01a 100755 --- a/src/fqe/bitstring.py +++ b/src/fqe/bitstring.py @@ -22,9 +22,10 @@ from numpy import ndarray as Nparray from scipy import special +import fqe from fqe.lib.bitstring import _lexicographic_bitstring_generator, _count_bits, \ _get_occupation -from fqe.settings import use_accelerated_code, c_string_max_norb +from fqe.settings import c_string_max_norb def gbit_index(str0: int) -> Generator[int, None, None]: @@ -68,7 +69,7 @@ def integer_index(string: int) -> 'Nparray': Nparray: a list of integers indicating the index of each occupied \ orbital """ - if use_accelerated_code: + if fqe.settings.use_accelerated_code: return _get_occupation(string) else: return numpy.array(list(gbit_index(int(string)))).astype(numpy.int32) @@ -107,7 +108,7 @@ def lexicographic_bitstring_generator(nele: int, norb: int) -> 'Nparray': if nele > norb: raise ValueError("nele cannot be larger than norb") - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: out = numpy.zeros((int(special.comb(norb, nele)),), dtype=numpy.uint64) _lexicographic_bitstring_generator(out, norb, nele) return out @@ -127,7 +128,7 @@ def count_bits(string: int) -> int: Returns: int: the number of bits equal to 1 """ - if use_accelerated_code: + if fqe.settings.use_accelerated_code: return _count_bits(string) else: return bin(int(string)).count('1') diff --git a/src/fqe/fci_graph.py b/src/fqe/fci_graph.py index dc282fc..789c074 100755 --- a/src/fqe/fci_graph.py +++ b/src/fqe/fci_graph.py @@ -22,13 +22,14 @@ import numpy from numpy import ndarray as Nparray +import fqe from fqe.bitstring import integer_index, lexicographic_bitstring_generator, \ get_bit, set_bit, unset_bit, count_bits_above, count_bits_between from fqe.lib.fci_graph import _build_mapping_strings, _map_deexc, \ _calculate_string_address, \ _c_map_to_deexc_alpha_icol, \ _make_mapping_each, _calculate_Z_matrix -from fqe.settings import use_accelerated_code, c_string_max_norb +from fqe.settings import c_string_max_norb Spinmap = Dict[Tuple[int, ...], Nparray] @@ -54,7 +55,7 @@ def map_to_deexc(mappings: Spinmap, states: int, norbs: int, index = numpy.zeros((states,), dtype=numpy.uint32) for (i, j), values in mappings.items(): idx = i * norbs + j - if use_accelerated_code and norbs <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norbs <= c_string_max_norb: _map_deexc(dexc, values, index, idx) else: for state, target, parity in values: @@ -82,7 +83,7 @@ def _get_Z_matrix(norb: int, nele: int) -> Nparray: if Z.size == 0: return Z - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: _calculate_Z_matrix(Z, norb, nele) else: for k in range(1, nele): @@ -217,7 +218,7 @@ def _build_mapping(self, strings: Nparray, nele: int, """ norb = self._norb - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: return _build_mapping_strings(strings, _get_Z_matrix(norb, nele), nele, norb) else: @@ -315,7 +316,7 @@ def _build_strings(self, nele: int, norb = self._norb blist = lexicographic_bitstring_generator(nele, norb) - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: Z = _get_Z_matrix(norb, nele) string_list = _calculate_string_address(Z, nele, norb, blist) else: @@ -470,7 +471,7 @@ def _map_to_deexc_alpha_icol(self): length2, ), dtype=numpy.int32) astrings = self.string_alpha_all() - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: _c_map_to_deexc_alpha_icol(exc, diag, index, astrings, norb, self._alpha_map) else: @@ -538,7 +539,7 @@ def make_mapping_each(self, result: 'Nparray', alpha: bool, dag: List[int], strings = self.string_beta_all() length = self.lenb() - if use_accelerated_code: + if fqe.settings.use_accelerated_code: count = _make_mapping_each(result, strings, length, numpy.array(dag, dtype=numpy.int32), numpy.array(undag, dtype=numpy.int32)) diff --git a/src/fqe/fci_graph_set.py b/src/fqe/fci_graph_set.py index 6ad7655..b5884a2 100644 --- a/src/fqe/fci_graph_set.py +++ b/src/fqe/fci_graph_set.py @@ -23,7 +23,8 @@ import scipy from scipy import special -from fqe.settings import use_accelerated_code, c_string_max_norb +import fqe +from fqe.settings import c_string_max_norb from fqe.fci_graph import FciGraph from fqe.util import alpha_beta_electrons from fqe.bitstring import integer_index, count_bits_between @@ -171,7 +172,7 @@ def _postprocess(spinmap, dnv, index0, index1): if dna != 0: (iasec, jasec) = (isec, jsec) if dna < 0 else (jsec, isec) - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: ndowna, nupa = _make_mapping_each_set(iasec.string_alpha_all(), abs(dna), norb, iasec.nalpha()) @@ -192,7 +193,7 @@ def _postprocess(spinmap, dnv, index0, index1): if dnb != 0: (ibsec, jbsec) = (isec, jsec) if dnb < 0 else (jsec, isec) - if use_accelerated_code and norb <= c_string_max_norb: + if fqe.settings.use_accelerated_code and norb <= c_string_max_norb: ndownb, nupb = _make_mapping_each_set(ibsec.string_beta_all(), abs(dnb), norb, ibsec.nbeta()) diff --git a/src/fqe/wavefunction.py b/src/fqe/wavefunction.py index f5bccae..4081506 100644 --- a/src/fqe/wavefunction.py +++ b/src/fqe/wavefunction.py @@ -28,7 +28,7 @@ import pickle import numpy from scipy import linalg -from scipy.special import factorial, jv +from scipy.special import jv from fqe.fqe_decorators import wrap_apply, wrap_apply_generated_unitary from fqe.fqe_decorators import wrap_time_evolve, wrap_rdm @@ -555,12 +555,12 @@ def apply_generated_unitary(self, if algo == 'taylor': ham_arrays = hamil.iht(time) - time_evol = copy.deepcopy(base) work = copy.deepcopy(base) + coeff = 1.0 for order in range(1, max_expansion): work = work.apply(ham_arrays) - coeff = 1.0 / factorial(order) + coeff /= order time_evol.ax_plus_y(coeff, work) if work.norm() * numpy.abs(coeff) < accuracy: break @@ -1231,7 +1231,7 @@ def _evolve_individual_nbody(self, raise ValueError( 'Operators in _evolve_individual_nbody is not Hermitian' ) - else: + elif hamil.nterms() == 1: [ (coeff0, alpha0, beta0), ] = hamil.terms() @@ -1247,6 +1247,9 @@ def _evolve_individual_nbody(self, 'Operators in _evolve_individual_nbody is not Hermitian' ) coeff0 *= 0.5 + else: + assert hamil.nterms() == 0 + return self if inplace else copy.deepcopy(self) daga = [] dagb = [] diff --git a/tests/adapt_vqe_test.py b/tests/adapt_vqe_test.py index 3b81cb9..4b4a00c 100644 --- a/tests/adapt_vqe_test.py +++ b/tests/adapt_vqe_test.py @@ -101,4 +101,4 @@ def test_adapt(): verbose=True) adapt.adapt_vqe(fqe_wf) assert np.isclose(adapt.energies[-1], -8.957417182801091) - assert np.isclose(adapt.energies[0], -8.95741717733075) \ No newline at end of file + assert np.isclose(adapt.energies[0], -8.95741717733075) diff --git a/tests/bitstring_test.py b/tests/bitstring_test.py index d340b81..b375328 100755 --- a/tests/bitstring_test.py +++ b/tests/bitstring_test.py @@ -21,9 +21,10 @@ import fqe.settings -def test_bit_integer_index_val(): +def test_bit_integer_index_val(c_or_python): """The index of bits should start at 0. """ + fqe.settings.use_accelerated_code = c_or_python _gbiti = bitstring.gbit_index(1) assert next(_gbiti) == 0 _gbiti = bitstring.gbit_index(11) diff --git a/tests/fci_graph_test.py b/tests/fci_graph_test.py index a5bac8b..d0cd0af 100755 --- a/tests/fci_graph_test.py +++ b/tests/fci_graph_test.py @@ -15,10 +15,14 @@ """ import numpy -import fqe import pytest from scipy import special + +import fqe from fqe import fci_graph +from fqe.lib.fci_graph import _c_map_to_deexc_alpha_icol +from fqe.settings import CodePath + from tests.unittest_data.fci_graph_data import loader from tests.comparisons import compare_Spinmap @@ -311,6 +315,14 @@ def test_map_to_deexc_alpha_icol(c_or_python, norb, nalpha, nbeta): assert numpy.array_equal(rexc, exc) assert numpy.array_equal(rdiag, diag) + if c_or_python == CodePath.PYTHON: + return + + tmp = numpy.zeros((norb, norb)) + test = {(i, i): numpy.zeros((1, 1)) for i in range(norb)} + with pytest.raises(ValueError): + _c_map_to_deexc_alpha_icol(tmp, tmp, tmp, tmp, norb, test) + @pytest.mark.parametrize("nalpha,nbeta,norb", cases) def test_get_block_mappings(norb, nalpha, nbeta): diff --git a/tests/linalg_test.py b/tests/linalg_test.py index fd234b9..3536251 100644 --- a/tests/linalg_test.py +++ b/tests/linalg_test.py @@ -16,10 +16,22 @@ import pytest import numpy import fqe +from fqe.lib import c_double_complex from fqe.lib.linalg import _zimatadd, _transpose from fqe.settings import CodePath +def test_ctype_double_complex(c_or_python): + if c_or_python == CodePath.PYTHON: + # No test needed for python + return + + rval = 1.0 + ival = -1.0 + test = c_double_complex(rval, ival) + assert test.value == rval + 1.j * ival + + def test_zimatadd(c_or_python): """Testing zimatadd""" if c_or_python == CodePath.PYTHON: diff --git a/tests/unittest_data/wavefunction/020206.pickle b/tests/unittest_data/wavefunction/020206.pickle index f84b6c1..8c2b084 100644 Binary files a/tests/unittest_data/wavefunction/020206.pickle and b/tests/unittest_data/wavefunction/020206.pickle differ diff --git a/tests/unittest_data/wavefunction/03-307.pickle b/tests/unittest_data/wavefunction/03-307.pickle index 8be731a..80b84ba 100644 Binary files a/tests/unittest_data/wavefunction/03-307.pickle and b/tests/unittest_data/wavefunction/03-307.pickle differ diff --git a/tests/unittest_data/wavefunction/03XX07.pickle b/tests/unittest_data/wavefunction/03XX07.pickle index 1fb1782..0c62260 100644 Binary files a/tests/unittest_data/wavefunction/03XX07.pickle and b/tests/unittest_data/wavefunction/03XX07.pickle differ diff --git a/tests/unittest_data/wavefunction/04XX06.pickle b/tests/unittest_data/wavefunction/04XX06.pickle index 94429e4..aa93414 100644 Binary files a/tests/unittest_data/wavefunction/04XX06.pickle and b/tests/unittest_data/wavefunction/04XX06.pickle differ diff --git a/tests/unittest_data/wavefunction/070108.pickle b/tests/unittest_data/wavefunction/070108.pickle index 676813b..e6978d5 100644 Binary files a/tests/unittest_data/wavefunction/070108.pickle and b/tests/unittest_data/wavefunction/070108.pickle differ diff --git a/tests/unittest_data/wavefunction/080006.pickle b/tests/unittest_data/wavefunction/080006.pickle index c644279..4e9ec36 100644 Binary files a/tests/unittest_data/wavefunction/080006.pickle and b/tests/unittest_data/wavefunction/080006.pickle differ diff --git a/tests/unittest_data/wavefunction/XX0004.pickle b/tests/unittest_data/wavefunction/XX0004.pickle index 0d81aa8..5604190 100644 Binary files a/tests/unittest_data/wavefunction/XX0004.pickle and b/tests/unittest_data/wavefunction/XX0004.pickle differ diff --git a/tests/unittest_data/wavefunction/XX0102.pickle b/tests/unittest_data/wavefunction/XX0102.pickle index a9331ec..bf280fc 100644 Binary files a/tests/unittest_data/wavefunction/XX0102.pickle and b/tests/unittest_data/wavefunction/XX0102.pickle differ diff --git a/tests/unittest_data/wavefunction/generate_data.py b/tests/unittest_data/wavefunction/generate_data.py index e917a56..b01e976 100644 --- a/tests/unittest_data/wavefunction/generate_data.py +++ b/tests/unittest_data/wavefunction/generate_data.py @@ -65,6 +65,7 @@ def generate_data(param): cr = rng.uniform(-0.5, 0.5, size=value.coeff.shape) ci = rng.uniform(-0.5, 0.5, size=value.coeff.shape) value.coeff = cr + 1.j * ci + wfn.normalize() data['wfn'] = wfn if not spin_conserving: @@ -115,8 +116,8 @@ def generate_data(param): # array ham full = not (spin_conserving and number_conserving) - h1e = build_H1(norbs, full=full) - h2e = build_H2(norbs, full=full) + h1e = build_H1(norbs, full=full) / 10 + h2e = build_H2(norbs, full=full) / 20 e_0 = rng.uniform() + rng.uniform() * 1j if not full: hamil = get_restricted_hamiltonian(( @@ -139,6 +140,7 @@ def generate_data(param): except ValueError: # Don't generate an output evolved = None + #except RuntimeError data['apply_array'] = { 'hamil': hamil, diff --git a/tests/wavefunction_test.py b/tests/wavefunction_test.py index 8688b77..ddf80f4 100755 --- a/tests/wavefunction_test.py +++ b/tests/wavefunction_test.py @@ -15,18 +15,16 @@ """ # pylint: disable=protected-access -import os -import sys import copy -import numpy -import pytest - +import os from io import StringIO +import sys -from scipy.special import binom - +import numpy from openfermion import FermionOperator from openfermion.utils import hermitian_conjugated +import pytest +from scipy.special import binom import fqe from fqe.wavefunction import Wavefunction @@ -116,7 +114,7 @@ def test_empty_copy(): test1 = test.empty_copy() assert test1._norb == test._norb assert len(test1._civec) == len(test._civec) - assert (2, 0) in test1._civec.keys() + assert (2, 0) in test1._civec assert not numpy.any(test1._civec[(2, 0)].coeff) @@ -160,6 +158,10 @@ def test_apply_type_error(): with pytest.raises(TypeError): wfn.time_evolve(0.1, hamil) + hamil2 = restricted_hamiltonian.RestrictedHamiltonian((data,)) + with pytest.raises(TypeError, match="expansion must be an int"): + wfn.apply_generated_unitary(0.1, 'taylor', hamil2, expansion=0.5) + def test_apply_value_error(): data = numpy.zeros((2, 2), dtype=numpy.complex128) @@ -299,6 +301,60 @@ def test_apply_empty_nbody(): assert (wfn - out1).norm() < 1e-13 +def test_expansion_failure(): + """Test exceptions when the max expansion limit is reached. + """ + norb = 4 + nalpha = 2 + nbeta = 1 + nele = nalpha + nbeta + time = 0.1 + + h1e = build_hamiltonian.build_H1(norb, full=True, asymmetric=True) + + eig, _ = numpy.linalg.eigh(h1e) + hamil = fqe.get_general_hamiltonian((h1e,)) + + wfn = Wavefunction([[nele, nalpha - nbeta, norb]]) + wfn.set_wfn(strategy='random') + + err = "maximum chebyshev expansion limit reached" + with pytest.raises(RuntimeError, match=err): + wfn.apply_generated_unitary(time, + 'chebyshev', + hamil, + expansion=1, + accuracy=1e-9, + spec_lim=(eig[0], eig[-1])) + + err = "maximum taylor expansion limit reached" + with pytest.raises(RuntimeError, match=err): + wfn.apply_generated_unitary(time, + 'taylor', + hamil, + expansion=1, + accuracy=1e-9) + + +def test_empty_nbody_evolve(): + """Check that '_evolve_individual_nbody' works with an empty + FermionOperator/SparseHamiltonian + """ + norb = 4 + nele = 4 + time = 0.1 + ops = FermionOperator() + sham = fqe.get_sparse_hamiltonian(ops, conserve_spin=False) + + wfn = fqe.get_number_conserving_wavefunction(nele, norb) + wfn.set_wfn(strategy='random') + wfn.normalize() + + out = wfn._evolve_individual_nbody(time, sham) + + assert (out - wfn).norm() < 1.e-8 + + def test_nbody_evolve(): """Check that 'individual_nbody' is consistent with a general 1-electron op """ @@ -626,7 +682,6 @@ def test_apply(c_or_python, param, kind): assert Wavefunction_isclose(out, reference_data['wfn_out']) -@pytest.mark.skip(reason="Test seems to be failing due to bad reference_data") @pytest.mark.parametrize("param,kind", [(c, k) for c in all_cases for k in [ 'apply_array', 'apply_sparse', 'apply_diagonal', 'apply_quadratic', 'apply_dc'