started complex h_apply

src/determinants/h_apply.irp.f

@@ -6,6 +6,7 @@ type H_apply_buffer_type
   integer                        :: sze
   integer(bit_kind), pointer     :: det(:,:,:)
   double precision , pointer     :: coef(:,:)
+  complex*16 , pointer           :: coef_complex(:,:)
   double precision , pointer     :: e2(:,:)
 end type H_apply_buffer_type
 
@@ -32,11 +33,16 @@ type(H_apply_buffer_type), pointer :: H_apply_buffer(:)
     H_apply_buffer(iproc)%sze = sze
     allocate (                                                       &
         H_apply_buffer(iproc)%det(N_int,2,sze),                      &
-        H_apply_buffer(iproc)%coef(sze,N_states),                    &
         H_apply_buffer(iproc)%e2(sze,N_states)                       &
         )
-    H_apply_buffer(iproc)%det  = 0_bit_kind
+    if (is_complex) then
+      allocate(H_apply_buffer(iproc)%coef_complex(sze,N_states))
+      H_apply_buffer(iproc)%coef_complex = (0.d0,0.d0)
+    else
+      allocate(H_apply_buffer(iproc)%coef(sze,N_states))
       H_apply_buffer(iproc)%coef = 0.d0
+    endif
+    H_apply_buffer(iproc)%det  = 0_bit_kind
     H_apply_buffer(iproc)%e2   = 0.d0
     call omp_init_lock(H_apply_buffer_lock(1,iproc))
     !$OMP END PARALLEL
@@ -59,6 +65,7 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   integer, intent(in)            :: new_size, iproc
   integer(bit_kind), pointer     :: buffer_det(:,:,:)
   double precision,  pointer     :: buffer_coef(:,:)
+  complex*16,  pointer           :: buffer_coef_complex(:,:)
   double precision,  pointer     :: buffer_e2(:,:)
   integer                        :: i,j,k
   integer                        :: Ndet
@@ -74,9 +81,14 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   ASSERT (iproc < nproc)
 
   allocate ( buffer_det(N_int,2,new_size),                           &
-      buffer_coef(new_size,N_states),                                &
       buffer_e2(new_size,N_states) )
+  if (is_complex) then
+    allocate(buffer_coef_complex(new_size,N_states))
+    buffer_coef_complex = (0.d0,0.d0)
+  else
+    allocate(buffer_coef(new_size,N_states))
     buffer_coef = 0.d0
+  endif
   buffer_e2 = 0.d0
   do i=1,min(new_size,H_apply_buffer(iproc)%N_det)
     do k=1,N_int
@@ -89,6 +101,15 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   deallocate(H_apply_buffer(iproc)%det)
   H_apply_buffer(iproc)%det => buffer_det
 
+  if (is_complex) then
+    do k=1,N_states
+      do i=1,min(new_size,H_apply_buffer(iproc)%N_det)
+        buffer_coef_complex(i,k) = H_apply_buffer(iproc)%coef_complex(i,k)
+      enddo
+    enddo
+    deallocate(H_apply_buffer(iproc)%coef_complex)
+    H_apply_buffer(iproc)%coef_complex => buffer_coef_complex
+  else
   do k=1,N_states
     do i=1,min(new_size,H_apply_buffer(iproc)%N_det)
       buffer_coef(i,k) = H_apply_buffer(iproc)%coef(i,k)
@@ -96,6 +117,7 @@ subroutine resize_H_apply_buffer(new_size,iproc)
   enddo
   deallocate(H_apply_buffer(iproc)%coef)
   H_apply_buffer(iproc)%coef => buffer_coef
+  endif
 
   do k=1,N_states
     do i=1,min(new_size,H_apply_buffer(iproc)%N_det)
@@ -119,6 +141,7 @@ subroutine copy_H_apply_buffer_to_wf
   END_DOC
   integer(bit_kind), allocatable :: buffer_det(:,:,:)
   double precision, allocatable  :: buffer_coef(:,:)
+  complex*16, allocatable        :: buffer_coef_complex(:,:)
   integer                        :: i,j,k
   integer                        :: N_det_old
 
@@ -128,7 +151,12 @@ subroutine copy_H_apply_buffer_to_wf
   ASSERT (N_int > 0)
   ASSERT (N_det > 0)
 
-  allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) )
+  allocate ( buffer_det(N_int,2,N_det))
+  if (is_complex) then
+    allocate(buffer_coef_complex(N_det,N_states))
+  else
+    allocate(buffer_coef(N_det,N_states))
+  endif
 
   ! Backup determinants
   j=0
@@ -142,6 +170,17 @@ subroutine copy_H_apply_buffer_to_wf
   N_det_old = j
 
   ! Backup coefficients
+  if (is_complex) then
+    do k=1,N_states
+      j=0
+      do i=1,N_det
+        if (pruned(i)) cycle  ! Pruned determinants
+        j += 1
+        buffer_coef_complex(j,k) = psi_coef_complex(i,k)
+      enddo
+      ASSERT ( j == N_det_old )
+    enddo
+  else
   do k=1,N_states
     j=0
     do i=1,N_det
@@ -151,6 +190,7 @@ subroutine copy_H_apply_buffer_to_wf
     enddo
     ASSERT ( j == N_det_old )
   enddo
+  endif
 
   ! Update N_det
   N_det = N_det_old
@@ -170,14 +210,57 @@ subroutine copy_H_apply_buffer_to_wf
     ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
     ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
   enddo
+  if (is_complex) then
+    do k=1,N_states
+      do i=1,N_det_old
+        psi_coef_complex(i,k) = buffer_coef_complex(i,k)
+      enddo
+    enddo
+  else
   do k=1,N_states
     do i=1,N_det_old
       psi_coef(i,k) = buffer_coef(i,k)
     enddo
   enddo
+  endif
 
   ! Copy new buffers
 
+  if (is_complex) then
+    !$OMP PARALLEL DEFAULT(SHARED)                                     &
+        !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old)                   &
+        !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef_complex,N_states,psi_det_size)
+    j=0
+    !$ j=omp_get_thread_num()
+    do k=0,j-1
+      N_det_old += H_apply_buffer(k)%N_det
+    enddo
+    do i=1,H_apply_buffer(j)%N_det
+      do k=1,N_int
+        psi_det(k,1,i+N_det_old) = H_apply_buffer(j)%det(k,1,i)
+        psi_det(k,2,i+N_det_old) = H_apply_buffer(j)%det(k,2,i)
+      enddo
+      ASSERT (sum(popcnt(psi_det(:,1,i+N_det_old))) == elec_alpha_num)
+      ASSERT (sum(popcnt(psi_det(:,2,i+N_det_old))) == elec_beta_num )
+    enddo
+    do k=1,N_states
+      do i=1,H_apply_buffer(j)%N_det
+        psi_coef_complex(i+N_det_old,k) = H_apply_buffer(j)%coef_complex(i,k)
+      enddo
+    enddo
+    !$OMP BARRIER
+    H_apply_buffer(j)%N_det = 0
+    !$OMP END PARALLEL
+    SOFT_TOUCH N_det psi_det psi_coef_complex
+
+    logical :: found_duplicates
+    call remove_duplicates_in_psi_det(found_duplicates)
+    do k=1,N_states
+      call normalize(psi_coef_complex(1,k),N_det)
+    enddo
+    SOFT_TOUCH N_det psi_det psi_coef_complex
+  else
+
   !$OMP PARALLEL DEFAULT(SHARED)                                     &
       !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old)                   &
       !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states,psi_det_size)
@@ -211,6 +294,7 @@ subroutine copy_H_apply_buffer_to_wf
   enddo
   SOFT_TOUCH N_det psi_det psi_coef
   
+  endif
 end
 
 subroutine remove_duplicates_in_psi_det(found_duplicates)
@@ -275,6 +359,29 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
   !$OMP END DO
   !$OMP END PARALLEL
   
+  if (is_complex) then
+    if (found_duplicates) then
+      k=0
+      do i=1,N_det
+        if (.not.duplicate(i)) then
+          k += 1
+          psi_det(:,:,k) = psi_det_sorted_bit (:,:,i)
+          psi_coef_complex(k,:)  = psi_coef_sorted_bit_complex(i,:)
+        else
+          if (sum(cdabs(psi_coef_sorted_bit_complex(i,:))) /= 0.d0 ) then
+            psi_coef_complex(k,:)  = psi_coef_sorted_bit_complex(i,:)
+          endif
+        endif
+      enddo
+      N_det = k
+      psi_det_sorted_bit(:,:,1:N_det) = psi_det(:,:,1:N_det)
+      psi_coef_sorted_bit_complex(1:N_det,:) = psi_coef_complex(1:N_det,:)
+      TOUCH N_det psi_det psi_coef_complex psi_det_sorted_bit psi_coef_sorted_bit_complex c0_weight
+    endif
+    psi_det = psi_det_sorted
+    psi_coef_complex = psi_coef_sorted_complex
+    SOFT_TOUCH psi_det psi_coef_complex psi_det_sorted_bit psi_coef_sorted_bit_complex
+  else
   if (found_duplicates) then
     k=0
     do i=1,N_det
@@ -296,6 +403,7 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
   psi_det = psi_det_sorted
   psi_coef = psi_coef_sorted
   SOFT_TOUCH psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit
+  endif
   deallocate (duplicate,bit_tmp)
 end
 
@@ -329,11 +437,19 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num)
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num)
   enddo
+  if (is_complex) then
+    do j=1,N_states
+      do i=1,N_selected
+        H_apply_buffer(iproc)%coef_complex(i+H_apply_buffer(iproc)%N_det,j) = (0.d0,0.d0)
+      enddo
+    enddo
+  else
   do j=1,N_states
     do i=1,N_selected
       H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = 0.d0
     enddo
   enddo
+  endif
   H_apply_buffer(iproc)%N_det = new_size
   do i=1,H_apply_buffer(iproc)%N_det
     ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
@@ -343,6 +459,11 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
 end
 
 subroutine push_pt2(zmq_socket_push,pt2,norm_pert,H_pert_diag,i_generator,N_st,task_id)
+  !todo: modify for complex
+  if (is_complex) then
+    print*,irp_here,' not implemented for complex'
+    stop -1
+  endif
   use f77_zmq
   implicit none
   BEGIN_DOC
@@ -404,6 +525,11 @@ IRP_ENDIF
 end
 
 subroutine pull_pt2(zmq_socket_pull,pt2,norm_pert,H_pert_diag,i_generator,N_st,n,task_id)
+  !todo: modify for complex
+  if (is_complex) then
+    print*,irp_here,' not implemented for complex'
+    stop -1
+  endif
   use f77_zmq
   implicit none
   BEGIN_DOC

src/determinants/h_apply_nozmq.template.f

@@ -17,8 +17,11 @@ subroutine $subroutine($params_main)
   double precision, allocatable  :: fock_diag_tmp(:,:)
 
   $initialization
+  if (is_complex) then
+    PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators_complex
+  else
   PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators
-
+  endif
 
   call wall_time(wall_0)
 

src/determinants/occ_pattern.irp.f

@@ -513,7 +513,7 @@ subroutine make_s2_eigenfunction
       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)
+        call fill_h_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread)
         N_det_new = 0
       endif
     enddo
@@ -528,7 +528,7 @@ subroutine make_s2_eigenfunction
   !$OMP END PARALLEL
 
   if (update) then
-    call copy_H_apply_buffer_to_wf
+    call copy_h_apply_buffer_to_wf
     if (is_complex) then
       TOUCH N_det psi_coef_complex psi_det psi_occ_pattern N_occ_pattern
     else

src/generators_full/generators.irp.f

@@ -23,19 +23,34 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
 END_PROVIDER
 
 BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! For Single reference wave functions, the generator is the
  ! Hartree-Fock determinant
  END_DOC
  psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
- psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
+END_PROVIDER
+
+BEGIN_PROVIDER [ complex*16, psi_coef_generators_complex, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_coef_generators_complex(1:N_det,1:N_states) = psi_coef_sorted_complex(1:N_det,1:N_states)
+END_PROVIDER
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
 &BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order,     (psi_det_size)  ]
 
  implicit none
@@ -44,10 +59,26 @@ END_PROVIDER
  ! Hartree-Fock determinant
  END_DOC
  psi_det_sorted_gen = psi_det_sorted
- psi_coef_sorted_gen = psi_coef_sorted
  psi_det_sorted_gen_order = psi_det_sorted_order
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_coef_sorted_gen = psi_coef_sorted
+END_PROVIDER
+
+BEGIN_PROVIDER [ complex*16, psi_coef_sorted_gen_complex, (psi_det_size,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! For Single reference wave functions, the generator is the
+ ! Hartree-Fock determinant
+ END_DOC
+ psi_coef_sorted_gen_complex = psi_coef_sorted_complex
+END_PROVIDER
 
 BEGIN_PROVIDER [integer, degree_max_generators]
  implicit none

src/utils_complex/qp2-pbc-diff.txt

@@ -22,6 +22,11 @@ determinants:
   (done) filter_connected.irp.f
   (done) fock_diag.irp.f
   (****) h_apply.irp.f
+           added coef_complex to h_apply_buffer_type
+           either coef or coef_complex will remain unallocated
+           (if this causes problems (it shouldn't), maybe just allocate unused one with size 1?)
+           check {push,pull}_pt2
+             pt2, norm_pert, h_pert_diag types? (should be real? documentation?)
   (****) h_apply_nozmq.template.f
   (****) h_apply.template.f
   (****) h_apply_zmq.template.f