From 84531d8021763996777b709f1ace6fa2efb1938c Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 17 Mar 2020 17:57:56 -0500 Subject: [PATCH] working on kpts --- src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f | 14 +- src/bitmask/core_inact_act_virt.irp.f | 441 +++++++++++++++++++++ src/hartree_fock/scf.irp.f | 21 +- src/mo_basis/mos_cplx.irp.f | 274 +++++++++++++ src/mo_basis/utils_cplx.irp.f | 267 +++++++++++++ src/mo_guess/mo_ortho_lowdin_cplx.irp.f | 61 +++ src/scf_utils/diagonalize_fock_cplx.irp.f | 62 +++ src/scf_utils/fock_matrix_cplx.irp.f | 374 +++++++++++++++++ src/scf_utils/huckel_cplx.irp.f | 49 +++ src/utils_complex/qp2-pbc-diff.txt | 4 +- 10 files changed, 1551 insertions(+), 16 deletions(-) diff --git a/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f b/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f index 1cc4ba9d..01a02f02 100644 --- a/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f +++ b/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f @@ -12,7 +12,7 @@ integer :: ibegin,j,k integer :: prev prev = 0 - ao_cart_to_sphe_coefi_kpts(:,:) = (0.d0,0.d0) + ao_cart_to_sphe_coef_kpts(:,:) = (0.d0,0.d0) ! Assume order provided by ao_power_index i = 1 ao_cart_to_sphe_num_per_kpt = 0 @@ -79,13 +79,13 @@ BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_overlap_kpts, (ao_cart_to_sphe_num_ call zgemm('T','N',ao_cart_to_sphe_num_per_kpt,ao_num_per_kpt,ao_num_per_kpt, (1.d0,0.d0), & ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1), & - ao_overlap_kpts(1,1,k),size(ao_overlap_kpts,1), (0.d0,0.d0), & + ao_overlap_kpts(:,:,k),size(ao_overlap_kpts,1), (0.d0,0.d0), & S, size(S,1)) call zgemm('N','N',ao_cart_to_sphe_num_per_kpt,ao_cart_to_sphe_num_per_kpt,ao_num_per_kpt, (1.d0,0.d0), & S, size(S,1), & ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1), (0.d0,0.d0), & - ao_cart_to_sphe_overlap_kpts(1,1,k),size(ao_cart_to_sphe_overlap_kpts,1)) + ao_cart_to_sphe_overlap_kpts(:,:,k),size(ao_cart_to_sphe_overlap_kpts,1)) enddo deallocate(S) @@ -106,7 +106,7 @@ BEGIN_PROVIDER [ complex*16, ao_ortho_cano_coef_inv_kpts, (ao_num_per_kpt,ao_num enddo END_PROVIDER - BEGIN_PROVIDER [ complex*16, ao_ortho_canonical_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt)] + BEGIN_PROVIDER [ complex*16, ao_ortho_canonical_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num)] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num_per_kpt, (kpt_num) ] &BEGIN_PROVIDER [ integer, ao_ortho_canonical_num_per_kpt_max ] implicit none @@ -155,14 +155,14 @@ END_PROVIDER ao_cart_to_sphe_num_per_kpt, S, size(S,1), ao_ortho_canonical_num_per_kpt(k)) call zgemm('N','N', ao_num_per_kpt, ao_ortho_canonical_num_per_kpt(k), ao_cart_to_sphe_num_per_kpt, (1.d0,0.d0), & - ao_cart_to_sphe_coef_kpts(:,:,k), size(ao_cart_to_sphe_coef_kpts,1), & + ao_cart_to_sphe_coef_kpts, size(ao_cart_to_sphe_coef_kpts,1), & S, size(S,1), & (0.d0,0.d0), ao_ortho_canonical_coef_kpts(:,:,k), size(ao_ortho_canonical_coef_kpts,1)) enddo deallocate(S) endif - ao_ortho_canonical_num_per_kpt_max = max(ao_ortho_canonical_num_per_kpt) + ao_ortho_canonical_num_per_kpt_max = maxval(ao_ortho_canonical_num_per_kpt) END_PROVIDER BEGIN_PROVIDER [complex*16, ao_ortho_canonical_overlap_kpts, (ao_ortho_canonical_num_per_kpt_max,ao_ortho_canonical_num_per_kpt_max,kpt_num)] @@ -176,7 +176,7 @@ BEGIN_PROVIDER [complex*16, ao_ortho_canonical_overlap_kpts, (ao_ortho_canonical do k=1,kpt_num do j=1, ao_ortho_canonical_num_per_kpt_max do i=1, ao_ortho_canonical_num_per_kpt_max - ao_ortho_canonical_overlap_complex(i,j,k) = (0.d0,0.d0) + ao_ortho_canonical_overlap_kpts(i,j,k) = (0.d0,0.d0) enddo enddo enddo diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index d30e989f..c0484057 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -413,3 +413,444 @@ END_PROVIDER print *, list_inact_act(1:n_inact_act_orb) END_PROVIDER +!============================================! +! ! +! kpts ! +! ! +!============================================! +BEGIN_PROVIDER [ integer, n_core_orb_kpts, (kpt_num)] + implicit none + BEGIN_DOC + ! Number of core MOs + END_DOC + integer :: i,k,kshift + + do k=1,kpt_num + n_core_orb_kpts(k) = 0 + kshift = (1-k)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+kshift) == 'Core')then + n_core_orb_kpts(k) += 1 + endif + enddo + enddo + +! call write_int(6,n_core_orb, 'Number of core MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_inact_orb_kpts, (kpt_num)] + implicit none + BEGIN_DOC + ! Number of inactive MOs + END_DOC + integer :: i,k,kshift + + do k=1,kpt_num + n_inact_orb_kpts(k) = 0 + kshift = (1-k)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+kshift) == 'Inactive')then + n_inact_orb_kpts(k) += 1 + endif + enddo + enddo + +! call write_int(6,n_inact_orb, 'Number of inactive MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_act_orb_kpts, (kpt_num)] + implicit none + BEGIN_DOC + ! Number of active MOs + END_DOC + integer :: i,k,kshift + + do k=1,kpt_num + n_act_orb_kpts(k) = 0 + kshift = (1-k)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+kshift) == 'Active')then + n_act_orb_kpts(k) += 1 + endif + enddo + enddo + +! call write_int(6,n_act_orb, 'Number of active MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_virt_orb_kpts, (kpt_num)] + implicit none + BEGIN_DOC + ! Number of virtual MOs + END_DOC + integer :: i,k,kshift + + do k=1,kpt_num + n_virt_orb_kpts(k) = 0 + kshift = (1-k)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+kshift) == 'Virtual')then + n_virt_orb_kpts(k) += 1 + endif + enddo + enddo + +! call write_int(6,n_virt_orb, 'Number of virtual MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_del_orb_kpts, (kpt_num)] + implicit none + BEGIN_DOC + ! Number of deleted MOs + END_DOC + integer :: i,k,kshift + + do k=1,kpt_num + n_del_orb_kpts(k) = 0 + kshift = (1-k)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+kshift) == 'Deleted')then + n_del_orb_kpts(k) += 1 + endif + enddo + enddo + +! call write_int(6,n_del_orb, 'Number of deleted MOs') + +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_core_inact_orb_kpts, (kpt_num) ] + !todo: finish implementation for kpts (will need kpt_mask) + implicit none + BEGIN_DOC + ! n_core + n_inact + END_DOC + integer :: i,k + do k=1,kpt_num + n_core_inact_orb_kpts(k) = 0 + do i = 1, N_int + n_core_inact_orb_kpts(k) += popcnt(iand(kpt_mask(i,k),reunion_of_core_inact_bitmask(i,1))) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer, n_inact_act_orb_kpts, (kpt_num) ] + implicit none + BEGIN_DOC + ! n_inact + n_act + END_DOC + integer :: k + do k=1,kpt_num + n_inact_act_orb_kpts(k) = (n_inact_orb_kpts(k)+n_act_orb_kpts(k)) + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_core_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_core. + ! it is at least 1 + END_DOC + dim_list_core_orb_kpts = max(maxval(n_core_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_inact_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_inact. + ! it is at least 1 + END_DOC + dim_list_inact_orb_kpts = max(maxval(n_inact_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_core_inact_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_core. + ! it is at least 1 + END_DOC + dim_list_core_inact_orb_kpts = max(maxval(n_core_inact_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_act_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_act. + ! it is at least 1 + END_DOC + dim_list_act_orb_kpts = max(maxval(n_act_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_virt_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_virt. + ! it is at least 1 + END_DOC + dim_list_virt_orb_kpts = max(maxval(n_virt_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_del_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_del. + ! it is at least 1 + END_DOC + dim_list_del_orb_kpts = max(maxval(n_del_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, n_core_inact_act_orb_kpts, (kpt_num) ] + implicit none + BEGIN_DOC + ! Number of core inactive and active MOs + END_DOC + integer :: k + do k=1,kpt_num + n_core_inact_act_orb_kpts(k) = (n_core_orb_kpts(k) + n_inact_orb_kpts(k) + n_act_orb_kpts(k)) + enddo +END_PROVIDER + + +!todo: finish below for kpts +! +! BEGIN_PROVIDER [ integer(bit_kind), core_bitmask , (N_int,2) ] +! implicit none +! BEGIN_DOC +! ! Bitmask identifying the core MOs +! END_DOC +! core_bitmask = 0_bit_kind +! if(n_core_orb > 0)then +! call list_to_bitstring( core_bitmask(1,1), list_core, n_core_orb, N_int) +! call list_to_bitstring( core_bitmask(1,2), list_core, n_core_orb, N_int) +! endif +! END_PROVIDER +! +! BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask, (N_int,2) ] +! implicit none +! BEGIN_DOC +! ! Bitmask identifying the inactive MOs +! END_DOC +! inact_bitmask = 0_bit_kind +! if(n_inact_orb > 0)then +! call list_to_bitstring( inact_bitmask(1,1), list_inact, n_inact_orb, N_int) +! call list_to_bitstring( inact_bitmask(1,2), list_inact, n_inact_orb, N_int) +! endif +! END_PROVIDER +! +! BEGIN_PROVIDER [ integer(bit_kind), act_bitmask , (N_int,2) ] +! implicit none +! BEGIN_DOC +! ! Bitmask identifying the active MOs +! END_DOC +! act_bitmask = 0_bit_kind +! if(n_act_orb > 0)then +! call list_to_bitstring( act_bitmask(1,1), list_act, n_act_orb, N_int) +! call list_to_bitstring( act_bitmask(1,2), list_act, n_act_orb, N_int) +! endif +! END_PROVIDER +! +! BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask , (N_int,2) ] +! implicit none +! BEGIN_DOC +! ! Bitmask identifying the virtual MOs +! END_DOC +! virt_bitmask = 0_bit_kind +! if(n_virt_orb > 0)then +! call list_to_bitstring( virt_bitmask(1,1), list_virt, n_virt_orb, N_int) +! call list_to_bitstring( virt_bitmask(1,2), list_virt, n_virt_orb, N_int) +! endif +! END_PROVIDER +! +! BEGIN_PROVIDER [ integer(bit_kind), del_bitmask , (N_int,2) ] +! implicit none +! BEGIN_DOC +! ! Bitmask identifying the deleted MOs +! END_DOC +! +! del_bitmask = 0_bit_kind +! +! if(n_del_orb > 0)then +! call list_to_bitstring( del_bitmask(1,1), list_del, n_del_orb, N_int) +! call list_to_bitstring( del_bitmask(1,2), list_del, n_del_orb, N_int) +! endif +! +! END_PROVIDER +! +! +! +! +! +! BEGIN_PROVIDER [ integer, list_core , (dim_list_core_orb) ] +!&BEGIN_PROVIDER [ integer, list_core_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of MO indices which are in the core. +! END_DOC +! integer :: i, n +! list_core = 0 +! list_core_reverse = 0 +! +! n=0 +! do i = 1, mo_num +! if(mo_class(i) == 'Core')then +! n += 1 +! list_core(n) = i +! list_core_reverse(i) = n +! endif +! enddo +! print *, 'Core MOs:' +! print *, list_core(1:n_core_orb) +! +!END_PROVIDER +! +! BEGIN_PROVIDER [ integer, list_inact , (dim_list_inact_orb) ] +!&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of MO indices which are inactive. +! END_DOC +! integer :: i, n +! list_inact = 0 +! list_inact_reverse = 0 +! +! n=0 +! do i = 1, mo_num +! if (mo_class(i) == 'Inactive')then +! n += 1 +! list_inact(n) = i +! list_inact_reverse(i) = n +! endif +! enddo +! print *, 'Inactive MOs:' +! print *, list_inact(1:n_inact_orb) +! +!END_PROVIDER +! +! BEGIN_PROVIDER [ integer, list_virt , (dim_list_virt_orb) ] +!&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of MO indices which are virtual +! END_DOC +! integer :: i, n +! list_virt = 0 +! list_virt_reverse = 0 +! +! n=0 +! do i = 1, mo_num +! if (mo_class(i) == 'Virtual')then +! n += 1 +! list_virt(n) = i +! list_virt_reverse(i) = n +! endif +! enddo +! print *, 'Virtual MOs:' +! print *, list_virt(1:n_virt_orb) +! +!END_PROVIDER +! +! BEGIN_PROVIDER [ integer, list_del , (dim_list_del_orb) ] +!&BEGIN_PROVIDER [ integer, list_del_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of MO indices which are deleted. +! END_DOC +! integer :: i, n +! list_del = 0 +! list_del_reverse = 0 +! +! n=0 +! do i = 1, mo_num +! if (mo_class(i) == 'Deleted')then +! n += 1 +! list_del(n) = i +! list_del_reverse(i) = n +! endif +! enddo +! print *, 'Deleted MOs:' +! print *, list_del(1:n_del_orb) +! +!END_PROVIDER +! +! BEGIN_PROVIDER [ integer, list_act , (dim_list_act_orb) ] +!&BEGIN_PROVIDER [ integer, list_act_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of MO indices which are in the active. +! END_DOC +! integer :: i, n +! list_act = 0 +! list_act_reverse = 0 +! +! n=0 +! do i = 1, mo_num +! if (mo_class(i) == 'Active')then +! n += 1 +! list_act(n) = i +! list_act_reverse(i) = n +! endif +! enddo +! print *, 'Active MOs:' +! print *, list_act(1:n_act_orb) +! +!END_PROVIDER +! +! +! +! BEGIN_PROVIDER [ integer, list_core_inact , (dim_list_core_inact_orb) ] +!&BEGIN_PROVIDER [ integer, list_core_inact_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of indices of the core and inactive MOs +! END_DOC +! integer :: i,itmp +! call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), list_core_inact, itmp, N_int) +! list_core_inact_reverse = 0 +! ASSERT (itmp == n_core_inact_orb) +! do i = 1, n_core_inact_orb +! list_core_inact_reverse(list_core_inact(i)) = i +! enddo +! print *, 'Core and Inactive MOs:' +! print *, list_core_inact(1:n_core_inact_orb) +!END_PROVIDER +! +! +! BEGIN_PROVIDER [ integer, list_core_inact_act , (n_core_inact_act_orb) ] +!&BEGIN_PROVIDER [ integer, list_core_inact_act_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of indices of the core inactive and active MOs +! END_DOC +! integer :: i,itmp +! call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), list_core_inact_act, itmp, N_int) +! list_core_inact_act_reverse = 0 +! ASSERT (itmp == n_core_inact_act_orb) +! do i = 1, n_core_inact_act_orb +! list_core_inact_act_reverse(list_core_inact_act(i)) = i +! enddo +! print *, 'Core, Inactive and Active MOs:' +! print *, list_core_inact_act(1:n_core_inact_act_orb) +!END_PROVIDER +! +! +! BEGIN_PROVIDER [ integer, list_inact_act , (n_inact_act_orb) ] +!&BEGIN_PROVIDER [ integer, list_inact_act_reverse, (mo_num) ] +! implicit none +! BEGIN_DOC +! ! List of indices of the inactive and active MOs +! END_DOC +! integer :: i,itmp +! call bitstring_to_list(reunion_of_inact_act_bitmask(1,1), list_inact_act, itmp, N_int) +! list_inact_act_reverse = 0 +! ASSERT (itmp == n_inact_act_orb) +! do i = 1, n_inact_act_orb +! list_inact_act_reverse(list_inact_act(i)) = i +! enddo +! print *, 'Inactive and Active MOs:' +! print *, list_inact_act(1:n_inact_act_orb) +!END_PROVIDER +! diff --git a/src/hartree_fock/scf.irp.f b/src/hartree_fock/scf.irp.f index 2b93a1df..9d1671ec 100644 --- a/src/hartree_fock/scf.irp.f +++ b/src/hartree_fock/scf.irp.f @@ -46,21 +46,25 @@ subroutine create_guess logical :: exists PROVIDE ezfio_filename if (is_complex) then - call ezfio_has_mo_basis_mo_coef_complex(exists) +! call ezfio_has_mo_basis_mo_coef_complex(exists) + call ezfio_has_mo_basis_mo_coef_kpts(exists) else call ezfio_has_mo_basis_mo_coef(exists) endif if (.not.exists) then if (mo_guess_type == "HCore") then if (is_complex) then - mo_coef_complex = ao_ortho_lowdin_coef_complex - TOUCH mo_coef_complex + !mo_coef_complex = ao_ortho_lowdin_coef_complex + mo_coef_kpts = ao_ortho_lowdin_coef_kpts + TOUCH mo_coef_kpts mo_label = 'Guess' - call mo_as_eigvectors_of_mo_matrix_complex(mo_one_e_integrals_complex, & - size(mo_one_e_integrals_complex,1), & - size(mo_one_e_integrals_complex,2), & + !call mo_as_eigvectors_of_mo_matrix_complex(mo_one_e_integrals_kpts, & + call mo_as_eigvectors_of_mo_matrix_kpts(mo_one_e_integrals_kpts, & + size(mo_one_e_integrals_kpts,1), & + size(mo_one_e_integrals_kpts,2), & + size(mo_one_e_integrals_kpts,3), & mo_label,1,.false.) - SOFT_TOUCH mo_coef_complex mo_label + SOFT_TOUCH mo_coef_kpts mo_label else mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef @@ -73,7 +77,8 @@ subroutine create_guess endif else if (mo_guess_type == "Huckel") then if (is_complex) then - call huckel_guess_complex + !call huckel_guess_complex + call huckel_guess_kpts else call huckel_guess endif diff --git a/src/mo_basis/mos_cplx.irp.f b/src/mo_basis/mos_cplx.irp.f index a4f3f9ed..6b7e14c7 100644 --- a/src/mo_basis/mos_cplx.irp.f +++ b/src/mo_basis/mos_cplx.irp.f @@ -93,6 +93,7 @@ BEGIN_PROVIDER [ complex*16, mo_coef_complex_kpts, (ao_num_per_kpt, mo_num_per_k END_PROVIDER + BEGIN_PROVIDER [ complex*16, mo_coef_transp_complex, (mo_num,ao_num) ] &BEGIN_PROVIDER [ complex*16, mo_coef_transp_complex_conjg, (mo_num,ao_num) ] implicit none @@ -198,3 +199,276 @@ subroutine ao_ortho_cano_to_ao_cplx(A_ao,LDA_ao,A,LDA) deallocate(T) end +!============================================! +! ! +! kpts ! +! ! +!============================================! + + +BEGIN_PROVIDER [ complex*16, mo_coef_kpts, (ao_num_per_kpt, mo_num_per_kpt, kpt_num) ] + implicit none + BEGIN_DOC + ! Molecular orbital coefficients on |AO| basis set + ! + ! mo_coef_kpts(i,j,k) = coefficient of the i-th |AO| on the jth |MO| in kth kpt + ! + ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) + END_DOC + integer :: i, j, k + logical :: exists + PROVIDE ezfio_filename + + if (mpi_master) then + ! Coefs + call ezfio_has_mo_basis_mo_coef_complex(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_kpts with MPI' + endif + IRP_ENDIF + + if (exists) then + if (mpi_master) then + call ezfio_get_mo_basis_mo_coef_kpts(mo_coef_kpts) + write(*,*) 'Read mo_coef_kpts' + endif + IRP_IF MPI + call MPI_BCAST( mo_coef_kpts, kpt_num*mo_num_per_kpt*ao_num_per_kpt, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_kpts with MPI' + endif + IRP_ENDIF + else + ! Orthonormalized AO basis + + do k=1,kpt_num + do i=1,mo_num_per_kpt + do j=1,ao_num_per_kpt + mo_coef_kpts(j,i,k) = ao_ortho_canonical_coef_kpts(j,i,k) + enddo + enddo + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, mo_coef_in_ao_ortho_basis_kpts, (ao_num_per_kpt, mo_num_per_kpt, kpt_num) ] + implicit none + BEGIN_DOC + ! |MO| coefficients in orthogonalized |AO| basis + ! + ! $C^{-1}.C_{mo}$ + END_DOC + integer :: k + do k=1,kpt_num + + call zgemm('N','N',ao_num_per_kpt,mo_num_per_kpt,ao_num_per_kpt,(1.d0,0.d0), & + ao_ortho_cano_coef_inv_kpts(:,:,k), size(ao_ortho_cano_coef_inv_kpts,1),& + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), (0.d0,0.d0), & + mo_coef_in_ao_ortho_basis_kpts(:,:,k), size(mo_coef_in_ao_ortho_basis_kpts,1)) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, mo_coef_transp_kpts, (mo_num_per_kpt,ao_num_per_kpt,kpt_num) ] +&BEGIN_PROVIDER [ complex*16, mo_coef_transp_kpts_conjg, (mo_num_per_kpt,ao_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! |MO| coefficients on |AO| basis set + END_DOC + integer :: i, j, k + + do k=1,kpt_num + do j=1,ao_num_per_kpt + do i=1,mo_num_per_kpt + mo_coef_transp_kpts(i,j,k) = mo_coef_kpts(j,i,k) + mo_coef_transp_kpts_conjg(i,j,k) = dconjg(mo_coef_kpts(j,i,k)) + enddo + enddo + enddo + +END_PROVIDER + +subroutine ao_to_mo_kpts(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + !todo: check this + BEGIN_DOC + ! Transform A from the AO basis to the MO basis + ! where A is complex in the AO basis + ! + ! C^\dagger.A_ao.C + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num_per_kpt,kpt_num) + complex*16, intent(out) :: A_mo(LDA_mo,mo_num_per_kpt,kpt_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num_per_kpt,mo_num_per_kpt) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + integer :: k + + do k=1,kpt_num + call zgemm('N','N', ao_num_per_kpt, mo_num_per_kpt, ao_num_per_kpt, & + (1.d0,0.d0), A_ao,LDA_ao, & + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('C','N', mo_num_per_kpt, mo_num_per_kpt, ao_num_per_kpt, & + (1.d0,0.d0), mo_coef_kpts(:,:,k),size(mo_coef_kpts,1), & + T, ao_num_per_kpt, & + (0.d0,0.d0), A_mo(:,:,k), size(A_mo,1)) + enddo + + deallocate(T) +end + +subroutine ao_to_mo_noconjg_kpts(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + BEGIN_DOC + ! Transform A from the AO basis to the MO basis + ! where A is complex in the AO basis + ! + ! C^T.A_ao.C + ! needed for 4idx tranform in four_idx_novvvv + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num_per_kpt,kpt_num) + complex*16, intent(out) :: A_mo(LDA_mo,mo_num_per_kpt,kpt_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num_per_kpt,mo_num_per_kpt) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + integer :: k + do k=1,kpt_num + call zgemm('N','N', ao_num_per_kpt, mo_num_per_kpt, ao_num_per_kpt, & + (1.d0,0.d0), A_ao,LDA_ao, & + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('T','N', mo_num_per_kpt, mo_num_per_kpt, ao_num_per_kpt, & + (1.d0,0.d0), mo_coef_kpts(:,:,k),size(mo_coef_kpts,1), & + T, ao_num_per_kpt, & + (0.d0,0.d0), A_mo(:,:,k), size(A_mo,1)) + enddo + deallocate(T) +end + + +subroutine ao_ortho_cano_to_ao_kpts(A_ao,LDA_ao,A,LDA) + implicit none + !todo: check this; no longer using assumed-size arrays + BEGIN_DOC + ! Transform A from the |AO| basis to the orthogonal |AO| basis + ! + ! $C^{-1}.A_{ao}.C^{\dagger-1}$ + END_DOC + integer, intent(in) :: LDA_ao,LDA + complex*16, intent(in) :: A_ao(LDA_ao,ao_num_per_kpt,kpt_num) + complex*16, intent(out) :: A(LDA,ao_num_per_kpt,kpt_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num_per_kpt,ao_num_per_kpt) ) + + integer :: k + do k=1,kpt_num + call zgemm('C','N', ao_num_per_kpt, ao_num_per_kpt, ao_num_per_kpt, & + (1.d0,0.d0), & + ao_ortho_cano_coef_inv_kpts(:,:,k), size(ao_ortho_cano_coef_inv_kpts,1),& + A_ao(:,:,k),size(A_ao,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num_per_kpt, ao_num_per_kpt, ao_num_per_kpt, (1.d0,0.d0), & + T, size(T,1), & + ao_ortho_cano_coef_inv_kpts(:,:,k),size(ao_ortho_cano_coef_inv_kpts,1),& + (0.d0,0.d0), A(:,:,k), size(A,1)) + enddo + + deallocate(T) +end + + +!============================================! +! ! +! elec kpts ! +! ! +!============================================! + + BEGIN_PROVIDER [ integer, elec_alpha_num_kpts, (kpt_num) ] +&BEGIN_PROVIDER [ integer, elec_beta_num_kpts, (kpt_num) ] + !todo: reorder? if not integer multiple, use some list of kpts to determine filling order + implicit none + + integer :: i,k,kpt + + PROVIDE elec_alpha_num elec_beta_num + + do k=1,kpt_num + elec_alpha_num_kpts(k) = 0 + elec_beta_num_kpts(k) = 0 + enddo + kpt=1 + do i=1,elec_beta_num + elec_alpha_num_kpts(kpt) += 1 + elec_beta_num_kpts(kpt) += 1 + kpt += 1 + if (kpt > kpt_num) then + kpt = 1 + endif + enddo + do i=elec_beta_num+1,elec_alpha_num + elec_alpha_num_kpts(kpt) += 1 + kpt += 1 + if (kpt > kpt_num) then + kpt = 1 + endif + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_occ_kpts, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! |MO| occupation numbers + END_DOC + PROVIDE ezfio_filename elec_beta_num_kpts elec_alpha_num_kpts + if (mpi_master) then + logical :: exists + call ezfio_has_mo_basis_mo_occ_kpts(exists) + if (exists) then + call ezfio_get_mo_basis_mo_occ_kpts(mo_occ_kpts) + else + mo_occ_kpts = 0.d0 + integer :: i,k + do k=1,kpt_num + do i=1,elec_beta_num_kpts(k) + mo_occ_kpts(i,k) = 2.d0 + enddo + do i=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) + mo_occ_kpts(i,k) = 1.d0 + enddo + enddo + endif + write(*,*) 'Read mo_occ_kpts' + 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( mo_occ_kpts, mo_num_per_kpt*kpt_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_occ_kpts with MPI' + endif + IRP_ENDIF + +END_PROVIDER diff --git a/src/mo_basis/utils_cplx.irp.f b/src/mo_basis/utils_cplx.irp.f index a967cec4..58230cdd 100644 --- a/src/mo_basis/utils_cplx.irp.f +++ b/src/mo_basis/utils_cplx.irp.f @@ -245,4 +245,271 @@ subroutine mo_coef_new_as_svd_vectors_of_mo_matrix_eig_complex(matrix,lda,m,n,mo end +!============================================! +! ! +! kpts ! +! ! +!============================================! + +subroutine mo_as_eigvectors_of_mo_matrix_kpts(matrix,n,m,nk,label,sign,output) + !TODO: test this + implicit none + integer,intent(in) :: n,m,nk, sign + character*(64), intent(in) :: label + complex*16, intent(in) :: matrix(n,m,nk) + logical, intent(in) :: output + + integer :: i,j,k + double precision, allocatable :: eigvalues(:) + complex*16, allocatable :: mo_coef_new(:,:), R(:,:), A(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R + + call write_time(6) + if (m /= mo_num_per_kpt) then + print *, irp_here, ': Error : m/= mo_num_per_kpt' + stop 1 + endif + if (nk /= kpt_num) then + print *, irp_here, ': Error : nk/= kpt_num' + stop 1 + endif + allocate(A(n,m),R(n,m),mo_coef_new(ao_num_per_kpt,m),eigvalues(m)) + do k=1,nk + if (sign == -1) then + do j=1,m + do i=1,n + A(i,j) = -matrix(i,j,k) + enddo + enddo + else + do j=1,m + do i=1,n + A(i,j) = matrix(i,j,k) + enddo + enddo + endif + mo_coef_new = mo_coef_kpts(:,:,k) + + call lapack_diag_complex(eigvalues,R,A,n,m) + if (sign == -1) then + do i=1,m + eigvalues(i) = -eigvalues(i) + enddo + endif + if (output) then + do i=1,m + write (6,'(2(I8),1X,F16.10)') k,i,eigvalues(i) + enddo + write (6,'(A)') '======== ================' + write (6,'(A)') '' + write (6,'(A)') 'Fock Matrix' + write (6,'(A)') '-----------' + do i=1,n + write(*,'(200(E24.15))') A(i,:) + enddo + endif + + call zgemm('N','N',ao_num_per_kpt,m,m,(1.d0,0.d0), & + mo_coef_new,size(mo_coef_new,1),R,size(R,1),(0.d0,0.d0), & + mo_coef_kpts(:,:,k),size(mo_coef_kpts,1)) + enddo + deallocate(A,mo_coef_new,R,eigvalues) + call write_time(6) + + mo_label = label + if (output) then + write (6,'(A)') 'MOs are now **'//trim(label)//'**' + write (6,'(A)') '' + write (6,'(A)') 'Eigenvalues' + write (6,'(A)') '-----------' + write (6,'(A)') '' + write (6,'(A)') '======== ================' + endif +end + +subroutine mo_as_svd_vectors_of_mo_matrix_kpts(matrix,lda,m,n,label) + !TODO: implement + print *, irp_here, ' not implemented for kpts' + stop 1 + implicit none + integer,intent(in) :: lda,m,n + character*(64), intent(in) :: label + complex*16, intent(in) :: matrix(lda,n) + + integer :: i,j + double precision :: accu + double precision, allocatable :: D(:) + complex*16, allocatable :: mo_coef_new(:,:), U(:,:), A(:,:), Vt(:,:) + !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' + stop 1 + endif + + allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num,m),D(m),Vt(lda,n)) + + do j=1,n + do i=1,m + A(i,j) = matrix(i,j) + enddo + enddo + mo_coef_new = mo_coef_complex + + call svd_complex(A,lda,U,lda,D,Vt,lda,m,n) + + write (6,'(A)') 'MOs are now **'//trim(label)//'**' + write (6,'(A)') '' + 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 + 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)) + deallocate(A,mo_coef_new,U,Vt,D) + call write_time(6) + + mo_label = label +end + + +subroutine mo_as_svd_vectors_of_mo_matrix_eig_kpts(matrix,lda,m,n,eig,label) + !TODO: implement + print *, irp_here, ' not implemented for kpts' + stop 1 + implicit none + integer,intent(in) :: lda,m,n + character*(64), intent(in) :: label + complex*16, intent(in) :: matrix(lda,n) + double precision, intent(out) :: eig(m) + + integer :: i,j + 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' + stop 1 + endif + + allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num,m),D(m),Vt(lda,n)) + + do j=1,n + do i=1,m + A(i,j) = matrix(i,j) + enddo + enddo + mo_coef_new = mo_coef_complex + + call svd_complex(A,lda,U,lda,D,Vt,lda,m,n) + + write (6,'(A)') 'MOs are now **'//trim(label)//'**' + write (6,'(A)') '' + 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 + 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 + +end + + +subroutine mo_coef_new_as_svd_vectors_of_mo_matrix_eig_kpts(matrix,lda,m,n,mo_coef_before,eig,mo_coef_new) + !TODO: implement + print *, irp_here, ' not implemented for kpts' + stop 1 + implicit none + BEGIN_DOC +! You enter with matrix in the MO basis defined with the mo_coef_before. +! +! You SVD the matrix and set the eigenvectors as mo_coef_new ordered by increasing singular values + END_DOC + integer,intent(in) :: lda,m,n + complex*16, intent(in) :: matrix(lda,n),mo_coef_before(ao_num,m) + double precision, intent(out) :: eig(m) + complex*16, intent(out) :: mo_coef_new(ao_num,m) + + integer :: i,j + double precision :: accu + double precision, allocatable :: D(:) + complex*16, allocatable :: mo_coef_tmp(:,:), U(:,:), A(:,:), Vt(:,:), work(:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, Vt, A + + call write_time(6) + if (m /= mo_num) then + print *, irp_here, ': Error : m/= mo_num' + stop 1 + endif + + allocate(A(lda,n),U(lda,n),D(m),Vt(lda,n),mo_coef_tmp(ao_num,mo_num)) + + do j=1,n + do i=1,m + A(i,j) = matrix(i,j) + enddo + enddo + mo_coef_tmp = mo_coef_before + + call svd_complex(A,lda,U,lda,D,Vt,lda,m,n) + + write (6,'(A)') '' + 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 + enddo + write (6,'(A)') '======== ================ ================' + write (6,'(A)') '' + + call zgemm('N','N',ao_num,m,m,(1.d0,0.d0),mo_coef_tmp,size(mo_coef_new,1),U,size(U,1),(0.d0,0.d0),mo_coef_new,size(mo_coef_new,1)) + + do i=1,m + eig(i) = D(i) + enddo + + deallocate(A,U,Vt,D,mo_coef_tmp) + call write_time(6) + +end diff --git a/src/mo_guess/mo_ortho_lowdin_cplx.irp.f b/src/mo_guess/mo_ortho_lowdin_cplx.irp.f index 5e1dacbe..3a2750cd 100644 --- a/src/mo_guess/mo_ortho_lowdin_cplx.irp.f +++ b/src/mo_guess/mo_ortho_lowdin_cplx.irp.f @@ -46,3 +46,64 @@ BEGIN_PROVIDER [complex*16, ao_ortho_lowdin_overlap_complex, (ao_num,ao_num)] enddo enddo END_PROVIDER + +!============================================! +! ! +! kpts ! +! ! +!============================================! + +BEGIN_PROVIDER [complex*16, ao_ortho_lowdin_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num)] + implicit none + BEGIN_DOC +! matrix of the coefficients of the mos generated by the +! orthonormalization by the S^{-1/2} canonical transformation of the aos +! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital + END_DOC + integer :: i,j,k,l + complex*16, allocatable :: tmp_matrix(:,:) + allocate (tmp_matrix(ao_num,ao_num)) + do k=1,kpt_num + tmp_matrix(:,:) = (0.d0,0.d0) + do j=1, ao_num + tmp_matrix(j,j) = (1.d0,0.d0) + enddo + call ortho_lowdin_complex(ao_overlap_kpts(:,:,k),ao_num_per_kpt,ao_num_per_kpt,tmp_matrix,ao_num_per_kpt,ao_num_per_kpt) + do i=1, ao_num_per_kpt + do j=1, ao_num_per_kpt + ao_ortho_lowdin_coef_kpts(j,i,k) = tmp_matrix(i,j) + enddo + enddo + enddo + deallocate(tmp_matrix) +END_PROVIDER + +BEGIN_PROVIDER [complex*16, ao_ortho_lowdin_overlap_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num)] + implicit none + BEGIN_DOC +! overlap matrix of the ao_ortho_lowdin +! supposed to be the Identity + END_DOC + integer :: i,j,k,l,kk + complex*16 :: c + do kk=1,kpt_num + do j=1, ao_num_per_kpt + do i=1, ao_num_per_kpt + ao_ortho_lowdin_overlap_kpts(i,j,kk) = (0.d0,0.d0) + enddo + enddo + enddo + do kk=1,kpt_num + do k=1, ao_num_per_kpt + do j=1, ao_num_per_kpt + c = (0.d0,0.d0) + do l=1, ao_num_per_kpt + c += dconjg(ao_ortho_lowdin_coef_kpts(j,l,kk)) * ao_overlap_kpts(k,l,kk) + enddo + do i=1, ao_num_per_kpt + ao_ortho_lowdin_overlap_kpts(i,j,kk) += ao_ortho_lowdin_coef_kpts(i,k,kk) * c + enddo + enddo + enddo + enddo +END_PROVIDER diff --git a/src/scf_utils/diagonalize_fock_cplx.irp.f b/src/scf_utils/diagonalize_fock_cplx.irp.f index 645dbcf9..de38c767 100644 --- a/src/scf_utils/diagonalize_fock_cplx.irp.f +++ b/src/scf_utils/diagonalize_fock_cplx.irp.f @@ -51,3 +51,65 @@ BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_complex, (ao_num,mo_num END_PROVIDER +!============================================! +! ! +! kpts ! +! ! +!============================================! +BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts, (ao_num_per_kpt,mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! Eigenvectors of the Fock matrix in the |MO| basis obtained with level shift. + END_DOC + + integer :: i,j,k + integer :: n + complex*16, allocatable :: F(:,:) + double precision, allocatable :: diag(:) + + + allocate( F(mo_num_per_kpt,mo_num_per_kpt) ) + allocate (diag(mo_num_per_kpt) ) + + do k=1,kpt_num + do j=1,mo_num + do i=1,mo_num + !F(i,j) = fock_matrix_mo_complex(i,j) + F(i,j) = fock_matrix_mo_kpts(i,j,k) + enddo + enddo + + if(frozen_orb_scf)then + integer :: iorb,jorb + !todo: core/act per kpt + do i = 1, n_core_orb + iorb = list_core(i) + do j = 1, n_act_orb + jorb = list_act(j) + F(iorb,jorb) = (0.d0,0.d0) + F(jorb,iorb) = (0.d0,0.d0) + enddo + enddo + endif + + ! Insert level shift here + !todo: elec per kpt + do i = elec_beta_num_per_kpt(k)+1, elec_alpha_num_per_kpt(k) + F(i,i) += 0.5d0*level_shift + enddo + + do i = elec_alpha_num_per_kpt(k)+1, mo_num_per_kpt + F(i,i) += level_shift + enddo + + n = mo_num_per_kpt + call lapack_diagd_diag_in_place_complex(diag,F,n,n) + + call zgemm('N','N',ao_num_per_kpt,mo_num_per_kpt,mo_num_per_kpt, (1.d0,0.d0), & + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), F, size(F,1), & + (0.d0,0.d0), eigenvectors_Fock_matrix_mo_kpts(:,:,k), size(eigenvectors_Fock_matrix_mo_kpts,1)) + enddo + deallocate(F, diag) + + +END_PROVIDER diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index 61d23467..94508570 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -359,3 +359,377 @@ END_PROVIDER enddo END_PROVIDER + +!============================================! +! ! +! kpts ! +! ! +!============================================! + + BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo_kpts, (mo_num_per_kpt,kpt_num)] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is :: + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + integer :: i,j,n,k + !todo: fix for kpts? (okay for simple cases) + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_mo_kpts = Fock_matrix_mo_alpha_kpts + else + do k=1,kpt_num + do j=1,elec_beta_num_kpts(k) + ! F-K + do i=1,elec_beta_num_kpts(k) !CC + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k))& + - (Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + ! F+K/2 + do i=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) !CA + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k))& + + 0.5d0*(Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + ! F + do i=elec_alpha_num_kpts(k)+1, mo_num_per_kpt !CV + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k)) + enddo + enddo + + do j=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) + ! F+K/2 + do i=1,elec_beta_num_kpts(k) !AC + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k))& + + 0.5d0*(Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + ! F + do i=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) !AA + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k)) + enddo + ! F-K/2 + do i=elec_alpha_num_kpts(k)+1, mo_num_per_kpt !AV + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k))& + - 0.5d0*(Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + enddo + + do j=elec_alpha_num_kpts(k)+1, mo_num_per_kpt + ! F + do i=1,elec_beta_num_kpts(k) !VC + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k)) + enddo + ! F-K/2 + do i=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) !VA + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k))& + - 0.5d0*(Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + ! F+K + do i=elec_alpha_num_kpts(k)+1,mo_num_per_kpt !VV + Fock_matrix_mo_kpts(i,j,k) = 0.5d0*(Fock_matrix_mo_alpha_kpts(i,j,k)+Fock_matrix_mo_beta_kpts(i,j,k)) & + + (Fock_matrix_mo_beta_kpts(i,j,k) - Fock_matrix_mo_alpha_kpts(i,j,k)) + enddo + enddo + enddo + + endif + do k=1,kpt_num + do i = 1, mo_num_per_kpt + Fock_matrix_diag_mo_kpts(i,k) = dble(Fock_matrix_mo_kpts(i,i,k)) + if (dabs(dimag(Fock_matrix_mo_kpts(i,i,k))) .gt. 1.0d-12) then + !stop 'diagonal elements of Fock matrix should be real' + print *, 'diagonal elements of Fock matrix should be real',i,Fock_matrix_mo_kpts(i,i,k) + !stop -1 + endif + enddo + enddo + + + if(frozen_orb_scf)then + integer :: iorb,jorb + do k=1,kpt_num + ! for tags: list_core, n_core_orb, n_act_orb, list_act + do i = 1, n_core_orb_kpts(k) + iorb = list_core_kpts(i,k) + do j = 1, n_act_orb_kpts(k) + jorb = list_act_kpts(j,k) + fock_matrix_mo_kpts(iorb,jorb,k) = (0.d0,0.d0) + fock_matrix_mo_kpts(jorb,iorb,k) = (0.d0,0.d0) + enddo + enddo + enddo + endif + +END_PROVIDER + + + +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_alpha_complex,size(Fock_matrix_ao_alpha_complex,1), & + Fock_matrix_mo_alpha_complex,size(Fock_matrix_mo_alpha_complex,1)) +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_beta_complex,size(Fock_matrix_ao_beta_complex,1), & + Fock_matrix_mo_beta_complex,size(Fock_matrix_mo_beta_complex,1)) +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if(frozen_orb_scf)then + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + else + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(i,j) + enddo + enddo + else + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + endif + endif +END_PROVIDER + + + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_complex, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_complex , (ao_num, ao_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha and Beta Fock matrices in AO basis set + END_DOC + !TODO: finish implementing this: see complex qp1 (different mapping) + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + complex*16 :: integral, c0 + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:) + + ao_two_e_integral_alpha_complex = (0.d0,0.d0) + ao_two_e_integral_beta_complex = (0.d0,0.d0) + PROVIDE ao_two_e_integrals_in_map + + integer(omp_lock_kind) :: lck(ao_num) + integer(map_size_kind) :: i8 + integer :: ii(4), jj(4), kk(4), ll(4), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) + integer(key_kind) :: key1 + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & + !$OMP SCF_density_matrix_ao_beta_complex, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_num)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_1(ii,jj,kk,ll,key1) + ! i<=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + + if (shiftl(key1,1)==keys(k1)) then !imaginary part (even) + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = values(k1) + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_complex += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & + !$OMP SCF_density_matrix_ao_beta_complex, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) + + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_num)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map_2%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map_2,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_2(ii,jj,kk,ll,key1) + ! i>=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + if (shiftl(key1,1)==keys(k1)) then !imaginary part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + integral = values(k1) + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_complex += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_alpha_complex, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_beta_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Alpha Fock matrix in AO basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_alpha_complex(i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_alpha_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_beta_complex (i,j) + enddo + enddo + +END_PROVIDER diff --git a/src/scf_utils/huckel_cplx.irp.f b/src/scf_utils/huckel_cplx.irp.f index d6da7ffb..a41ed831 100644 --- a/src/scf_utils/huckel_cplx.irp.f +++ b/src/scf_utils/huckel_cplx.irp.f @@ -40,3 +40,52 @@ subroutine huckel_guess_complex deallocate(A) end +!============================================! +! ! +! kpts ! +! ! +!============================================! +subroutine huckel_guess_kpts + implicit none + BEGIN_DOC +! Build the MOs using the extended Huckel model + END_DOC + integer :: i,j,k + double precision :: accu + double precision :: c + character*(64) :: label + complex*16, allocatable :: A(:,:) + label = "Guess" + c = 0.5d0 * 1.75d0 + + allocate (A(ao_num, ao_num)) + do k=1,kpt_num + A = (0.d0,0.d0) + do j=1,ao_num_per_kpt + do i=1,ao_num_per_kpt + A(i,j) = c * ao_overlap_kpts(i,j,k) * (ao_one_e_integrals_diag_kpts(i,k) + ao_one_e_integrals_diag_kpts(j,k)) + enddo + A(j,j) = ao_one_e_integrals_diag_kpts(j,k) + dble(ao_two_e_integral_alpha_kpts(j,j,k)) + if (dabs(dimag(ao_two_e_integral_alpha_kpts(j,j,k))) .gt. 1.0d-10) then + stop 'diagonal elements of ao_bi_elec_integral_alpha should be real' + endif + enddo + +! Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) +! Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) + call zlacpy('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), & + Fock_matrix_ao_alpha_kpts(:,:,k), size(Fock_matrix_ao_alpha_kpts,1)) + call zlacpy('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), & + Fock_matrix_ao_beta_kpts(:,:,k), size(Fock_matrix_ao_beta_kpts, 1)) + enddo + +! TOUCH mo_coef + + !TOUCH fock_matrix_ao_alpha_complex fock_matrix_ao_beta_kpts + TOUCH fock_matrix_ao_alpha_kpts fock_matrix_ao_beta_kpts + mo_coef_kpts = eigenvectors_fock_matrix_mo_kpts + SOFT_TOUCH mo_coef_complex + call save_mos + deallocate(A) + +end diff --git a/src/utils_complex/qp2-pbc-diff.txt b/src/utils_complex/qp2-pbc-diff.txt index 7cadd3c4..e7345240 100644 --- a/src/utils_complex/qp2-pbc-diff.txt +++ b/src/utils_complex/qp2-pbc-diff.txt @@ -1,7 +1,9 @@ todo: -change everything to be blocked by kpt + change everything to be blocked by kpt + elec_alpha_num_per_kpt (maybe add to mo_basis?) + bitmasks per kpt? (or at least occ/act/virt num and list) ------------------------------------------------------------------------------------- old: