Working on periodic

src/ao_one_e_ints/ao_overlap.irp.f

@@ -70,6 +70,29 @@
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_overlap_imag, (ao_num, ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Imaginary part of the overlap
+ END_DOC
+ ao_overlap_imag = 0.d0
+END_PROVIDER
+
+BEGIN_PROVIDER [ complex*16, ao_overlap_complex, (ao_num, ao_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Overlap for complex AOs
+  END_DOC
+  integer                        :: i,j
+  do j=1,ao_num
+    do i=1,ao_num
+      ao_overlap_complex(i,j) = dcmplx( ao_overlap(i,j), ao_overlap_imag(i,j) )
+    enddo
+  enddo
+END_PROVIDER
+
+
+
 
 BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
   implicit none
@@ -86,6 +109,13 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
   double precision :: A_center(3), B_center(3)
   integer :: power_A(3), power_B(3)
   double precision :: lower_exp_val, dx
+  if (periodic) then
+    do j=1,ao_num
+      do i= 1,ao_num
+        ao_overlap_abs(i,j)= cdabs(ao_overlap_complex(i,j))
+      enddo
+    enddo
+  else
     dim1=100
     lower_exp_val = 40.d0
     !$OMP PARALLEL DO SCHEDULE(GUIDED)                               &
@@ -124,6 +154,7 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
       enddo
     enddo
     !$OMP END PARALLEL DO
+  endif
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
@@ -134,6 +165,15 @@ BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
  call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1))
 END_PROVIDER
 
+BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]
+ implicit none
+ BEGIN_DOC
+! Inverse of the overlap matrix
+ END_DOC
+ call get_pseudo_inverse_complex(ao_overlap_complex, &
+    size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1))
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
 
   BEGIN_DOC

src/ao_two_e_ints/map_integrals.irp.f

@@ -21,7 +21,7 @@ subroutine two_e_integrals_index(i,j,k,l,i1)
   implicit none
   BEGIN_DOC
 ! Gives a unique index for i,j,k,l using permtuation symmetry.
-! i <-> k, j <-> l, and (i,k) <-> (j,l)
+! i <-> k, j <-> l, and (i,k) <-> (j,l) for non-periodic systems
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer(key_kind), intent(out) :: i1
@@ -37,6 +37,8 @@ subroutine two_e_integrals_index(i,j,k,l,i1)
   i1 = i1+shiftr(i2*i2-i2,1)
 end
 
+
+
 subroutine two_e_integrals_index_reverse(i,j,k,l,i1)
   use map_module
   implicit none
@@ -126,6 +128,155 @@ subroutine two_e_integrals_index_reverse(i,j,k,l,i1)
 
 end
 
+
+
+subroutine ao_idx2_sq(i,j,ij)
+  implicit none
+  integer, intent(in)  :: i,j
+  integer, intent(out) :: ij
+  if (i<j) then
+    ij=(j-1)*(j-1)+2*i-mod(j+1,2)
+  else if (i>j) then
+    ij=(i-1)*(i-1)+2*j-mod(i,2)
+  else
+    ij=i*i
+  endif
+end
+  
+subroutine idx2_tri_int(i,j,ij)
+  implicit none
+  integer, intent(in)  :: i,j
+  integer, intent(out) :: ij
+  integer :: p,q
+  p = max(i,j)
+  q = min(i,j)
+  ij = q+ishft(p*p-p,-1)
+end
+  
+subroutine ao_idx2_tri_key(i,j,ij)
+  use map_module
+  implicit none
+  integer, intent(in)  :: i,j
+  integer(key_kind), intent(out) :: ij
+  integer(key_kind) :: p,q
+  p = max(i,j)
+  q = min(i,j)
+  ij = q+ishft(p*p-p,-1)
+end
+  
+subroutine two_e_integrals_index_2fold(i,j,k,l,i1)                                           
+  use map_module
+  implicit none
+  integer, intent(in)            :: i,j,k,l
+  integer(key_kind), intent(out) :: i1
+  integer                        :: ik,jl
+
+  call ao_idx2_sq(i,k,ik)
+  call ao_idx2_sq(j,l,jl)
+  call ao_idx2_tri_key(ik,jl,i1)
+end
+
+subroutine ao_idx2_sq_rev(i,k,ik)                                                                
+  BEGIN_DOC
+  ! reverse square compound index
+  END_DOC
+!  p = ceiling(dsqrt(dble(ik)))
+!  q = ceiling(0.5d0*(dble(ik)-dble((p-1)*(p-1))))
+!  if (mod(ik,2)==0) then
+!    k=p
+!    i=q
+!  else
+!    i=p
+!    k=q
+!  endif
+  integer, intent(in)           :: ik
+  integer, intent(out)          :: i,k
+  integer                       :: pq(0:1),i1,i2
+  pq(0) = ceiling(dsqrt(dble(ik)))
+  pq(1) = ceiling(0.5d0*(dble(ik)-dble((pq(0)-1)*(pq(0)-1))))
+  i1=mod(ik,2)
+  i2=mod(ik+1,2)
+
+  k=pq(i1)
+  i=pq(i2)
+end
+
+subroutine ao_idx2_tri_rev_key(i,k,ik)
+  use map_module
+  BEGIN_DOC
+  !return i<=k
+  END_DOC
+  integer(key_kind), intent(in) :: ik
+  integer, intent(out)          :: i,k
+  integer(key_kind) :: tmp_k
+  k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
+  tmp_k = k
+  i = int(ik - ishft(tmp_k*tmp_k-tmp_k,-1))
+end
+
+subroutine idx2_tri_rev_int(i,k,ik)
+  BEGIN_DOC
+  !return i<=k
+  END_DOC
+  integer, intent(in)           :: ik
+  integer, intent(out)          :: i,k
+  k = ceiling(0.5d0*(dsqrt(8.d0*dble(ik)+1.d0)-1.d0))
+  i = int(ik - ishft(k*k-k,-1))
+end
+
+subroutine two_e_integrals_index_reverse_2fold(i,j,k,l,i1)
+  use map_module
+  implicit none
+  integer, intent(out)           :: i(2),j(2),k(2),l(2)
+  integer(key_kind), intent(in)  :: i1
+  integer(key_kind)              :: i0
+  integer                        :: i2,i3
+  i = 0
+  call ao_idx2_tri_rev_key(i3,i2,i1)
+
+  call ao_idx2_sq_rev(j(1),l(1),i2)
+  call ao_idx2_sq_rev(i(1),k(1),i3)
+
+              !ijkl
+  i(2) = j(1) !jilk
+  j(2) = i(1)
+  k(2) = l(1)
+  l(2) = k(1)
+
+!  i(3) = k(1) !klij   complex conjugate
+!  j(3) = l(1)
+!  k(3) = i(1)
+!  l(3) = j(1)
+!
+!  i(4) = l(1) !lkji   complex conjugate
+!  j(4) = k(1)
+!  k(4) = j(1)
+!  l(4) = i(1)
+
+  integer :: ii
+  if ( (i(1)==i(2)).and. &
+       (j(1)==j(2)).and. &
+       (k(1)==k(2)).and. &
+       (l(1)==l(2)) ) then
+    i(2) = 0
+  endif
+! This has been tested with up to 1000 AOs, and all the reverse indices are
+! correct ! We can remove the test
+!  do ii=1,2
+!    if (i(ii) /= 0) then
+!      call two_e_integrals_index_2fold(i(ii),j(ii),k(ii),l(ii),i0)
+!      if (i1 /= i0) then
+!        print *,  i1, i0
+!        print *,  i(ii), j(ii), k(ii), l(ii)
+!        stop 'two_e_integrals_index_reverse_2fold failed'
+!      endif
+!    endif
+!  enddo
+end
+
+
+
+
  BEGIN_PROVIDER [ integer, ao_integrals_cache_min ]
 &BEGIN_PROVIDER [ integer, ao_integrals_cache_max ]
  implicit none
@@ -144,8 +295,11 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
  END_DOC
  PROVIDE ao_two_e_integrals_in_map
  integer                        :: i,j,k,l,ii
- integer(key_kind)              :: idx
+ integer(key_kind)              :: idx, idx2
  real(integral_kind)            :: integral
+ real(integral_kind)            :: tmp_re, tmp_im
+ integer(key_kind)              :: idx_re,idx_im
+
   !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
   do l=ao_integrals_cache_min,ao_integrals_cache_max
     do k=ao_integrals_cache_min,ao_integrals_cache_max
@@ -165,7 +319,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_cache, (0:64*64*64*64) ]
     enddo
   enddo
   !$OMP END PARALLEL DO
-
 END_PROVIDER
 
 
@@ -207,6 +360,113 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result)
   result = tmp
 end
 
+BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ]
+ implicit none
+ BEGIN_DOC
+ ! Cache of AO integrals for fast access
+ END_DOC
+ PROVIDE ao_two_e_integrals_in_map
+ integer                        :: i,j,k,l,ii
+ integer(key_kind)              :: idx1, idx2
+ real(integral_kind)            :: tmp_re, tmp_im
+ integer(key_kind)              :: idx_re,idx_im
+ complex(integral_kind)         :: integral
+
+
+ !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx1,idx2,tmp_re,tmp_im,idx_re,idx_im,ii,integral)
+ do l=ao_integrals_cache_min,ao_integrals_cache_max
+   do k=ao_integrals_cache_min,ao_integrals_cache_max
+     do j=ao_integrals_cache_min,ao_integrals_cache_max
+       do i=ao_integrals_cache_min,ao_integrals_cache_max
+         !DIR$ FORCEINLINE
+         call two_e_integrals_index_2fold(i,j,k,l,idx1)
+         !DIR$ FORCEINLINE
+         call two_e_integrals_index_2fold(k,l,i,j,idx2)
+         idx_re = min(idx1,idx2)
+         idx_im = max(idx1,idx2)
+         !DIR$ FORCEINLINE
+         call map_get(ao_integrals_map,idx_re,tmp_re)
+         if (idx_re /= idx_im) then
+           call map_get(ao_integrals_map,idx_im,tmp_im)
+           if (idx1 < idx2) then
+             integral = cmplx(tmp_re,tmp_im)
+           else
+             integral = cmplx(tmp_re,-tmp_im)
+           endif
+         else
+           tmp_im = 0.d0
+           integral = cmplx(tmp_re,tmp_im)
+         endif
+         
+         ii = l-ao_integrals_cache_min
+         ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
+         ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
+         ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
+         ao_integrals_cache_periodic(ii) = integral
+       enddo
+     enddo
+   enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
+
+complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map) result(result)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets one AO bi-electronic integral from the AO map
+  END_DOC
+  integer, intent(in)            :: i,j,k,l
+  integer(key_kind)              :: idx1,idx2
+  real(integral_kind)            :: tmp_re, tmp_im
+  integer(key_kind)              :: idx_re,idx_im
+  type(map_type), intent(inout)  :: map
+  integer                        :: ii
+  complex(integral_kind)         :: tmp
+  PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min
+  !DIR$ FORCEINLINE
+  if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then
+    tmp = (0.d0,0.d0)
+  else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then
+    tmp = (0.d0,0.d0)
+  else
+    ii = l-ao_integrals_cache_min
+    ii = ior(ii, k-ao_integrals_cache_min)
+    ii = ior(ii, j-ao_integrals_cache_min)
+    ii = ior(ii, i-ao_integrals_cache_min)
+    if (iand(ii, -64) /= 0) then
+         !DIR$ FORCEINLINE
+         call two_e_integrals_index_2fold(i,j,k,l,idx1)
+         !DIR$ FORCEINLINE
+         call two_e_integrals_index_2fold(k,l,i,j,idx2)
+         idx_re = min(idx1,idx2)
+         idx_im = max(idx1,idx2)
+         !DIR$ FORCEINLINE
+         call map_get(ao_integrals_map,idx_re,tmp_re)
+         if (idx_re /= idx_im) then
+           call map_get(ao_integrals_map,idx_im,tmp_im)
+           if (idx1 < idx2) then
+             tmp = cmplx(tmp_re,tmp_im)
+           else
+             tmp = cmplx(tmp_re,-tmp_im)
+           endif
+         else
+           tmp_im = 0.d0
+           tmp = cmplx(tmp_re,tmp_im)
+         endif
+    else
+      ii = l-ao_integrals_cache_min
+      ii = ior( shiftl(ii,6), k-ao_integrals_cache_min)
+      ii = ior( shiftl(ii,6), j-ao_integrals_cache_min)
+      ii = ior( shiftl(ii,6), i-ao_integrals_cache_min)
+      tmp = ao_integrals_cache_periodic(ii)
+    endif
+    result = tmp
+  endif
+end
+
 
 subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   use map_module
@@ -237,6 +497,36 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
 
 end
 
+
+subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val)
+  use map_module
+  BEGIN_DOC
+  ! Gets multiple AO bi-electronic integral from the AO map .
+  ! All i are retrieved for j,k,l fixed.
+  ! physicist convention : <ij|kl> 
+  END_DOC
+  implicit none
+  integer, intent(in)            :: j,k,l, sze
+  complex(integral_kind), intent(out) :: out_val(sze)
+
+  integer                        :: i
+  integer(key_kind)              :: hash
+  double precision               :: thresh
+  PROVIDE ao_two_e_integrals_in_map ao_integrals_map
+  thresh = ao_integrals_threshold
+
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  double precision :: get_ao_two_e_integral
+  do i=1,sze
+    out_val(i) = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
+  enddo
+
+end
+
 subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_zero_int)
   use map_module
   implicit none
@@ -407,81 +697,81 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
 end
 
 
-subroutine dump_ao_integrals(filename)
+!subroutine dump_ao_integrals(filename)
-  use map_module
+!  use map_module
-  implicit none
+!  implicit none
-  BEGIN_DOC
+!  BEGIN_DOC
-  ! Save to disk the |AO| integrals
+!  ! Save to disk the |AO| integrals
-  END_DOC
+!  END_DOC
-  character*(*), intent(in)      :: filename
+!  character*(*), intent(in)      :: filename
-  integer(cache_key_kind), pointer :: key(:)
+!  integer(cache_key_kind), pointer :: key(:)
-  real(integral_kind), pointer   :: val(:)
+!  real(integral_kind), pointer   :: val(:)
-  integer*8                      :: i,j, n
+!  integer*8                      :: i,j, n
-  if (.not.mpi_master) then
+!  if (.not.mpi_master) then
-    return
+!    return
-  endif
+!  endif
-  call ezfio_set_work_empty(.False.)
+!  call ezfio_set_work_empty(.False.)
-  open(unit=66,file=filename,FORM='unformatted')
+!  open(unit=66,file=filename,FORM='unformatted')
-  write(66) integral_kind, key_kind
+!  write(66) integral_kind, key_kind
-  write(66) ao_integrals_map%sorted, ao_integrals_map%map_size,    &
+!  write(66) ao_integrals_map%sorted, ao_integrals_map%map_size,    &
-      ao_integrals_map%n_elements
+!      ao_integrals_map%n_elements
-  do i=0_8,ao_integrals_map%map_size
+!  do i=0_8,ao_integrals_map%map_size
-    write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,&
+!    write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,&
-        ao_integrals_map%map(i)%n_elements
+!        ao_integrals_map%map(i)%n_elements
-  enddo
+!  enddo
-  do i=0_8,ao_integrals_map%map_size
+!  do i=0_8,ao_integrals_map%map_size
-    key => ao_integrals_map%map(i)%key
+!    key => ao_integrals_map%map(i)%key
-    val => ao_integrals_map%map(i)%value
+!    val => ao_integrals_map%map(i)%value
-    n = ao_integrals_map%map(i)%n_elements
+!    n = ao_integrals_map%map(i)%n_elements
-    write(66) (key(j), j=1,n), (val(j), j=1,n)
+!    write(66) (key(j), j=1,n), (val(j), j=1,n)
-  enddo
+!  enddo
-  close(66)
+!  close(66)
-
+!
-end
+!end
 
 
-integer function load_ao_integrals(filename)
+!integer function load_ao_integrals(filename)
-  implicit none
+!  implicit none
-  BEGIN_DOC
+!  BEGIN_DOC
-  ! Read from disk the |AO| integrals
+!  ! Read from disk the |AO| integrals
-  END_DOC
+!  END_DOC
-  character*(*), intent(in)      :: filename
+!  character*(*), intent(in)      :: filename
-  integer*8                      :: i
+!  integer*8                      :: i
-  integer(cache_key_kind), pointer :: key(:)
+!  integer(cache_key_kind), pointer :: key(:)
-  real(integral_kind), pointer   :: val(:)
+!  real(integral_kind), pointer   :: val(:)
-  integer                        :: iknd, kknd
+!  integer                        :: iknd, kknd
-  integer*8                      :: n, j
+!  integer*8                      :: n, j
-  load_ao_integrals = 1
+!  load_ao_integrals = 1
-  open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
+!  open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
-  read(66,err=98,end=98) iknd, kknd
+!  read(66,err=98,end=98) iknd, kknd
-  if (iknd /= integral_kind) then
+!  if (iknd /= integral_kind) then
-    print *,  'Wrong integrals kind in file :', iknd
+!    print *,  'Wrong integrals kind in file :', iknd
-    stop 1
+!    stop 1
-  endif
+!  endif
-  if (kknd /= key_kind) then
+!  if (kknd /= key_kind) then
-    print *,  'Wrong key kind in file :', kknd
+!    print *,  'Wrong key kind in file :', kknd
-    stop 1
+!    stop 1
-  endif
+!  endif
-  read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,&
+!  read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,&
-      ao_integrals_map%n_elements
+!      ao_integrals_map%n_elements
-  do i=0_8, ao_integrals_map%map_size
+!  do i=0_8, ao_integrals_map%map_size
-    read(66,err=99,end=99) ao_integrals_map%map(i)%sorted,          &
+!    read(66,err=99,end=99) ao_integrals_map%map(i)%sorted,          &
-        ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements
+!        ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements
-    call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size)
+!    call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size)
-  enddo
+!  enddo
-  do i=0_8, ao_integrals_map%map_size
+!  do i=0_8, ao_integrals_map%map_size
-    key => ao_integrals_map%map(i)%key
+!    key => ao_integrals_map%map(i)%key
-    val => ao_integrals_map%map(i)%value
+!    val => ao_integrals_map%map(i)%value
-    n = ao_integrals_map%map(i)%n_elements
+!    n = ao_integrals_map%map(i)%n_elements
-    read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
+!    read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
-  enddo
+!  enddo
-  call map_sort(ao_integrals_map)
+!  call map_sort(ao_integrals_map)
-  load_ao_integrals = 0
+!  load_ao_integrals = 0
-  return
+!  return
-  99 continue
+!  99 continue
-  call map_deinit(ao_integrals_map)
+!  call map_deinit(ao_integrals_map)
-  98 continue
+!  98 continue
-  stop 'Problem reading ao_integrals_map file in work/'
+!  stop 'Problem reading ao_integrals_map file in work/'
-
+!
-end
+!end
-
+!

src/utils/linear_algebra.irp.f

@@ -43,6 +43,690 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
 end
 
 
+subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Compute A = U.D.Vt
+  !   
+  ! LDx : leftmost dimension of x
+  !
+  ! Dimension of A is m x n
+  ! A,U,Vt are complex*16
+  ! D is double precision
+  END_DOC
+  
+  integer, intent(in)             :: LDA, LDU, LDVt, m, n
+  complex*16, intent(in)          :: A(LDA,n)
+  complex*16, intent(out)         :: U(LDU,m)
+  complex*16, intent(out)         :: Vt(LDVt,n)
+  double precision,intent(out)    :: D(min(m,n))
+  complex*16,allocatable          :: work(:)
+  double precision,allocatable    :: rwork(:)
+  integer                         :: info, lwork, i, j, k, lrwork
+  
+  complex*16,allocatable          :: A_tmp(:,:)
+  allocate (A_tmp(LDA,n)) 
+  A_tmp = A
+  lrwork = 5*min(m,n)
+  
+  ! Find optimal size for temp arrays
+  allocate(work(1),rwork(lrwork))
+  lwork = -1
+  call zgesvd('A','A', m, n, A_tmp, LDA,                             &
+      D, U, LDU, Vt, LDVt, work, lwork, rwork, info)
+  lwork = int(work(1))
+  deallocate(work)
+  
+  allocate(work(lwork))
+  call zgesvd('A','A', m, n, A_tmp, LDA,                             &
+      D, U, LDU, Vt, LDVt, work, lwork, rwork, info)
+  deallocate(work,rwork,A_tmp)
+  
+  if (info /= 0) then
+    print *,  info, ': SVD failed'
+    stop
+  endif
+      
+end
+  
+subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m)
+  implicit none
+  BEGIN_DOC
+  ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization.
+  !
+  ! overlap : overlap matrix 
+  !
+  ! LDA : leftmost dimension of overlap array
+  !
+  ! N : Overlap matrix is NxN (array is (LDA,N) )
+  !
+  ! C : Coefficients of the vectors to orthogonalize. On exit,
+  !     orthogonal vectors
+  !
+  ! LDC : leftmost dimension of C
+  !
+  ! m : Coefficients matrix is MxN, ( array is (LDC,N) )
+  !
+  END_DOC
+  
+  integer, intent(in)            :: lda, ldc, n
+  integer, intent(out)           :: m
+  complex*16, intent(in)   :: overlap(lda,n)
+  complex*16, intent(inout) :: C(ldc,n)
+  complex*16, allocatable  :: U(:,:)
+  complex*16, allocatable  :: Vt(:,:)
+  double precision, allocatable  :: D(:)
+  complex*16, allocatable  :: S(:,:)
+  !DIR$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
+  integer                        :: info, i, j
+  
+  if (n < 2) then
+    return
+  endif
+  
+  allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n))
+  
+  call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n)
+  
+  D(:) = dsqrt(D(:))
+  m=n
+  do i=1,n
+    if ( D(i) >= 1.d-6 ) then
+      D(i) = 1.d0/D(i)
+    else
+      m = i-1
+      print *,  'Removed Linear dependencies below:', 1.d0/D(m)
+      exit
+    endif
+  enddo
+  do i=m+1,n 
+    D(i) = 0.d0
+  enddo                                                                                       
+      
+  do i=1,m
+    if ( D(i) >= 1.d5 ) then
+      print *,  'Warning: Basis set may have linear dependence problems'
+    endif
+  enddo
+  
+  do j=1,n
+    do i=1,n
+      S(i,j) = U(i,j)*D(j)
+    enddo
+  enddo
+  
+  do j=1,n
+    do i=1,n
+      U(i,j) = C(i,j)
+    enddo
+  enddo
+  
+  call zgemm('N','N',n,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1))
+  deallocate (U, Vt, D, S)                                                                    
+    
+end   
+    
+  
+subroutine ortho_qr_complex(A,LDA,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Orthogonalization using Q.R factorization
+  !   
+  ! A : matrix to orthogonalize
+  !
+  ! LDA : leftmost dimension of A
+  !
+  ! n : Number of rows of A
+  !
+  ! m : Number of columns of A
+  !
+  END_DOC
+  integer, intent(in)            :: m,n, LDA
+  complex*16, intent(inout) :: A(LDA,n)
+  
+  integer                        :: lwork, info
+  integer, allocatable           :: jpvt(:)
+  complex*16, allocatable  :: tau(:), work(:)
+
+  allocate (jpvt(n), tau(n), work(1))
+  LWORK=-1
+  call  zgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
+  LWORK=2*int(WORK(1))
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+  call  zgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
+  call  zungqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO)
+  deallocate(WORK,jpvt,tau)
+end
+
+subroutine ortho_qr_unblocked_complex(A,LDA,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Orthogonalization using Q.R factorization
+  !
+  ! A : matrix to orthogonalize
+  !
+  ! LDA : leftmost dimension of A
+  !
+  ! n : Number of rows of A
+  !
+  ! m : Number of columns of A
+  !
+  END_DOC
+  integer, intent(in)            :: m,n, LDA
+  double precision, intent(inout) :: A(LDA,n)
+  
+  integer                        :: info
+  integer, allocatable           :: jpvt(:)
+  double precision, allocatable  :: tau(:), work(:)
+
+  print *, irp_here, ': TO DO'
+  stop -1
+
+!  allocate (jpvt(n), tau(n), work(n))
+!  call  dgeqr2( m, n, A, LDA, TAU, WORK, INFO )
+!  call dorg2r(m, n, n, A, LDA, tau, WORK, INFO)
+!  deallocate(WORK,jpvt,tau)
+end
+  
+subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m)
+  implicit none
+  BEGIN_DOC
+  ! Compute C_new=C_old.S^-1/2 orthogonalization.
+  !
+  ! overlap : overlap matrix 
+  !
+  ! LDA : leftmost dimension of overlap array
+  !
+  ! N : Overlap matrix is NxN (array is (LDA,N) )
+  !
+  ! C : Coefficients of the vectors to orthogonalize. On exit,
+  !     orthogonal vectors
+  !
+  ! LDC : leftmost dimension of C
+  !
+  ! M : Coefficients matrix is MxN, ( array is (LDC,N) )
+  !
+  END_DOC
+  
+  integer, intent(in)            :: LDA, ldc, n, m
+  complex*16, intent(in)   :: overlap(lda,n)
+  complex*16, intent(inout) :: C(ldc,n)
+  complex*16, allocatable  :: U(:,:)
+  complex*16, allocatable  :: Vt(:,:)
+  double precision, allocatable  :: D(:)
+  complex*16, allocatable  :: S(:,:)
+  integer                        :: info, i, j, k
+  
+  if (n < 2) then
+    return
+  endif
+
+  allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n))
+
+  call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n)
+
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP SHARED(S,U,D,Vt,n,C,m) &
+  !$OMP PRIVATE(i,j,k)
+
+  !$OMP DO
+  do i=1,n
+    if ( D(i) < 1.d-6 ) then
+      D(i) = 0.d0
+    else
+      D(i) = 1.d0/dsqrt(D(i))
+    endif
+    do j=1,n
+      S(j,i) = (0.d0,0.d0)
+    enddo
+  enddo
+  !$OMP END DO
+
+  do k=1,n
+    if (D(k) /= 0.d0) then
+      !$OMP DO
+      do j=1,n
+        do i=1,n
+          S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j)
+        enddo
+      enddo
+      !$OMP END DO NOWAIT
+    endif
+  enddo
+  
+  !$OMP BARRIER
+  !$OMP DO
+  do j=1,n
+    do i=1,m
+      U(i,j) = C(i,j)
+    enddo
+  enddo
+  !$OMP END DO
+  
+  !$OMP END PARALLEL
+
+  call zgemm('N','N',m,n,n,(1.d0,0.d0),U,size(U,1),S,size(S,1),(0.d0,0.d0),C,size(C,1))
+  
+  deallocate(U,Vt,S,D)
+end
+
+subroutine get_inverse_complex(A,LDA,m,C,LDC)
+  implicit none
+  BEGIN_DOC
+  ! Returns the inverse of the square matrix A
+  END_DOC
+  integer, intent(in)            :: m, LDA, LDC
+  complex*16, intent(in)   :: A(LDA,m)
+  complex*16, intent(out)  :: C(LDC,m)
+
+  integer                        :: info,lwork
+  integer, allocatable           :: ipiv(:)
+  complex*16,allocatable   :: work(:)
+  allocate (ipiv(m), work(m*m))
+  lwork = size(work)
+  C(1:m,1:m) = A(1:m,1:m)
+  call zgetrf(m,m,C,size(C,1),ipiv,info)
+  if (info /= 0) then
+    print *,  info
+    stop 'error in inverse (zgetrf)'
+  endif
+  call zgetri(m,C,size(C,1),ipiv,work,lwork,info)
+  if (info /= 0) then
+    print *,  info
+    stop 'error in inverse (zgetri)'
+  endif
+  deallocate(ipiv,work)
+end
+
+
+subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC)
+  implicit none
+  BEGIN_DOC
+  ! Find C = A^-1
+  END_DOC
+  integer, intent(in)            :: m,n, LDA, LDC
+  complex*16, intent(in)   :: A(LDA,n)
+  complex*16, intent(out)  :: C(LDC,m)
+  
+  double precision, allocatable  :: D(:), rwork(:)
+  complex*16, allocatable  :: U(:,:), Vt(:,:), work(:), A_tmp(:,:)
+  integer                        :: info, lwork
+  integer                        :: i,j,k
+  allocate (D(n),U(m,n),Vt(n,n),work(1),A_tmp(m,n),rwork(5*n))
+  do j=1,n
+    do i=1,m
+      A_tmp(i,j) = A(i,j)
+    enddo
+  enddo
+  lwork = -1
+  call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info)
+  if (info /= 0) then
+    print *,  info, ': SVD failed'
+    stop
+  endif
+  lwork = int(real(work(1)))
+  deallocate(work)
+  allocate(work(lwork))
+  call zgesvd('S','A', m, n, A_tmp, m,D,U,m,Vt,n,work,lwork,rwork,info)
+  if (info /= 0) then
+    print *,  info, ':: SVD failed'
+    stop 1
+  endif
+  
+  do i=1,n
+    if (D(i)/D(1) > 1.d-10) then
+      D(i) = 1.d0/D(i)
+    else
+      D(i) = 0.d0
+    endif
+  enddo
+  
+  C = (0.d0,0.d0)
+  do i=1,m
+    do j=1,n
+      do k=1,n
+        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
+      enddo
+    enddo
+  enddo
+  
+  deallocate(U,D,Vt,work,A_tmp,rwork)
+                                                                                              
+end
+
+subroutine lapack_diagd_diag_in_place_complex(eigvalues,eigvectors,nmax,n)
+  implicit none
+  BEGIN_DOC
+  ! Diagonalize matrix H(complex)
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+!  double precision, intent(out)  :: eigvectors(nmax,n)
+  complex*16, intent(inout)        :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+!  double precision, intent(in)   :: H(nmax,n)
+  complex*16,allocatable         :: work(:)
+  integer         ,allocatable   :: iwork(:)
+!  complex*16,allocatable         :: A(:,:)
+  double precision, allocatable  :: rwork(:)
+  integer                        :: lrwork, lwork, info, i,j,l,k, liwork
+
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  lwork = 2*n*n + 2*n
+  lrwork = 2*n*n + 5*n+ 1
+  liwork = 5*n + 3
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+
+  lwork = -1
+  liwork = -1
+  lrwork = -1
+  ! get optimal work size
+  call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( real(work(1)))
+  liwork = iwork(1)
+  lrwork = int(rwork(1))
+  deallocate (work,iwork,rwork)
+
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+  call ZHEEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  deallocate(work,iwork,rwork)
+
+
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'ZHEEVD Failed; calling ZHEEV'
+     lwork = 2*n - 1
+     lrwork = 3*n - 2
+     allocate(work(lwork),rwork(lrwork))
+     lwork = -1
+     call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info)
+     if (info < 0) then
+       print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value'
+       stop 2
+     endif
+     lwork = int(work(1))
+     deallocate(work)
+     allocate(work(lwork))
+     call ZHEEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,rwork,info)
+     if (info /= 0 ) then
+       write(*,*)'ZHEEV Failed'
+       stop 1
+     endif
+    deallocate(work,rwork)
+  end if
+
+end
+
+subroutine lapack_diagd_diag_complex(eigvalues,eigvectors,H,nmax,n)
+  implicit none                                                                               
+  BEGIN_DOC
+  ! Diagonalize matrix H(complex)
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+!  double precision, intent(out)  :: eigvectors(nmax,n)
+  complex*16, intent(out)        :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+!  double precision, intent(in)   :: H(nmax,n)
+  complex*16, intent(in)         :: H(nmax,n)
+  double precision, allocatable  :: eigenvalues(:)
+  complex*16,allocatable         :: work(:)
+  integer         ,allocatable   :: iwork(:)
+  complex*16,allocatable         :: A(:,:)
+  double precision, allocatable  :: rwork(:)
+  integer                        :: lrwork, lwork, info, i,j,l,k, liwork
+
+  allocate(A(nmax,n),eigenvalues(n))
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  A=H
+  lwork = 2*n*n + 2*n
+  lrwork = 2*n*n + 5*n+ 1
+  liwork = 5*n + 3
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+
+  lwork = -1
+  liwork = -1
+  lrwork = -1
+  ! get optimal work size
+  call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( real(work(1)))
+  liwork = iwork(1)
+  lrwork = int(rwork(1))
+  deallocate (work,iwork,rwork)
+
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+  call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  deallocate(work,iwork,rwork)
+
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'ZHEEVD Failed; calling ZHEEV'
+     lwork = 2*n - 1
+     lrwork = 3*n - 2
+     allocate(work(lwork),rwork(lrwork))
+     lwork = -1
+     call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info)
+     if (info < 0) then
+       print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value'
+       stop 2
+     endif
+     lwork = int(work(1))
+     deallocate(work)
+     allocate(work(lwork))
+     call ZHEEV('V','L',n,A,nmax,eigenvalues,work,lwork,rwork,info)
+     if (info /= 0 ) then
+       write(*,*)'ZHEEV Failed'
+       stop 1
+     endif
+    deallocate(work,rwork)
+  end if
+
+  eigvectors = (0.d0,0.d0)
+  eigvalues = 0.d0
+  do j = 1, n
+    eigvalues(j) = eigenvalues(j)
+    do i = 1, n
+      eigvectors(i,j) = A(i,j)
+    enddo
+  enddo
+  deallocate(A,eigenvalues)
+end
+
+subroutine lapack_diagd_complex(eigvalues,eigvectors,H,nmax,n)
+  implicit none
+  BEGIN_DOC
+  ! Diagonalize matrix H(complex)
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+!  double precision, intent(out)  :: eigvectors(nmax,n)
+  complex*16, intent(out)        :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+!  double precision, intent(in)   :: H(nmax,n)
+  complex*16, intent(in)         :: H(nmax,n)
+  double precision, allocatable  :: eigenvalues(:)
+  complex*16,allocatable         :: work(:)
+  integer         ,allocatable   :: iwork(:)
+  complex*16,allocatable         :: A(:,:)
+  double precision, allocatable  :: rwork(:)
+  integer                        :: lrwork, lwork, info, i,j,l,k, liwork
+
+  allocate(A(nmax,n),eigenvalues(n))
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  A=H
+  lwork = 2*n*n + 2*n
+  lrwork = 2*n*n + 5*n+ 1
+  liwork = 5*n + 3
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+
+  lwork = -1
+  liwork = -1
+  lrwork = -1
+  call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( work( 1 ) )
+  liwork = iwork(1)
+  lrwork = rwork(1)
+  deallocate (work,iwork,rwork)
+
+  allocate (work(lwork),iwork(liwork),rwork(lrwork))
+  call ZHEEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, lrwork, iwork, liwork, info )
+  deallocate(work,iwork,rwork)                                                                
+
+
+  if (info < 0) then
+    print *, irp_here, ': ZHEEVD: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'ZHEEVD Failed'
+     stop 1
+  end if
+
+  eigvectors = (0.d0,0.d0)
+  eigvalues = 0.d0
+  do j = 1, n
+    eigvalues(j) = eigenvalues(j)
+    do i = 1, n
+      eigvectors(i,j) = A(i,j)
+    enddo
+  enddo
+  deallocate(A,eigenvalues)
+end
+
+subroutine lapack_diag_complex(eigvalues,eigvectors,H,nmax,n)
+  implicit none                                                                               
+  BEGIN_DOC
+  ! Diagonalize matrix H (complex)
+  !
+  ! H is untouched between input and ouptut
+  !
+  ! eigevalues(i) = ith lowest eigenvalue of the H matrix
+  !
+  ! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  !
+  END_DOC
+  integer, intent(in)            :: n,nmax
+  complex*16, intent(out)        :: eigvectors(nmax,n)
+  double precision, intent(out)  :: eigvalues(n)
+  complex*16, intent(in)   :: H(nmax,n)
+  double precision,allocatable   :: eigenvalues(:)
+  complex*16,allocatable   :: work(:)
+  complex*16,allocatable   :: A(:,:)
+  double precision,allocatable   :: rwork(:)
+  integer                        :: lwork, info, i,j,l,k,lrwork
+
+  allocate(A(nmax,n),eigenvalues(n))
+! print*,'Diagonalization by jacobi'
+! print*,'n = ',n
+
+  A=H
+  !lwork = 2*n*n + 6*n+ 1
+  lwork = 2*n - 1
+  lrwork = 3*n - 2
+  allocate (work(lwork),rwork(lrwork))
+
+  lwork = -1
+  call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, info )
+  if (info < 0) then
+    print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  endif
+  lwork  = int( work( 1 ) )
+  deallocate (work)
+
+  allocate (work(lwork))
+  call ZHEEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
+    rwork, info )
+  deallocate(work,rwork)
+
+  if (info < 0) then
+    print *, irp_here, ': ZHEEV: the ',-info,'-th argument had an illegal value'
+    stop 2
+  else if( info > 0  ) then
+     write(*,*)'ZHEEV Failed : ', info
+     do i=1,n
+      do j=1,n
+        print *,  H(i,j)
+      enddo
+     enddo 
+     stop 1
+  end if
+  
+  eigvectors = (0.d0,0.d0)
+  eigvalues = 0.d0
+  do j = 1, n
+    eigvalues(j) = eigenvalues(j)
+    do i = 1, n
+      eigvectors(i,j) = A(i,j)
+    enddo
+  enddo
+  deallocate(A,eigenvalues)
+end 
+    
+subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda)
+ implicit none
+ BEGIN_DOC
+! performs u1 += u0 * matrix 
+ END_DOC
+ integer, intent(in)             :: sze,lda
+ complex*16, intent(in)    :: u0(sze)
+ complex*16, intent(inout) :: u1(sze)
+ complex*16, intent(in)    :: matrix(lda,sze)
+ integer :: i,j
+ integer                        :: incx,incy
+ incx = 1
+ incy = 1
+ !call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy)
+ call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy)
+end
+
 subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
   implicit none
   BEGIN_DOC
@@ -356,6 +1040,8 @@ subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC)
 
 end
 
+
+
 subroutine find_rotation(A,LDA,B,m,C,n)
   implicit none
   BEGIN_DOC

src/utils_periodic/import_integrals_ao_periodic.irp.f

@@ -0,0 +1,137 @@
+program print_integrals
+  print *, 'Number of AOs?'
+  read(*,*) ao_num
+  TOUCH ao_num
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j,k,l
+  double precision :: integral
+  double precision, allocatable :: A(:,:), B(:,:)
+  double precision :: tmp_re, tmp_im
+
+  integer             :: n_integrals 
+  integer(key_kind), allocatable   :: buffer_i(:) 
+  real(integral_kind), allocatable :: buffer_values(:)
+
+  call ezfio_set_ao_basis_ao_num(ao_num)
+
+  allocate (A(ao_num,ao_num), B(ao_num,ao_num) )
+  
+  A(1,1) = huge(1.d0)
+  iunit = getunitandopen('E.qp','r')
+  read (iunit,*,end=9) A(1,1)
+  9 continue
+  close(iunit)
+  if (A(1,1) /= huge(1.d0)) then
+    call ezfio_set_nuclei_nuclear_repulsion(A(1,1))
+    call ezfio_set_nuclei_io_nuclear_repulsion("Read")
+  endif
+
+  A = 0.d0
+  B = 0.d0
+  iunit = getunitandopen('T.qp','r')
+  do 
+    read (iunit,*,end=10) i,j, tmp_re, tmp_im
+    A(i,j) = tmp_re
+    B(i,j) = tmp_im
+    if (i.ne.j) then
+      A(j,i) =  tmp_re
+      B(j,i) = -tmp_im
+    endif
+  enddo
+  10 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_ao_integrals_kinetic_imag(B(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic("Read")
+
+  A = 0.d0
+  B = 0.d0
+  iunit = getunitandopen('S.qp','r')
+  do 
+    read (iunit,*,end=11) i,j, tmp_re, tmp_im
+    A(i,j) = tmp_re
+    B(i,j) = tmp_im
+    if (i.ne.j) then
+      A(j,i) =  tmp_re
+      B(j,i) = -tmp_im
+    endif
+  enddo
+  11 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_ao_integrals_overlap_imag(B(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap("Read")
+
+  A = 0.d0
+  B = 0.d0
+  iunit = getunitandopen('P.qp','r')
+  do 
+    read (iunit,*,end=14) i,j, tmp_re, tmp_im
+    A(i,j) = tmp_re
+    B(i,j) = tmp_im
+    if (i.ne.j) then
+      A(j,i) =  tmp_re
+      B(j,i) = -tmp_im
+    endif
+  enddo
+  14 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A(1:ao_num,1:ao_num))
+  call ezfio_set_ao_one_e_ints_ao_integrals_pseudo_imag(B(1:ao_num,1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo("Read")
+
+  A = 0.d0
+  B = 0.d0
+  iunit = getunitandopen('V.qp','r')
+  do 
+    read (iunit,*,end=12) i,j, tmp_re, tmp_im
+    A(i,j) = tmp_re
+    B(i,j) = tmp_im
+    if (i.ne.j) then
+      A(j,i) =  tmp_re
+      B(j,i) = -tmp_im
+    endif
+  enddo
+  12 continue
+  close(iunit)
+  call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_ao_integrals_n_e_imag(B(1:ao_num, 1:ao_num))
+  call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e("Read")
+
+!  allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
+!  iunit = getunitandopen('W.qp','r')
+!  n_integrals=0
+!  buffer_values = 0.d0
+!  do 
+!    read (iunit,*,end=13) i,j,k,l, integral
+!    n_integrals += 1
+!    call two_e_integrals_index(i, j, k, l, buffer_i(n_integrals) )
+!    buffer_values(n_integrals) = integral
+!    if (n_integrals == size(buffer_i)) then
+!      call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
+!      n_integrals = 0
+!    endif
+!  enddo
+!  13 continue
+!  close(iunit)
+!  
+!  if (n_integrals > 0) then
+!    call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
+!  endif
+!
+!  call map_sort(ao_integrals_map)
+!  call map_unique(ao_integrals_map)
+!
+!  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+!  call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+  
+end

src/utils_periodic/import_mo_coef_periodic.irp.f

@@ -0,0 +1,27 @@
+program import_mo_coef_periodic
+
+  PROVIDE ezfio_filename
+  call run
+end
+
+subroutine run
+  use map_module
+  implicit none
+  
+  integer :: iunit
+  integer :: getunitandopen
+
+  integer ::i,j
+  double precision :: int_re, int_im
+
+  
+  iunit = getunitandopen('C.qp','r')
+  do 
+    read (iunit,*,end=10) i,j, mo_coef(i,j), mo_coef_imag(i,j)
+  enddo
+  10 continue
+  close(iunit)
+  mo_label = "None"
+  call save_mos
+
+end