From b7493137624432e6ea1b1770c8ee7275699af3c0 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Fri, 3 Apr 2020 14:49:11 -0500 Subject: [PATCH] started kpts nos --- src/determinants/density_matrix.irp.f | 73 ++++++++++++---------- src/mo_basis/utils.irp.f | 6 +- src/mo_basis/utils_cplx.irp.f | 65 +++++++++++-------- src/scf_utils/roothaan_hall_scf_cplx.irp.f | 2 +- 4 files changed, 83 insertions(+), 63 deletions(-) diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index e645c390..7d9bf873 100644 --- a/src/determinants/density_matrix.irp.f +++ b/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 diff --git a/src/mo_basis/utils.irp.f b/src/mo_basis/utils.irp.f index 5d94e853..b4ccddb6 100644 --- a/src/mo_basis/utils.irp.f +++ b/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 diff --git a/src/mo_basis/utils_cplx.irp.f b/src/mo_basis/utils_cplx.irp.f index 13327f57..5ca82213 100644 --- a/src/mo_basis/utils_cplx.irp.f +++ b/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 diff --git a/src/scf_utils/roothaan_hall_scf_cplx.irp.f b/src/scf_utils/roothaan_hall_scf_cplx.irp.f index e074f884..87c33cb5 100644 --- a/src/scf_utils/roothaan_hall_scf_cplx.irp.f +++ b/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), &