Merge remote-tracking branch 'origin/dev' into cleaning_dft

REPLACE

@@ -1,5 +1,4 @@
 # This file contains all the renamings that occured between qp1 and qp2.
-#
 qp_name aa_operator_bielec -r aa_operator_two_e
 qp_name ac_operator_bielec -r ac_operator_two_e
 qp_name ao_bi_elec_integral_alpha -r ao_two_e_integral_alpha
@@ -222,6 +221,11 @@ qp_name potential_sr_xc_beta_ao_lda  --rename=potential_xc_beta_ao_sr_lda
 qp_name potential_sr_xc_beta_ao_lda --rename=potential_xc_beta_ao_sr_lda
 qp_name potential_sr_xc_beta_ao_pbe  --rename=potential_xc_beta_ao_sr_pbe
 qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe
+qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion
+qp_name nucl_elec_ref_bitmask_energy -r ref_bitmask_n_e_energy
+qp_name ref_bitmask_e_n_energy -r ref_bitmask_n_e_energy
+qp_name read_ao_integrals_e_n -r read_ao_integrals_n_e
+qp_name write_ao_integrals_e_n -r write_ao_integrals_n_e
 qp_name psi_energy_bielec -r psi_energy_two_e
 qp_name read_ao_integrals --rename="read_ao_two_e_integrals"
 qp_name read_ao_integrals --rename=read_ao_two_e_integrals
@@ -240,3 +244,12 @@ qp_name write_ao_integrals --rename=write_ao_two_e_integrals
 qp_name write_mo_integrals_erf -r write_mo_two_e_integrals_erf
 qp_name write_mo_integrals --rename="write_mo_two_e_integrals"
 qp_name write_mo_integrals --rename=write_mo_two_e_integrals
+qp_name io_mo_integrals_e_n -r io_mo_integrals_n_e
+qp_name write_mo_integrals_e_n -r write_mo_integrals_n_e
+qp_name read_mo_integrals_e_n -r read_mo_integrals_n_e
+qp_name mo_integrals_e_n -r mo_integrals_n_e
+qp_name ezfio_get_mo_one_e_ints_mo_integrals_e_n -r ezfio_get_mo_one_e_ints_mo_integrals_n_e
+qp_name ezfio_set_mo_one_e_ints_mo_integrals_e_n -r ezfio_set_mo_one_e_ints_mo_integrals_n_e
+qp_name ezfio_has_mo_one_e_ints_mo_integrals_e_n -r ezfio_has_mo_one_e_ints_mo_integrals_n_e
+qp_name ezfio_has_mo_one_e_ints_io_mo_integrals_e_n -r ezfio_has_mo_one_e_ints_io_mo_integrals_n_e
+qp_name ezfio_get_mo_one_e_ints_io_mo_integrals_e_n -r ezfio_get_mo_one_e_ints_io_mo_integrals_n_e

ocaml/Input_mo_basis.ml

@@ -10,6 +10,7 @@ module Mo_basis : sig
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array option;
         ao_md5          : MD5.t;
       } [@@deriving sexp]
   val read : unit -> t option
@@ -24,6 +25,7 @@ end = struct
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array option;
         ao_md5          : MD5.t;
       } [@@deriving sexp]
   let get_default = Qpackage.get_ezfio_default "mo_basis"
@@ -37,11 +39,18 @@ end = struct
 
 
   let reorder b ordering =
-    { b with mo_coef =
+    { b with
-      Array.map (fun mo ->
+      mo_coef = Array.map (fun mo ->
           Array.init (Array.length mo)
             (fun i -> mo.(ordering.(i)))
-      ) b.mo_coef 
+        ) b.mo_coef  ;
+      mo_coef_imag =
+        match b.mo_coef_imag with
+        | None -> None
+        | Some x -> Some ( Array.map (fun mo ->
+              Array.init (Array.length mo)
+                (fun i -> mo.(ordering.(i)))
+            ) x )
     }
 
   let read_ao_md5 () =
@@ -121,6 +130,20 @@ end = struct
         Array.sub a (j*ao_num) (ao_num) 
       )
 
+  let read_mo_coef_imag () =
+    if Ezfio.has_mo_basis_mo_coef_imag () then
+      let a =
+          Ezfio.get_mo_basis_mo_coef_imag ()
+                  |> Ezfio.flattened_ezfio
+                  |> Array.map MO_coef.of_float 
+      in
+      let mo_num = read_mo_num () |> MO_number.to_int in
+      let ao_num = (Array.length a)/mo_num in
+      Some (Array.init mo_num (fun j ->
+          Array.sub a (j*ao_num) (ao_num) 
+        ) )
+    else None
+
 
   let read () =
     if (Ezfio.has_mo_basis_mo_num ()) then
@@ -130,6 +153,7 @@ end = struct
           mo_class        = read_mo_class ();
           mo_occ          = read_mo_occ ();
           mo_coef         = read_mo_coef ();
+          mo_coef_imag    = read_mo_coef_imag ();
           ao_md5          = read_ao_md5 ();
         }
     else
@@ -137,6 +161,7 @@ end = struct
 
 
   let mo_coef_to_string mo_coef =
+  (*TODO : add imaginary part here *)
     let ao_num = Array.length mo_coef.(0)
     and mo_num = Array.length mo_coef in
     let rec print_five imin imax =
@@ -222,6 +247,7 @@ MO coefficients ::
 
 
   let to_string b =
+  (*TODO : add imaginary part here *)
     Printf.sprintf "
 mo_label        = \"%s\"
 mo_num          = %s
@@ -244,12 +270,12 @@ mo_coef         = %s
   let write_mo_num n =
     MO_number.to_int n
     |> Ezfio.set_mo_basis_mo_num
-  ;;
+
 
   let write_mo_label a =
     MO_label.to_string a
     |> Ezfio.set_mo_basis_mo_label
-  ;;
+
 
   let write_mo_class a =
     let mo_num = Array.length a in
@@ -257,7 +283,7 @@ mo_coef         = %s
     |> Array.to_list
     in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| mo_num |] ~data
     |> Ezfio.set_mo_basis_mo_class
-  ;;
+
 
   let write_mo_occ a =
     let mo_num = Array.length a in
@@ -265,12 +291,12 @@ mo_coef         = %s
     |> Array.to_list
     in Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| mo_num |] ~data
     |> Ezfio.set_mo_basis_mo_occ
-  ;;
+
 
   let write_md5 a =
     MD5.to_string a
     |> Ezfio.set_mo_basis_ao_md5
-  ;;
+
 
   let write_mo_coef a =
     let mo_num = Array.length a in
@@ -282,7 +308,24 @@ mo_coef         = %s
     |> List.concat
     in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data
     |> Ezfio.set_mo_basis_mo_coef
-  ;;
+
+
+  let write_mo_coef_imag a =
+    match a with
+    | None -> ()
+    | Some a -> 
+      begin
+        let mo_num = Array.length a in
+        let ao_num = Array.length a.(0) in
+        let data =
+          Array.map (fun mo -> Array.map MO_coef.to_float mo
+          |> Array.to_list) a
+        |> Array.to_list
+        |> List.concat
+        in Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; mo_num |] ~data
+        |> Ezfio.set_mo_basis_mo_coef_imag
+      end
+
 
   let write 
       { mo_num          : MO_number.t ;
@@ -290,6 +333,7 @@ mo_coef         = %s
         mo_class        : MO_class.t array;
         mo_occ          : MO_occ.t array;
         mo_coef         : (MO_coef.t array) array;
+        mo_coef_imag    : (MO_coef.t array) array option;
         ao_md5          : MD5.t;
       } =
       write_mo_num mo_num;
@@ -297,8 +341,9 @@ mo_coef         = %s
       write_mo_class mo_class;
       write_mo_occ mo_occ;
       write_mo_coef mo_coef;
+      write_mo_coef_imag mo_coef_imag;
       write_md5 ao_md5
-  ;;
+
 
 end
 

src/ao_one_e_ints/EZFIO.cfg

@@ -1,10 +1,16 @@
-[ao_integrals_e_n]
+[ao_integrals_n_e]
 type: double precision
 doc: Nucleus-electron integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
-[io_ao_integrals_e_n]
+[ao_integrals_n_e_imag]
+type: double precision
+doc: Imaginary part of the nucleus-electron integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
+[io_ao_integrals_n_e]
 type: Disk_access
 doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
@@ -17,6 +23,12 @@ doc: Kinetic energy integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_kinetic_imag]
+type: double precision
+doc: Imaginary part of the kinetic energy integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_kinetic]
 type: Disk_access
 doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ]
@@ -30,6 +42,12 @@ doc: Pseudopotential integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_pseudo_imag]
+type: double precision
+doc: Imaginary part of the pseudopotential integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_pseudo]
 type: Disk_access
 doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ]
@@ -43,6 +61,12 @@ doc: Overlap integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num)
 interface: ezfio
 
+[ao_integrals_overlap_imag]
+type: double precision
+doc: Imaginary part of the overlap integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num)
+interface: ezfio
+
 [io_ao_integrals_overlap]
 type: Disk_access
 doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ]
@@ -56,6 +80,12 @@ doc: Combined integrals in |AO| basis set
 size: (ao_basis.ao_num,ao_basis.ao_num) 
 interface: ezfio
 
+[ao_one_e_integrals_imag]
+type: double precision
+doc: Imaginary part of the combined integrals in |AO| basis set
+size: (ao_basis.ao_num,ao_basis.ao_num) 
+interface: ezfio
+
 [io_ao_one_e_integrals]
 type: Disk_access
 doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ]

src/ao_one_e_ints/ao_one_e_ints.irp.f

@@ -27,3 +27,24 @@
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)]
+  implicit none
+  integer :: i,j,n,l
+  BEGIN_DOC
+ ! One-electron Hamiltonian in the |AO| basis.
+  END_DOC
+
+  IF (read_ao_one_e_integrals) THEN
+     call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag)
+  ELSE
+     print *,  irp_here, ': Not yet implemented'
+     stop -1
+  ENDIF
+
+  IF (write_ao_one_e_integrals) THEN
+       call ezfio_set_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals_imag)
+       print *,  'AO one-e integrals written to disk'
+  ENDIF
+
+END_PROVIDER
+

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,15 +109,22 @@ 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 (is_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)                               &
         !$OMP DEFAULT(NONE)                                          &
-  !$OMP PRIVATE(A_center,B_center,power_A,power_B,&
+        !$OMP PRIVATE(A_center,B_center,power_A,power_B,             &
         !$OMP  overlap_x,overlap_y, overlap_z, overlap,              &
         !$OMP  alpha, beta,i,j,dx)                                   &
         !$OMP SHARED(nucl_coord,ao_power,ao_prim_num,                &
-  !$OMP  ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
+        !$OMP  ao_overlap_abs,ao_num,ao_coef_normalized_ordered_transp,ao_nucl,&
         !$OMP  ao_expo_ordered_transp,dim1,lower_exp_val)
     do j=1,ao_num
       A_center(1) = nucl_coord( ao_nucl(j), 1 )
@@ -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_one_e_ints/kin_ao_ints.irp.f

@@ -149,3 +149,25 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integrals, (ao_num,ao_num)]
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [double precision, ao_kinetic_integrals_imag, (ao_num,ao_num)]
+  implicit none
+  BEGIN_DOC
+  ! Kinetic energy integrals in the |AO| basis.
+  !
+  ! $\langle \chi_i |\hat{T}| \chi_j \rangle$
+  !
+  END_DOC
+  integer                        :: i,j,k,l
+
+  if (read_ao_integrals_kinetic) then
+    call ezfio_get_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
+    print *,  'AO kinetic integrals read from disk'
+  else
+    print *,  irp_here, ': Not yet implemented'
+  endif
+  if (write_ao_integrals_kinetic) then
+    call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(ao_kinetic_integrals_imag)
+    print *,  'AO kinetic integrals written to disk'
+  endif
+END_PROVIDER
+

src/ao_one_e_ints/pot_ao_ints.irp.f

@@ -12,8 +12,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
   integer                        :: i,j,k,l,n_pt_in,m
   double precision               :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
 
-  if (read_ao_integrals_e_n) then
+  if (read_ao_integrals_n_e) then
-    call ezfio_get_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e)
+    call ezfio_get_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
     print *,  'AO N-e integrals read from disk'
   else
 
@@ -76,13 +76,36 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
     !$OMP END DO
     !$OMP END PARALLEL
   endif
-  if (write_ao_integrals_e_n) then
+  if (write_ao_integrals_n_e) then
-    call ezfio_set_ao_one_e_ints_ao_integrals_e_n(ao_integrals_n_e)
+    call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
     print *,  'AO N-e integrals written to disk'
   endif
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, ao_integrals_n_e_imag, (ao_num,ao_num)]
+  BEGIN_DOC
+  !  Nucleus-electron interaction, in the |AO| basis set.
+  !
+  !  :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
+  END_DOC
+  implicit none
+  double precision               :: alpha, beta, gama, delta
+  integer                        :: num_A,num_B
+  double precision               :: A_center(3),B_center(3),C_center(3)
+  integer                        :: power_A(3),power_B(3)
+  integer                        :: i,j,k,l,n_pt_in,m
+  double precision               :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
+
+  if (read_ao_integrals_n_e) then
+    call ezfio_get_ao_one_e_ints_ao_integrals_n_e_imag(ao_integrals_n_e_imag)
+    print *,  'AO N-e integrals read from disk'
+  else
+   print *,  irp_here, ': Not yet implemented'
+  endif
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)]
   BEGIN_DOC
 ! Nucleus-electron interaction in the |AO| basis set, per atom A.
@@ -166,7 +189,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
   double precision               :: P_center(3)
   double precision               :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor
   double precision               :: I_n_special_exact,integrate_bourrin,I_n_bibi
-  double precision               :: V_e_n,const_factor,dist_integral,tmp
+  double precision               :: V_n_e,const_factor,dist_integral,tmp
   double precision               :: accu,epsilo,rint
   integer                        :: n_pt_out,lmax
   include 'utils/constants.include.F'
@@ -178,7 +201,7 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b
         (A_center(3)/=C_center(3))) then
     continue
   else
-    NAI_pol_mult = V_e_n(power_A(1),power_A(2),power_A(3),           &
+    NAI_pol_mult = V_n_e(power_A(1),power_A(2),power_A(3),           &
         power_B(1),power_B(2),power_B(3),alpha,beta)
     return
   endif
@@ -476,7 +499,7 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
   endif
 end
 
-double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
+double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
   implicit none
   BEGIN_DOC
 ! Primitve nuclear attraction between the two primitves centered on the same atom.
@@ -489,9 +512,9 @@ double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)
   double precision               :: alpha,beta
   double precision               :: V_r, V_phi, V_theta
   if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then
-    V_e_n = 0.d0
+    V_n_e = 0.d0
   else
-    V_e_n =   V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta)              &
+    V_n_e =   V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta)              &
         * V_phi(a_x+b_x,a_y+b_y)                                     &
         * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1)
   endif

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 = dcmplx(tmp_re,tmp_im)
+           else
+             integral = dcmplx(tmp_re,-tmp_im)
+           endif
+         else
+           tmp_im = 0.d0
+           integral = dcmplx(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 = dcmplx(tmp_re,tmp_im)
+           else
+             tmp = dcmplx(tmp_re,-tmp_im)
+           endif
+         else
+           tmp_im = 0.d0
+           tmp = dcmplx(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/bitmask/track_orb.irp.f

@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, mo_coef_begin_iteration, (ao_num,mo_num) ]
    BEGIN_DOC
    ! Void provider to store the coefficients of the |MO| basis at the beginning of the SCF iteration
    !
-   ! Usefull to track some orbitals
+   ! Useful to track some orbitals
    END_DOC
 END_PROVIDER
 

src/mo_basis/EZFIO.cfg

@@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO|
 interface: ezfio
 size: (ao_basis.ao_num,mo_basis.mo_num)
 
+[mo_coef_imag]
+type: double precision
+doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)
+
 [mo_label]
 type: character*(64)
 doc: Label characterizing the MOS (Local, Canonical, Natural, *etc*)

src/mo_basis/NEED

@@ -1,3 +1,3 @@
 ao_basis
-ao_one_e_ints
 electrons
+ao_one_e_ints

src/mo_basis/mos.irp.f

@@ -93,6 +93,59 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ]
   endif
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, mo_coef_imag, (ao_num,mo_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Molecular orbital coefficients on |AO| basis set
+  !
+  ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO|
+  !
+  ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc)
+  END_DOC
+  integer                        :: i, j
+  double precision, allocatable  :: buffer(:,:)
+  logical                        :: exists
+  PROVIDE ezfio_filename
+
+
+  if (mpi_master) then
+    ! Coefs
+    call ezfio_has_mo_basis_mo_coef_imag(exists)
+  endif
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_coef_imag with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag)
+      write(*,*) 'Read  mo_coef_imag'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_coef_imag with MPI'
+      endif
+    IRP_ENDIF
+  else
+    ! Orthonormalized AO basis
+    do i=1,mo_num
+      do j=1,ao_num
+        mo_coef_imag(j,i) = 0.d0
+      enddo
+    enddo
+  endif
+END_PROVIDER
+
 BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ]
  implicit none
  BEGIN_DOC

src/mo_one_e_ints/EZFIO.cfg

@@ -1,10 +1,10 @@
-[mo_integrals_e_n]
+[mo_integrals_n_e]
 type: double precision
 doc: Nucleus-electron integrals in |MO| basis set
 size: (mo_basis.mo_num,mo_basis.mo_num)
 interface: ezfio
 
-[io_mo_integrals_e_n]
+[io_mo_integrals_n_e]
 type: Disk_access
 doc: Read/Write |MO| electron-nucleus attraction integrals from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml

src/mo_one_e_ints/pot_mo_ints.irp.f

@@ -4,8 +4,8 @@ BEGIN_PROVIDER [double precision, mo_integrals_n_e, (mo_num,mo_num)]
 ! Nucleus-electron interaction on the |MO| basis
  END_DOC
 
-  if (read_mo_integrals_e_n) then
+  if (read_mo_integrals_n_e) then
-     call ezfio_get_mo_one_e_ints_mo_integrals_e_n(mo_integrals_n_e)
+     call ezfio_get_mo_one_e_ints_mo_integrals_n_e(mo_integrals_n_e)
     print *,  'MO N-e integrals read from disk'
   else
     call ao_to_mo(                                                   &
@@ -15,8 +15,8 @@ BEGIN_PROVIDER [double precision, mo_integrals_n_e, (mo_num,mo_num)]
         size(mo_integrals_n_e,1)                                &
         )
   endif
-  if (write_mo_integrals_e_n) then
+  if (write_mo_integrals_n_e) then
-     call ezfio_set_mo_one_e_ints_mo_integrals_e_n(mo_integrals_n_e)
+     call ezfio_set_mo_one_e_ints_mo_integrals_n_e(mo_integrals_n_e)
     print *,  'MO N-e integrals written to disk'
   endif
 

src/nuclei/EZFIO.cfg

@@ -31,3 +31,9 @@ default: None
 doc: Nuclear repulsion (Computed automaticaly or Read in the |EZFIO|)
 type:double precision
 interface: ezfio
+
+[is_periodic]
+type: logical
+doc: If true, the calculation uses periodic boundary conditions
+interface: ezfio, provider, ocaml
+default: false

src/tools/save_natorb.irp.f

@@ -19,7 +19,7 @@ program save_natorb
   call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None')
   call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None')
   call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None')
-  call ezfio_set_mo_one_e_ints_io_mo_integrals_e_n('None')
+  call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None')
   call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None')
 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  = max(int( work( 1 ) ),lwork)
+  liwork = iwork(1)
+  lrwork = max(int(rwork(1),4),lrwork)
+  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/NEED

@@ -0,0 +1,4 @@
+ao_two_e_ints
+ao_one_e_ints
+mo_two_e_ints
+

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