Now ApqIJcontainer seems to be set up.

src/determinants/configurations.irp.f

@@ -1,623 +0,0 @@
-use bitmasks
-subroutine configuration_of_det(d,o,Nint)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Transforms a determinant to a configuration
-  !
-  ! occ(:,1) : Single occupations
-  !
-  ! occ(:,2) : Double occupations
-  !
-  END_DOC
-  integer          ,intent(in)   :: Nint
-  integer(bit_kind),intent(in)   :: d(Nint,2)
-  integer(bit_kind),intent(out)  :: o(Nint,2)
-
-  integer                        :: k
-
-  do k=1,Nint
-    o(k,1) = ieor(d(k,1),d(k,2))
-    o(k,2) = iand(d(k,1),d(k,2))
-  enddo
-end
-
-
-subroutine configuration_to_dets_size(o,sze,n_alpha,Nint)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-!  Number of possible determinants for a given configuration
-  END_DOC
-  integer          ,intent(in)   :: Nint, n_alpha
-  integer(bit_kind),intent(in)   :: o(Nint,2)
-  integer, intent(out)           :: sze
-  integer                        :: amax,bmax,k
-  double precision, external     :: binom_func
-
-  bmax = 0
-  amax = n_alpha
-  do k=1,Nint
-    bmax += popcnt( o(k,1) )
-    amax -= popcnt( o(k,2) )
-  enddo
-  if (binom_int(bmax, amax) > huge(1)) then
-    print *, bmax, amax
-    print *,  irp_here, ': Too many determinants to generate'
-    stop 1
-  endif
-  sze = int(binom_int(bmax, amax),4)
-end
-
-
-subroutine configuration_to_dets(o,d,sze,n_alpha,Nint)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Generate all possible determinants for a given configuration
-  !
-  ! Input :
-  !    o   : configuration : (doubly occupied, singly occupied)
-  !    sze : Number of produced determinants, computed by `configuration_to_dets_size`
-  !    n_alpha : Number of $\alpha$ electrons
-  !    Nint    : N_int
-  !
-  ! Output:
-  !    d : determinants
-  !
-  END_DOC
-  integer          ,intent(in)   :: Nint
-  integer          ,intent(in)   :: n_alpha        ! Number of alpha electrons
-  integer         ,intent(inout) :: sze            ! Dimension of the output dets
-  integer(bit_kind),intent(in)   :: o(Nint,2)      ! Configurations
-  integer(bit_kind),intent(out)  :: d(Nint,2,sze)  ! Output determinants
-
-  integer                        :: i, k, n, ispin, ispin2
-
-  ! Extract list of singly occupied MOs as (int,pos) pairs
-  ! ------------------------------------------------------
-
-  integer           :: iint(2*n_alpha), ipos(2*n_alpha)
-  integer(bit_kind) :: v, t, tt, diff, v_prev
-  integer           :: n_alpha_in_single
-
-  n=0
-  n_alpha_in_single = n_alpha
-  do i=1,Nint
-    v = o(i,1)
-    do while(v /= 0_bit_kind)
-      n = n+1
-      iint(n) = i
-      ipos(n) = trailz(v)
-      v = iand(v,v-1)
-    enddo
-    n_alpha_in_single = n_alpha_in_single - popcnt( o(i,2) )
-  enddo
-
-  v = shiftl(1,n_alpha_in_single) - 1
-
-  ! Initialize first determinant
-  d(:,1,1) = o(:,2)
-  d(:,2,1) = o(:,2)
-
-  do k=1,n_alpha_in_single
-    d(iint(k),1,1) = ibset( d(iint(k),1,1), ipos(k) )
-  enddo
-
-  do k=n_alpha_in_single+1,n
-    d(iint(k),2,1) = ibset( d(iint(k),2,1), ipos(k) )
-  enddo
-
-  sze = int(binom_int(n,n_alpha_in_single),4)
-
-  if ( (shiftl(n_alpha_in_single,1) == n).and.n>0 ) then
-
-    ! Time reversal symmetry
-    d(:,1,2) = d(:,2,1)
-    d(:,2,2) = d(:,1,1)
-
-    do i=3,sze,2
-      ! Generate next permutation with Anderson's algorithm
-      v_prev = v
-      t = ior(v,v-1)
-      tt = t+1
-      v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) )
-
-      ! Find what has changed between v_prev and v
-      diff = ieor(v,v_prev)
-
-      ! Initialize with previous determinant
-      d(:,1,i) = d(:,1,i-2)
-      d(:,2,i) = d(:,2,i-2)
-
-      ! Swap bits only where they have changed from v_prev to v
-      do while (diff /= 0_bit_kind)
-        k = trailz(diff)+1
-        if (btest(v,k-1)) then
-          d(iint(k),1,i) = ibset( d(iint(k),1,i), ipos(k) )
-          d(iint(k),2,i) = ibclr( d(iint(k),2,i), ipos(k) )
-        else
-          d(iint(k),1,i) = ibclr( d(iint(k),1,i), ipos(k) )
-          d(iint(k),2,i) = ibset( d(iint(k),2,i), ipos(k) )
-        endif
-        diff = iand(diff,diff-1_bit_kind)
-      enddo
-
-      ! Time reversal symmetry
-      d(:,1,i+1) = d(:,2,i)
-      d(:,2,i+1) = d(:,1,i)
-
-    enddo
-
-  else
-
-    do i=2,sze
-      ! Generate next permutation with Anderson's algorithm
-      v_prev = v
-      t = ior(v,v-1)
-      tt = t+1
-      v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) )
-
-      ! Find what has changed between v_prev and v
-      diff = ieor(v,v_prev)
-
-      ! Initialize with previous determinant
-      d(:,1,i) = d(:,1,i-1)
-      d(:,2,i) = d(:,2,i-1)
-
-      ! Swap bits only where they have changed from v_prev to v
-      do while (diff /= 0_bit_kind)
-        k = trailz(diff)+1
-        if (btest(v,k-1)) then
-          d(iint(k),1,i) = ibset( d(iint(k),1,i), ipos(k) )
-          d(iint(k),2,i) = ibclr( d(iint(k),2,i), ipos(k) )
-        else
-          d(iint(k),1,i) = ibclr( d(iint(k),1,i), ipos(k) )
-          d(iint(k),2,i) = ibset( d(iint(k),2,i), ipos(k) )
-        endif
-        diff = iand(diff,diff-1_bit_kind)
-      enddo
-
-    enddo
-
-  endif
-
-end
-
-
- BEGIN_PROVIDER [ integer(bit_kind), psi_configuration, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ integer, N_configuration ]
- implicit none
- BEGIN_DOC
-  ! Array of the configurations present in the wave function.
-  !
-  ! psi_configuration(:,1,j) = j-th configuration of the wave function : represents all the single occupations
-  !
-  ! psi_configuration(:,2,j) = j-th configuration of the wave function : represents all the double occupations
-  !
-  ! The occ patterns are sorted by :c:func:`configuration_search_key`
- END_DOC
- integer :: i,j,k
-
- ! create
- do i = 1, N_det
-  do k = 1, N_int
-   psi_configuration(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i))
-   psi_configuration(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i))
-  enddo
- enddo
-
- ! Sort
- integer, allocatable           :: iorder(:)
- integer*8, allocatable         :: bit_tmp(:)
- integer*8, external            :: configuration_search_key
- integer(bit_kind), allocatable :: tmp_array(:,:,:)
- logical,allocatable            :: duplicate(:)
- logical :: dup
-
-
- allocate ( iorder(N_det), duplicate(N_det), bit_tmp(N_det), tmp_array(N_int,2,N_det) )
-
- do i=1,N_det
-   iorder(i) = i
-   bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
- enddo
-
- call i8sort(bit_tmp,iorder,N_det)
-
-
- !$OMP PARALLEL DEFAULT(shared) PRIVATE(i,j,k,dup)
-
- !$OMP DO
- do i=1,N_det
-  do k=1,N_int
-    tmp_array(k,1,i) = psi_configuration(k,1,iorder(i))
-    tmp_array(k,2,i) = psi_configuration(k,2,iorder(i))
-  enddo
-  duplicate(i) = .False.
- enddo
- !$OMP END DO
-
- ! Find duplicates
- !$OMP DO
- do i=1,N_det-1
-  if (duplicate(i)) then
-    cycle
-  endif
-  j = i+1
-  do while (bit_tmp(j)==bit_tmp(i))
-    if (duplicate(j)) then
-      j+=1
-      if (j>N_det) then
-        exit
-      endif
-      cycle
-    endif
-    dup = .True.
-    do k=1,N_int
-      dup = dup .and. (tmp_array(k,1,i) == tmp_array(k,1,j)) &
-                .and. (tmp_array(k,2,i) == tmp_array(k,2,j))
-    enddo
-    if (dup) then
-      duplicate(j) = .True.
-    endif
-    j = j+1
-    if (j>N_det) then
-      exit
-    endif
-  enddo
- enddo
- !$OMP END DO
- !$OMP END PARALLEL
-
- ! Copy filtered result
- N_configuration=0
- do i=1,N_det
-  if (duplicate(i)) then
-    cycle
-  endif
-  N_configuration += 1
-  do k=1,N_int
-    psi_configuration(k,1,N_configuration) = tmp_array(k,1,i)
-    psi_configuration(k,2,N_configuration) = tmp_array(k,2,i)
-  enddo
- enddo
-
-!- Check
-!  print *,  'Checking for duplicates in configuration'
-!  do i=1,N_configuration
-!   do j=i+1,N_configuration
-!     duplicate(1) = .True.
-!     do k=1,N_int
-!       if (psi_configuration(k,1,i) /= psi_configuration(k,1,j)) then
-!         duplicate(1) = .False.
-!         exit
-!       endif
-!       if (psi_configuration(k,2,i) /= psi_configuration(k,2,j)) then
-!         duplicate(1) = .False.
-!         exit
-!       endif
-!     enddo
-!     if (duplicate(1)) then
-!       call debug_det(psi_configuration(1,1,i),N_int)
-!       call debug_det(psi_configuration(1,1,j),N_int)
-!       stop 'DUPLICATE'
-!     endif
-!   enddo
-!  enddo
-!  print *,  'No duplicates'
-!-
- deallocate(iorder,duplicate,bit_tmp,tmp_array)
-
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ]
-  implicit none
-  BEGIN_DOC
- ! Returns the index in psi_configuration of the first cfg with
- ! the requested seniority
- END_DOC
- integer :: i, k, s, sold
- cfg_seniority_index(:) = -1
- sold = -1
- do i=1,N_configuration
-   s = 0
-   do k=1,N_int
-     if (psi_configuration(k,1,i) == 0_bit_kind) cycle
-     s = s + popcnt(psi_configuration(k,1,i))
-   enddo
-   if (s /= sold) then
-     sold = s
-     cfg_seniority_index(s) = i
-   endif
- enddo
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
- implicit none
- BEGIN_DOC
- ! Returns the index of the configuration for each determinant
- END_DOC
- integer :: i,j,k,r,l
- integer*8 :: key, key2
- integer(bit_kind) :: occ(N_int,2)
- logical :: found
- integer*8, allocatable :: bit_tmp(:)
- integer*8, external            :: configuration_search_key
-
- allocate(bit_tmp(0:N_configuration))
- bit_tmp(0) = 0
- do i=1,N_configuration
-   bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
- enddo
-
- !$OMP PARALLEL DO DEFAULT(SHARED) &
- !$OMP PRIVATE(i,k,j,r,l,key,found,occ)
- do i=1,N_det
-    do k = 1, N_int
-      occ(k,1) = ieor(psi_det(k,1,i),psi_det(k,2,i))
-      occ(k,2) = iand(psi_det(k,1,i),psi_det(k,2,i))
-    enddo
-
-    key = configuration_search_key(occ,N_int)
-
-    ! Binary search
-    l = 0
-    r = N_configuration+1
-    j = shiftr(r-l,1)
-    do while (j>=1)
-      j = j+l
-      if (bit_tmp(j) == key) then
-        do while (bit_tmp(j) == bit_tmp(j-1))
-          j = j-1
-        enddo
-        do while (bit_tmp(j) == key)
-          found = .True.
-          do k=1,N_int
-            found = found .and. (psi_configuration(k,1,j) == occ(k,1)) &
-                          .and. (psi_configuration(k,2,j) == occ(k,2))
-          enddo
-          if (found) then
-            det_to_configuration(i) = j
-            exit
-          endif
-          j = j+1
-        enddo
-        if (found) exit
-      else if (bit_tmp(j) > key) then
-        r = j
-      else
-        l = j
-      endif
-      j = shiftr(r-l,1)
-    enddo
-
- enddo
- !$OMP END PARALLEL DO
- deallocate(bit_tmp)
-END_PROVIDER
-
-
-BEGIN_PROVIDER [ double precision, psi_configuration_Hii, (N_configuration) ]
- implicit none
- BEGIN_DOC
- ! $\langle I|H|I \rangle$ where $|I\rangle$ is a configuration.
- ! This is the minimum $H_{ii}$, where the $|i\rangle$ are the
- ! determinants of $|I\rangle$.
- END_DOC
- integer :: j, i
-
- psi_configuration_Hii(:) = huge(1.d0)
- do i=1,N_det
-  j = det_to_configuration(i)
-  psi_configuration_Hii(j) = min(psi_configuration_Hii(j), psi_det_Hii(i))
- enddo
-
-END_PROVIDER
-
-
-BEGIN_PROVIDER [ double precision, weight_configuration, (N_configuration,N_states) ]
- implicit none
- BEGIN_DOC
- ! Weight of the configurations in the wave function
- END_DOC
- integer :: i,j,k
- weight_configuration = 0.d0
- do i=1,N_det
-  j = det_to_configuration(i)
-  do k=1,N_states
-    weight_configuration(j,k) += psi_coef(i,k) * psi_coef(i,k)
-  enddo
- enddo
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, weight_configuration_average, (N_configuration) ]
- implicit none
- BEGIN_DOC
- ! State-average weight of the configurations in the wave function
- END_DOC
- integer :: i,j,k
- weight_configuration_average(:) = 0.d0
- do i=1,N_det
-  j = det_to_configuration(i)
-  do k=1,N_states
-    weight_configuration_average(j) += psi_coef(i,k) * psi_coef(i,k) * state_average_weight(k)
-  enddo
- enddo
-END_PROVIDER
-
- BEGIN_PROVIDER [ integer(bit_kind), psi_configuration_sorted, (N_int,2,N_configuration) ]
-&BEGIN_PROVIDER [ double precision, weight_configuration_average_sorted, (N_configuration) ]
-&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order, (N_configuration) ]
-&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order_reverse, (N_configuration) ]
- implicit none
- BEGIN_DOC
- ! Configurations sorted by weight
- END_DOC
- integer                        :: i,j,k
- integer, allocatable           :: iorder(:)
- allocate ( iorder(N_configuration) )
- do i=1,N_configuration
-   weight_configuration_average_sorted(i) = -weight_configuration_average(i)
-   iorder(i) = i
- enddo
- call dsort(weight_configuration_average_sorted,iorder,N_configuration)
- do i=1,N_configuration
-   do j=1,N_int
-     psi_configuration_sorted(j,1,i) = psi_configuration(j,1,iorder(i))
-     psi_configuration_sorted(j,2,i) = psi_configuration(j,2,iorder(i))
-   enddo
-   psi_configuration_sorted_order(iorder(i)) = i
-   psi_configuration_sorted_order_reverse(i) = iorder(i)
-   weight_configuration_average_sorted(i) = -weight_configuration_average_sorted(i)
- enddo
-
- deallocate(iorder)
-
-END_PROVIDER
-
-
-subroutine make_s2_eigenfunction
-  implicit none
-  integer                        :: i,j,k
-  integer                        :: smax, s
-  integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:)
-  integer                        :: N_det_new, ithread, omp_get_thread_num
-  integer, parameter             :: bufsze = 1000
-  logical, external              :: is_in_wavefunction
-  logical                        :: update
-
-  update=.False.
-  call write_int(6,N_configuration,'Number of configurations')
-
-  !$OMP PARALLEL DEFAULT(NONE) &
-  !$OMP  SHARED(N_configuration, psi_configuration, elec_alpha_num,N_int,update) &
-  !$OMP  PRIVATE(s,ithread, d, det_buffer, smax, N_det_new,i,j,k)
-  N_det_new = 0
-  call configuration_to_dets_size(psi_configuration(1,1,1),s,elec_alpha_num,N_int)
-  allocate (d(N_int,2,s+64), det_buffer(N_int,2,bufsze) )
-  smax = s
-  ithread=0
-  !$ ithread = omp_get_thread_num()
-  !$OMP DO SCHEDULE (dynamic,1000)
-  do i=1,N_configuration
-    call configuration_to_dets_size(psi_configuration(1,1,i),s,elec_alpha_num,N_int)
-    s += 1
-    if (s > smax) then
-      deallocate(d)
-      allocate ( d(N_int,2,s+64) )
-      smax = s
-    endif
-    call configuration_to_dets(psi_configuration(1,1,i),d,s,elec_alpha_num,N_int)
-    do j=1,s
-      if ( is_in_wavefunction(d(1,1,j), N_int) ) then
-        cycle
-      endif
-      update = .true.
-      N_det_new += 1
-      det_buffer(:,:,N_det_new) = d(:,:,j)
-      if (N_det_new == bufsze) then
-        call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread)
-        N_det_new = 0
-      endif
-    enddo
-  enddo
-  !$OMP END DO NOWAIT
-
-  if (N_det_new > 0) then
-    call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,ithread)
-  endif
-  !$OMP BARRIER
-  deallocate(d,det_buffer)
-  !$OMP END PARALLEL
-
-  if (update) then
-    call copy_H_apply_buffer_to_wf
-    TOUCH N_det psi_coef psi_det psi_configuration N_configuration
-  endif
-  call write_time(6)
-
-end
-
-
-
-BEGIN_PROVIDER [ integer, dominant_cfg, (N_states) ]
- implicit none
- BEGIN_DOC
- ! Configuration of the determinants with the largest weight, for each state
- END_DOC
- integer :: k
- do k=1,N_states
-   dominant_cfg(k) = det_to_configuration(dominant_det(k))
- enddo
-END_PROVIDER
-
-
-BEGIN_PROVIDER [ integer, N_dominant_dets_of_cfgs ]
- implicit none
- BEGIN_DOC
- ! Number of determinants in all dominant determinants
- END_DOC
- integer                        :: k, sze
-
- N_dominant_dets_of_cfgs = 0
- do k=1,N_states
-   call configuration_to_dets_size(  &
-          psi_configuration(1,1,dominant_cfg(k)), &
-          sze, elec_alpha_num, N_int)
-   N_dominant_dets_of_cfgs += sze
- enddo
-END_PROVIDER
-
-BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_dets_of_cfgs) ]
- implicit none
- BEGIN_DOC
- ! Configuration of the determinants with the largest weight, for each state
- END_DOC
- integer :: i,k,sze
- i=1
- do k=1,N_states
-   sze = N_dominant_dets_of_cfgs
-   call configuration_to_dets( &
-          psi_configuration(1,1,dominant_cfg(k)), &
-          dominant_dets_of_cfgs(1,1,i), &
-          sze,elec_alpha_num,N_int)
-   i += sze
- enddo
-END_PROVIDER
-
- BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
-&BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]
-
- implicit none
- BEGIN_DOC
- ! psi_configuration_to_psi_det_data(k) -> i : i is the index of the
- ! determinant in psi_det
- !
- ! psi_configuration_to_psi_det(1:2,k) gives the first and last index of the
- ! determinants of configuration k in array psi_configuration_to_psi_det_data.
- END_DOC
-
- integer :: i, k, iorder
- integer, allocatable :: confs(:)
- allocate (confs(N_det))
-
- do i=1,N_det
-   psi_configuration_to_psi_det_data(i) = i
-   confs(i) = det_to_configuration(i)
- enddo
-
- call isort(confs, psi_configuration_to_psi_det_data, N_det)
- k=1
- psi_configuration_to_psi_det(1,1) = 1
- do i=2,N_det
-   if (confs(i) /= confs(i-1)) then
-     psi_configuration_to_psi_det(2,k) = i-1
-     k = k+1
-     psi_configuration_to_psi_det(1,k) = i
-   endif
- enddo
- psi_configuration_to_psi_det(2,k) = N_det
-
-END_PROVIDER
-

src/determinants/configurations.org

@@ -1,194 +0,0 @@
-# -*- mode:org -*-
-#+TITLE: CFG-CI
-#+AUTHOR: Vijay Gopal Chilkuri
-#+FILE: configurations.org
-#+EMAIL: vijay.gopal.c@gmail.com
-#+OPTIONS: toc:t
-#+LATEX_CLASS: article
-#+LATEX_HEADER: \usepackage{tabularx}
-#+LATEX_HEADER: \usepackage{braket}
-#+LATEX_HEADER: \usepackage{minted}
-
-* Configuration based CI
-
-Here we write the main functions that perform the functions necessary for
-the Configuration based CI.
-
-There are three main functions required for doing the CI
-
-- Convert wavefunction from determinant basis to configuration state function (CSF) basis
-
-- Apply the Hamiltonian to the wavefunction in CSF basis
-
-- Convert the converged wavefunction back to determinant basis
-
-** TODO[0/3] Convert basis from DET to CSF
-
-The conversion of basis is done by going via bonded functions (BFs).
-Importantly, all the CSFs of a chosen configuration (CFG) are kept.
-
-The advantage is that the sigma-vector can be performed efficiently
-via BLAS level 3 operations.
-
-
-- [ ] Write a function to calculate the maximum dimensions required
-
-  Prototype array contains the \( <I|\hat{E}_{pq}|J> \) for all possible
-  CFGs \( I, J\) and all \(4\) types of excitations for all possible model
-  orbitals \(p,q\). Note that the orbital ids \(p,q\) here do not refer to
-  the actual MO ids, they simply refer to the orbitals involved in that spefcific
-  SOMO, for e.g. an excitation of the type [2 2 2 1 1 1 1 0] -> [ 2 2 1 1 1 1 1]
-  implies an excitation from orbital \(3\) to orbital \(8\) which are the real MO ids.
-  However, the prototype only concerns the SOMOs like so [2 1 1 1 1 0] -> [ 1 1 1 1 1 1]
-  therefore \(p,q\) are model space ids \(1,6\).
-
-  #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f
-
-  BEGIN_PROVIDER [ integer*8, NSOMOMax]
- &BEGIN_PROVIDER [ integer*8, NCSFMax]
- &BEGIN_PROVIDER [ integer*8, NMO]
-  implicit none
-  BEGIN_DOC
-  ! Documentation for NSOMOMax
-  ! The maximum number of SOMOs for the current calculation.
-  ! required for the calculation of prototype arrays.
-  END_DOC
-  NSOMOMax = 8
-  NCSFMax = 14 ! TODO: NCSFs for MS=0
-  NMO = NSOMOMax ! TODO: remove this
-  END_PROVIDER
-  #+end_src
-
-  The prototype matrix AIJpqMatrixList has the following dimensions
-  \(\left(NSOMOMax, NSOMOMax, 4, NSOMOMax, NSOMOMax,NCSFMAx,NCSFMax\right)\) where the first two
-  indices represent the somos in \(I,J\) followed by the type of excitation and
-  finally the two model space orbitals \(p,q\).
-
-  The dimensions for each Isomo, Jsomo pair are precalculated and stored in the AIJpqMatrixDimsList
-  variable which is provided here.
-
-
-  #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f
-  BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
-  use cfunctions
-  implicit none
-  BEGIN_DOC
-  ! Documentation for AIJpqMatrixList
-  ! The prototype matrix containing the <I|E_{pq}|J>
-  ! matrices for each I,J somo pair and orb ids.
-  END_DOC
-  integer i,j,k,l
-  integer*8 Isomo, Jsomo
-  Isomo = 0
-  Jsomo = 0
-  integer*8 rows, cols
-  rows = -1
-  cols = -1
-  integer*8 MS
-  MS = 0
-  print *,"NSOMOMax = ",NSOMOMax
-  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
-  do i = 2, NSOMOMax, 2
-     Isomo = ISHFT(1,i)-1
-     do j = i-2,i+2, 2
-        Jsomo = ISHFT(1,j)-1
-        if(j .GT. NSOMOMax .OR. j .LE. 0) then
-           cycle
-        end if
-        do k = 1,NSOMOMax
-           do l = k,NSOMOMax
-              call getApqIJMatrixDims(Isomo,           &
-                   Jsomo, &
-                   MS,                       &
-                   rows,                     &
-                   cols)
-              print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
-              ! i -> j
-              AIJpqMatrixDimsList(i,j,1,k,l,1) = rows
-              AIJpqMatrixDimsList(i,j,1,k,l,2) = cols
-              AIJpqMatrixDimsList(i,j,1,l,k,1) = rows
-              AIJpqMatrixDimsList(i,j,1,l,k,2) = cols
-              ! j -> i
-              AIJpqMatrixDimsList(j,i,1,k,l,1) = rows
-              AIJpqMatrixDimsList(j,i,1,k,l,2) = cols
-              AIJpqMatrixDimsList(j,i,1,l,k,1) = rows
-              AIJpqMatrixDimsList(j,i,1,l,k,2) = cols
-           end do
-        end do
-     end do
-  end do
-  END_PROVIDER
-
-  #+end_src
-
-- [ ] Read the transformation matrix based on the number of SOMOs
-
-  We go through all the possible SOMOs and build the matrix-elements \(<I|E_{pq}|I>\) and
-  store it in the AIJpq container.
-
-  #+begin_src f90 :main no :tangle configurations_sigma_vector.irp.f
-  BEGIN_PROVIDER [ real*8, AIJpqContainer, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,NSOMOMax,NSOMOMax)]
-  use cfunctions
-  implicit none
-  BEGIN_DOC
-  ! Documentation for AIJpqMatrixList
-  ! The prototype matrix containing the <I|E_{pq}|J>
-  ! matrices for each I,J somo pair and orb ids.
-  END_DOC
-  integer i,j,k,l, orbp, orbq, ri, ci
-  orbp = 0
-  orbq = 0
-  integer*8 Isomo, Jsomo
-  Isomo = 0
-  Jsomo = 0
-  integer*8 rows, cols
-  rows = -1
-  cols = -1
-  integer*8 MS
-  MS = 0
-  real*8,dimension(:),allocatable :: meMatrix
-  print *,"NSOMOMax = ",NSOMOMax
-  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
-  do i = 2, NSOMOMax, 2
-     Isomo = ISHFT(1,i)-1
-     do j = i-2,i+2, 2
-        Jsomo = ISHFT(1,j)-1
-        if(j .GT. NSOMOMax .OR. j .LE. 0) then
-           cycle
-        end if
-        do k = 1,NSOMOMax
-           do l = k,NSOMOMax
-              call getApqIJMatrixDims(Isomo,           &
-                   Jsomo, &
-                   MS,                       &
-                   rows,                     &
-                   cols)
-
-              ! allocate matrix
-              allocate(meMatrix(rows*cols))
-
-              ! fill matrix
-              call getApqIJMatrixDriver(Isomo,           &
-                   Jsomo, &
-                   orbp,                     &
-                   orbq,                     &
-                   MS,                       &
-                   NMO,                      &
-                   meMatrix,                 &
-                   rows,                     &
-                   cols)
-             print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
-              ! i -> j
-             do ri = 1,rows
-                 do ci = 1,cols
-                    AIJpqContainer(i,j,1,k,l,ri,ci) = meMatrix(ri*(cols-1) + ci)
-                 end do
-              end do
-           end do
-        end do
-     end do
-  end do
-  END_PROVIDER
-  #+end_src
-
-- [ ] Perform the conversion by matrix-vector BLAS level 2 call

src/determinants/configurations_sigma_vector.irp.f

@@ -1,126 +0,0 @@
-  BEGIN_PROVIDER [ integer*8, NSOMOMax]
- &BEGIN_PROVIDER [ integer*8, NCSFMax]
- &BEGIN_PROVIDER [ integer*8, NMO]
-  implicit none
-  BEGIN_DOC
-  ! Documentation for NSOMOMax
-  ! The maximum number of SOMOs for the current calculation.
-  ! required for the calculation of prototype arrays.
-  END_DOC
-  NSOMOMax = 8
-  NCSFMax = 14 ! TODO: NCSFs for MS=0
-  NMO = NSOMOMax ! TODO: remove this
-  END_PROVIDER
-
-  BEGIN_PROVIDER [ double precision, AIJpqMatrixDimsList, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2)]
-  use cfunctions
-  implicit none
-  BEGIN_DOC
-  ! Documentation for AIJpqMatrixList
-  ! The prototype matrix containing the <I|E_{pq}|J>
-  ! matrices for each I,J somo pair and orb ids.
-  END_DOC
-  integer i,j,k,l
-  integer*8 Isomo, Jsomo
-  Isomo = 0
-  Jsomo = 0
-  integer*8 rows, cols
-  rows = -1
-  cols = -1
-  integer*8 MS
-  MS = 0
-  print *,"NSOMOMax = ",NSOMOMax
-  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
-  do i = 2, NSOMOMax, 2
-     Isomo = ISHFT(1,i)-1
-     do j = i-2,i+2, 2
-        Jsomo = ISHFT(1,j)-1
-        if(j .GT. NSOMOMax .OR. j .LE. 0) then
-           cycle
-        end if
-        do k = 1,NSOMOMax
-           do l = k,NSOMOMax
-              call getApqIJMatrixDims(Isomo,           &
-                   Jsomo, &
-                   MS,                       &
-                   rows,                     &
-                   cols)
-              print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
-              ! i -> j
-              AIJpqMatrixDimsList(i,j,1,k,l,1) = rows
-              AIJpqMatrixDimsList(i,j,1,k,l,2) = cols
-              AIJpqMatrixDimsList(i,j,1,l,k,1) = rows
-              AIJpqMatrixDimsList(i,j,1,l,k,2) = cols
-              ! j -> i
-              AIJpqMatrixDimsList(j,i,1,k,l,1) = rows
-              AIJpqMatrixDimsList(j,i,1,k,l,2) = cols
-              AIJpqMatrixDimsList(j,i,1,l,k,1) = rows
-              AIJpqMatrixDimsList(j,i,1,l,k,2) = cols
-           end do
-        end do
-     end do
-  end do
-  END_PROVIDER
-
-  BEGIN_PROVIDER [ real*8, AIJpqContainer, (NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,NSOMOMax,NSOMOMax)]
-  use cfunctions
-  implicit none
-  BEGIN_DOC
-  ! Documentation for AIJpqMatrixList
-  ! The prototype matrix containing the <I|E_{pq}|J>
-  ! matrices for each I,J somo pair and orb ids.
-  END_DOC
-  integer i,j,k,l, orbp, orbq, ri, ci
-  orbp = 0
-  orbq = 0
-  integer*8 Isomo, Jsomo
-  Isomo = 0
-  Jsomo = 0
-  integer*8 rows, cols
-  rows = -1
-  cols = -1
-  integer*8 MS
-  MS = 0
-  real*8,dimension(:),allocatable :: meMatrix
-  print *,"NSOMOMax = ",NSOMOMax
-  !allocate(AIJpqMatrixDimsList(NSOMOMax,NSOMOMax,4,NSOMOMax,NSOMOMax,2))
-  do i = 2, NSOMOMax, 2
-     Isomo = ISHFT(1,i)-1
-     do j = i-2,i+2, 2
-        Jsomo = ISHFT(1,j)-1
-        if(j .GT. NSOMOMax .OR. j .LE. 0) then
-           cycle
-        end if
-        do k = 1,NSOMOMax
-           do l = k,NSOMOMax
-              call getApqIJMatrixDims(Isomo,           &
-                   Jsomo, &
-                   MS,                       &
-                   rows,                     &
-                   cols)
-
-              ! allocate matrix
-              allocate(meMatrix(rows*cols))
-
-              ! fill matrix
-              call getApqIJMatrixDriver(Isomo,           &
-                   Jsomo, &
-                   orbp,                     &
-                   orbq,                     &
-                   MS,                       &
-                   NMO,                      &
-                   meMatrix,                 &
-                   rows,                     &
-                   cols)
-             print *, i,j,k,l,">",Isomo,Jsomo,">",rows, cols
-              ! i -> j
-             do ri = 1,rows
-                 do ci = 1,cols
-                    AIJpqContainer(i,j,1,k,l,ri,ci) = meMatrix(ri*(cols-1) + ci)
-                 end do
-              end do
-           end do
-        end do
-     end do
-  end do
-  END_PROVIDER