Merge pull request #131 from qsimulate/more-orbitals

More orbitals

Author: ncrubin

src/fqe/bitstring.py

@@ -24,7 +24,7 @@
 
 from fqe.lib.bitstring import _lexicographic_bitstring_generator, _count_bits, \
                               _get_occupation
-import fqe.settings
+from fqe.settings import use_accelerated_code, c_string_max_norb
 
 
 def gbit_index(str0: int) -> Generator[int, None, None]:
@@ -68,7 +68,7 @@ def integer_index(string: int) -> 'Nparray':
         Nparray: a list of integers indicating the index of each occupied \
             orbital
     """
-    if fqe.settings.use_accelerated_code:
+    if use_accelerated_code:
         return _get_occupation(string)
     else:
         return numpy.array(list(gbit_index(int(string)))).astype(numpy.int32)
@@ -107,7 +107,7 @@ def lexicographic_bitstring_generator(nele: int, norb: int) -> 'Nparray':
     if nele > norb:
         raise ValueError("nele cannot be larger than norb")
 
-    if fqe.settings.use_accelerated_code:
+    if 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 +127,7 @@ def count_bits(string: int) -> int:
     Returns:
         int: the number of bits equal to 1
     """
-    if fqe.settings.use_accelerated_code:
+    if use_accelerated_code:
         return _count_bits(string)
     else:
         return bin(int(string)).count('1')
@@ -144,7 +144,7 @@ def get_bit(string: int, pos: int) -> int:
     Returns:
         int: 0 if the bit is 0, 2**pos if the bit is 1
     """
-    return int(string) & (1 << pos)
+    return int(string) & (1 << int(pos))
 
 
 def set_bit(string: int, pos: int) -> int:
@@ -158,7 +158,7 @@ def set_bit(string: int, pos: int) -> int:
     Returns:
         int: string with the pos bit set to 1
     """
-    return int(string) | (1 << pos)
+    return int(string) | (1 << int(pos))
 
 
 def unset_bit(string: int, pos: int) -> int:
@@ -172,7 +172,7 @@ def unset_bit(string: int, pos: int) -> int:
     Returns:
         int: string with the pos bit set to 0
     """
-    return int(string) & ~(1 << pos)
+    return int(string) & ~(1 << int(pos))
 
 
 def count_bits_above(string: int, pos: int) -> int:
@@ -186,7 +186,7 @@ def count_bits_above(string: int, pos: int) -> int:
     Returns:
         int: the number of 1 bits above pos
     """
-    return count_bits(int(string) & ~((1 << (pos + 1)) - 1))
+    return count_bits(int(string) & ~((1 << (int(pos) + 1)) - 1))
 
 
 def count_bits_below(string: int, pos: int) -> int:
@@ -200,7 +200,7 @@ def count_bits_below(string: int, pos: int) -> int:
     Returns:
         int: the number of 1 bits below pos
     """
-    return count_bits(int(string) & ((1 << pos) - 1))
+    return count_bits(int(string) & ((1 << int(pos)) - 1))
 
 
 def count_bits_between(string: int, pos1: int, pos2: int) -> int:

src/fqe/fci_graph.py

@@ -28,7 +28,7 @@
                               _calculate_string_address, \
                               _c_map_to_deexc_alpha_icol, \
                               _make_mapping_each, _calculate_Z_matrix
-import fqe.settings
+from fqe.settings import use_accelerated_code, c_string_max_norb
 
 Spinmap = Dict[Tuple[int, ...], Nparray]
 
@@ -54,7 +54,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 fqe.settings.use_accelerated_code:
+        if use_accelerated_code and norbs <= c_string_max_norb:
             _map_deexc(dexc, values, index, idx)
         else:
             for state, target, parity in values:
@@ -82,7 +82,7 @@ def _get_Z_matrix(norb: int, nele: int) -> Nparray:
     if Z.size == 0:
         return Z
 
-    if fqe.settings.use_accelerated_code:
+    if use_accelerated_code and norb <= c_string_max_norb:
         _calculate_Z_matrix(Z, norb, nele)
     else:
         for k in range(1, nele):
@@ -217,7 +217,7 @@ def _build_mapping(self, strings: Nparray, nele: int,
         """
         norb = self._norb
 
-        if fqe.settings.use_accelerated_code:
+        if use_accelerated_code and norb <= c_string_max_norb:
             return _build_mapping_strings(strings, _get_Z_matrix(norb, nele),
                                           nele, norb)
         else:
@@ -315,7 +315,7 @@ def _build_strings(self, nele: int,
         norb = self._norb
         blist = lexicographic_bitstring_generator(nele, norb)
 
-        if fqe.settings.use_accelerated_code:
+        if 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 +470,7 @@ def _map_to_deexc_alpha_icol(self):
             length2,
         ), dtype=numpy.int32)
         astrings = self.string_alpha_all()
-        if fqe.settings.use_accelerated_code:
+        if 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 +538,7 @@ def make_mapping_each(self, result: 'Nparray', alpha: bool, dag: List[int],
             strings = self.string_beta_all()
             length = self.lenb()
 
-        if fqe.settings.use_accelerated_code:
+        if use_accelerated_code:
             count = _make_mapping_each(result, strings, length,
                                        numpy.array(dag, dtype=numpy.int32),
                                        numpy.array(undag, dtype=numpy.int32))

src/fqe/fci_graph_set.py

@@ -23,7 +23,7 @@
 import scipy
 from scipy import special
 
-import fqe.settings
+from fqe.settings import use_accelerated_code, 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 +171,7 @@ def _postprocess(spinmap, dnv, index0, index1):
 
         if dna != 0:
             (iasec, jasec) = (isec, jsec) if dna < 0 else (jsec, isec)
-            if fqe.settings.use_accelerated_code:
+            if 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 +192,7 @@ def _postprocess(spinmap, dnv, index0, index1):
 
         if dnb != 0:
             (ibsec, jbsec) = (isec, jsec) if dnb < 0 else (jsec, isec)
-            if fqe.settings.use_accelerated_code:
+            if use_accelerated_code and norb <= c_string_max_norb:
                 ndownb, nupb = _make_mapping_each_set(ibsec.string_beta_all(),
                                                       abs(dnb), norb,
                                                       ibsec.nbeta())

src/fqe/fqe_data.py

@@ -802,7 +802,7 @@ def _apply_array_spatial12_lowfilling_fast(self, h1e: 'Nparray',
         nlt = norb * (norb + 1) // 2
 
         h2ecomp = numpy.zeros((nlt, nlt), dtype=self._dtype)
-        h2etemp = numpy.ascontiguousarray(h2e)
+        h2etemp = numpy.ascontiguousarray(h2e, dtype=self._dtype)
         _make_Hcomp(norb, nlt, h2etemp, h2ecomp)
 
         if nalpha - 2 >= 0:

src/fqe/lib/bitstring.c

@@ -33,12 +33,12 @@ void lexicographic_bitstring_generator(uint64_t *out,
     out[0] = 0;
     return;
   }
-  int64_t combo = (1ll << nele) - 1;
-  while (combo < (1ll << norb)) {
+  uint64_t combo = (1ull << nele) - 1;
+  while (combo < (1ull << norb)) {
     *out++ = combo;
-    const int64_t x = combo & -combo;
-    const int64_t y = combo + x;
-    const int64_t z = (combo & ~y);
+    const uint64_t x = combo & -combo;
+    const uint64_t y = combo + x;
+    const uint64_t z = (combo & ~y);
     combo = z / x;
     combo >>= 1;
     combo |= y;

src/fqe/lib/bitstring.h

@@ -27,11 +27,19 @@
 inline int gbit_index(uint64_t *str, int *bit_index) {
 #ifdef __GNUC__
   const int pos = __builtin_ffsll(*str);
-  *str >>= pos;
+  if (pos == 64) {
+    *str = 0ull;
+  } else {
+    *str >>= pos;
+  }
   *bit_index += pos;
   return pos;
 #else
-  if (*bit_index == -1) { *str = *str << 1; }
+  if (*bit_index == -1) {
+    //*str = *str << 1;
+    *bit_index += 1;
+    if (*str & 1) { return 1; }
+  }
   while  (*str) {
     *str = *str >> 1;
     *bit_index += 1;

src/fqe/settings.py

@@ -38,3 +38,14 @@ class CodePath(Flag):
 """
 A switch to check if accelerated code is used. Default is true if accelerated code is available
 """
+
+global_max_norb = 64
+"""
+The global maximum number of orbitals that can be handled by FQE.
+"""
+
+c_string_max_norb = 63
+"""
+The max number of orbitals correctly handled by the C codepath
+for generating determinant strings.
+"""