started kpts nos

src/determinants/density_matrix.irp.f

@@ -248,44 +248,49 @@ BEGIN_PROVIDER [ double precision, one_e_spin_density_mo, (mo_num,mo_num) ]
 END_PROVIDER
 
 subroutine set_natural_mos
-   implicit none
-   BEGIN_DOC
-   ! Set natural orbitals, obtained by diagonalization of the one-body density matrix
-   ! in the |MO| basis
-   END_DOC
-   character*(64)                 :: label
-   double precision, allocatable  :: tmp(:,:)
+  implicit none
+  BEGIN_DOC
+  ! Set natural orbitals, obtained by diagonalization of the one-body density matrix
+  ! in the |MO| basis
+  END_DOC
+  character*(64)                 :: label
+  double precision, allocatable  :: tmp(:,:)
 
-   label = "Natural"
-    integer :: i,j,iorb,jorb
-    if (is_complex) then
-    do i = 1, n_virt_orb
-     iorb = list_virt(i)
-     do j = 1, n_core_inact_act_orb
-      jorb = list_core_inact_act(j)
-      if(cdabs(one_e_dm_mo_complex(iorb,jorb)).ne. 0.d0)then
-        print*,'AHAHAH'
-        print*,iorb,jorb,one_e_dm_mo_complex(iorb,jorb)
-        stop
-      endif
-     enddo
+  label = "Natural"
+  integer :: i,j,iorb,jorb,k
+  if (is_complex) then
+    !todo: implement for kpts
+    do k=1,kpt_num
+      do i = 1, n_virt_orb_kpts(k)
+        iorb = list_virt_kpts(i,k)
+        do j = 1, n_core_inact_act_orb_kpts(k)
+          jorb = list_core_inact_act_kpts(j,k)
+          if(cdabs(one_e_dm_mo_kpts(iorb,jorb,k)).ne. 0.d0)then
+            print*,'AHAHAH'
+            print*,iorb,jorb,k,one_e_dm_mo_kpts(iorb,jorb,k)
+            stop
+          endif
+        enddo
+      enddo
     enddo
-   call mo_as_svd_vectors_of_mo_matrix_eig_complex(one_e_dm_mo_complex,size(one_e_dm_mo_complex,1),mo_num,mo_num,mo_occ,label)
-    else
+!    call mo_as_svd_vectors_of_mo_matrix_eig_complex(one_e_dm_mo_complex,size(one_e_dm_mo_complex,1),mo_num,mo_num,mo_occ,label)
+    call mo_as_svd_vectors_of_mo_matrix_eig_kpts(one_e_dm_mo_kpts,size(one_e_dm_mo_kpts,1),mo_num_per_kpt,mo_num_per_kpt,kpt_num,mo_occ_kpts,label)
+    soft_touch mo_occ_kpts
+  else
     do i = 1, n_virt_orb
-     iorb = list_virt(i)
-     do j = 1, n_core_inact_act_orb
-      jorb = list_core_inact_act(j)
-      if(one_e_dm_mo(iorb,jorb).ne. 0.d0)then
-        print*,'AHAHAH'
-        print*,iorb,jorb,one_e_dm_mo(iorb,jorb)
-        stop
-      endif
-     enddo
+      iorb = list_virt(i)
+      do j = 1, n_core_inact_act_orb
+        jorb = list_core_inact_act(j)
+        if(one_e_dm_mo(iorb,jorb).ne. 0.d0)then
+          print*,'AHAHAH'
+          print*,iorb,jorb,one_e_dm_mo(iorb,jorb)
+          stop
+        endif
+      enddo
     enddo
-   call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label)
-    endif
-   soft_touch mo_occ
+    call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label)
+    soft_touch mo_occ
+  endif
 
 end
 subroutine save_natural_mos

src/mo_basis/utils.irp.f

@@ -12,7 +12,7 @@ subroutine save_mos
   call ezfio_set_mo_basis_mo_label(mo_label)
   call ezfio_set_mo_basis_ao_md5(ao_md5)
   if (is_complex) then
-    allocate ( buffer_c(ao_num,mo_num))
+    !allocate ( buffer_c(ao_num,mo_num))
     allocate ( buffer_k(ao_num_per_kpt,mo_num_per_kpt,kpt_num))
     buffer_k = (0.d0,0.d0)
     do k=1,kpt_num
@@ -53,6 +53,8 @@ subroutine save_mos_no_occ
  !call ezfio_set_mo_basis_mo_label(mo_label)
  !call ezfio_set_mo_basis_ao_md5(ao_md5)
   if (is_complex) then
+    print*,irp_here, ' not implemented for kpts'
+    stop -1
     allocate ( buffer_c(ao_num,mo_num))
     buffer_c = (0.d0,0.d0)
     do j = 1, mo_num
@@ -88,6 +90,8 @@ subroutine save_mos_truncated(n)
   call ezfio_set_mo_basis_mo_label(mo_label)
   call ezfio_set_mo_basis_ao_md5(ao_md5)
   if (is_complex) then
+    print*,irp_here, ' not implemented for kpts'
+    stop -1
     allocate ( buffer_c(ao_num,mo_num))
     buffer_c = (0.d0,0.d0)
     do j = 1, n

src/mo_basis/utils_cplx.irp.f

@@ -384,63 +384,74 @@ subroutine mo_as_svd_vectors_of_mo_matrix_kpts(matrix,lda,m,n,label)
 end
 
 
-subroutine mo_as_svd_vectors_of_mo_matrix_eig_kpts(matrix,lda,m,n,eig,label)
+subroutine mo_as_svd_vectors_of_mo_matrix_eig_kpts(matrix,lda,m,n,nk,eig,label)
   !TODO: implement
-  print *, irp_here, ' not implemented for kpts'
-  stop 1
+  !print *, irp_here, ' not implemented for kpts'
+  !stop 1
   implicit none
-  integer,intent(in)             :: lda,m,n
+  integer,intent(in)             :: lda,m,n,nk
   character*(64), intent(in)     :: label
-  complex*16, intent(in)   :: matrix(lda,n)
-  double precision, intent(out)  :: eig(m)
+  complex*16, intent(in)   :: matrix(lda,n,nk)
+  double precision, intent(out)  :: eig(m,nk)
 
-  integer :: i,j
+  integer :: i,j,k
   double precision :: accu
   double precision, allocatable  :: D(:)
   complex*16, allocatable  :: mo_coef_new(:,:), U(:,:), A(:,:), Vt(:,:), work(:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
 
   call write_time(6)
-  if (m /= mo_num) then
-    print *, irp_here, ': Error : m/= mo_num'
+  if (m /= mo_num_per_kpt) then
+    print *, irp_here, ': Error : m/= mo_num_per_kpt'
     stop 1
   endif
 
-  allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num,m),D(m),Vt(lda,n))
+  
+  allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num_per_kpt,m),D(m),Vt(lda,n))
+
+  do k=1,nk
+    do j=1,n
+      do i=1,m
+        A(i,j) = matrix(i,j,k)
+      enddo
+    enddo
+    mo_coef_new = mo_coef_kpts(1,1,k)
+
+    call svd_complex(A,lda,U,lda,D,Vt,lda,m,n)
+
+
+
+    call zgemm('N','N',ao_num_per_kpt,m,m,  &
+               (1.d0,0.d0),mo_coef_new,size(mo_coef_new,1),U,size(U,1),&
+               (0.d0,0.d0),mo_coef_kpts(1,1,k),size(mo_coef_kpts,1))
 
-  do j=1,n
     do i=1,m
-      A(i,j) = matrix(i,j)
+      eig(i,k) = D(i)
     enddo
   enddo
-  mo_coef_new = mo_coef_complex
 
-  call svd_complex(A,lda,U,lda,D,Vt,lda,m,n)
+  deallocate(A,mo_coef_new,U,Vt,D)
 
   write (6,'(A)') 'MOs are now **'//trim(label)//'**'
   write (6,'(A)')  ''
-  write (6,'(A)') 'Eigenvalues'
+  write (6,'(A)') 'Eigenvalues '
   write (6,'(A)')  '-----------'
   write (6,'(A)') ''
   write (6,'(A)')  '======== ================ ================'
   write (6,'(A)')  '   MO       Eigenvalue       Cumulative   '
   write (6,'(A)')  '======== ================ ================'
-
-  accu = 0.d0
-  do i=1,m
-    accu = accu + D(i)
-    write (6,'(I8,1X,F16.10,1X,F16.10)')  i,D(i), accu
+ 
+  do k=1,nk
+    accu = 0.d0
+    do i=1,m
+      accu = accu + eig(i,k)
+      write (6,'(I8,1X,F16.10,1X,F16.10)')  i,eig(i,k), accu
+    enddo
+  write (6,'(A)')  '-------- ---------------- ----------------'
   enddo
   write (6,'(A)')  '======== ================ ================'
   write (6,'(A)')  ''
 
-  call zgemm('N','N',ao_num,m,m,(1.d0,0.d0),mo_coef_new,size(mo_coef_new,1),U,size(U,1),(0.d0,0.d0),mo_coef_complex,size(mo_coef_complex,1))
-
-  do i=1,m
-    eig(i) = D(i)
-  enddo
-
-  deallocate(A,mo_coef_new,U,Vt,D)
   call write_time(6)
 
   mo_label = label

src/scf_utils/roothaan_hall_scf_cplx.irp.f

@@ -344,7 +344,7 @@ END_DOC
   logical, external              :: qp_stop
   complex*16, allocatable :: mo_coef_save(:,:,:)
 
-  PROVIDE ao_md5 mo_occ level_shift
+  PROVIDE ao_md5 mo_occ_kpts level_shift
 
   allocate(mo_coef_save(ao_num_per_kpt,mo_num_per_kpt,kpt_num),                          &
       Fock_matrix_DIIS (ao_num_per_kpt,ao_num_per_kpt,max_dim_DIIS,kpt_num),                 &