From d44a22f3d8167bc807c34794e11e8f6493b89d50 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 12 Mar 2020 16:09:00 -0500 Subject: [PATCH] fix bitmask for kpoint ordered MOs --- src/bitmask/bitmasks.irp.f | 32 ++++++++++++-- src/scf_utils/fock_matrix_cplx.irp.f | 2 +- .../scf_density_matrix_ao_cplx.irp.f | 42 +++++++++++++++++-- 3 files changed, 67 insertions(+), 9 deletions(-) diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index 91617397..e86d747d 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -78,14 +78,38 @@ BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)] END_DOC integer :: i,j,n integer :: occ(elec_alpha_num) + integer :: occb(elec_beta_num) HF_bitmask = 0_bit_kind - do i=1,elec_alpha_num - occ(i) = i - enddo + if (is_complex) then + integer :: kpt,korb + kpt=1 + korb=1 + do i=1,elec_beta_num + occ(i) = korb + (kpt-1) * ao_num_per_kpt + occb(i) = korb + (kpt-1) * ao_num_per_kpt + kpt += 1 + if (kpt > kpt_num) then + kpt = 1 + korb += 1 + endif + enddo + do i=elec_beta_num+1,elec_alpha_num + occ(i) = korb + (kpt-1) * ao_num_per_kpt + kpt += 1 + if (kpt > kpt_num) then + kpt = 1 + korb += 1 + endif + enddo + else + do i=1,elec_alpha_num + occ(i) = i + enddo + endif call list_to_bitstring( HF_bitmask(1,1), occ, elec_alpha_num, N_int) ! elec_alpha_num <= elec_beta_num, so occ is already OK. - call list_to_bitstring( HF_bitmask(1,2), occ, elec_beta_num, N_int) + call list_to_bitstring( HF_bitmask(1,2), occb, elec_beta_num, N_int) END_PROVIDER diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index 577fe5c2..61d23467 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -80,7 +80,7 @@ if (dabs(dimag(Fock_matrix_mo_complex(i,i))) .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_complex(i,i) - stop -1 + !stop -1 endif enddo diff --git a/src/scf_utils/scf_density_matrix_ao_cplx.irp.f b/src/scf_utils/scf_density_matrix_ao_cplx.irp.f index 6e22e209..a6a66863 100644 --- a/src/scf_utils/scf_density_matrix_ao_cplx.irp.f +++ b/src/scf_utils/scf_density_matrix_ao_cplx.irp.f @@ -3,12 +3,29 @@ BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_alpha_complex, (ao_num,ao_num BEGIN_DOC ! $C.C^t$ over $\alpha$ MOs END_DOC + + complex*16, allocatable :: mo_coef_alpha_tmp(:,:) + integer :: occ(N_int*bit_kind_size) + integer :: na, i + + call bitstring_to_list(hf_bitmask(1,1), occ, na, n_int) + allocate(mo_coef_alpha_tmp(ao_num,na)) + do i=1,na + mo_coef_alpha_tmp(:,i) = mo_coef_complex(:,occ(i)) + enddo + call zgemm('N','C',ao_num,ao_num,elec_alpha_num,(1.d0,0.d0), & - mo_coef_complex, size(mo_coef_complex,1), & - mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & + mo_coef_alpha_tmp, size(mo_coef_alpha_tmp,1), & + mo_coef_alpha_tmp, size(mo_coef_alpha_tmp,1), (0.d0,0.d0), & scf_density_matrix_ao_alpha_complex, size(scf_density_matrix_ao_alpha_complex,1)) + deallocate(mo_coef_alpha_tmp) + !call zgemm('N','C',ao_num,ao_num,elec_alpha_num,(1.d0,0.d0), & + ! mo_coef_complex, size(mo_coef_complex,1), & + ! mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & + ! scf_density_matrix_ao_alpha_complex, size(scf_density_matrix_ao_alpha_complex,1)) + END_PROVIDER BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_beta_complex, (ao_num,ao_num) ] @@ -17,11 +34,28 @@ BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_beta_complex, (ao_num,ao_num ! $C.C^t$ over $\beta$ MOs END_DOC + complex*16, allocatable :: mo_coef_beta_tmp(:,:) + integer :: occ(N_int*bit_kind_size) + integer :: nb, i + + call bitstring_to_list(hf_bitmask(1,2), occ, nb, n_int) + allocate(mo_coef_beta_tmp(ao_num,nb)) + do i=1,nb + mo_coef_beta_tmp(:,i) = mo_coef_complex(:,occ(i)) + enddo + + call zgemm('N','C',ao_num,ao_num,elec_beta_num,(1.d0,0.d0), & - mo_coef_complex, size(mo_coef_complex,1), & - mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & + mo_coef_beta_tmp, size(mo_coef_beta_tmp,1), & + mo_coef_beta_tmp, size(mo_coef_beta_tmp,1), (0.d0,0.d0), & scf_density_matrix_ao_beta_complex, size(scf_density_matrix_ao_beta_complex,1)) + deallocate(mo_coef_beta_tmp) + !call zgemm('N','C',ao_num,ao_num,elec_beta_num,(1.d0,0.d0), & + ! mo_coef_complex, size(mo_coef_complex,1), & + ! mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & + ! scf_density_matrix_ao_beta_complex, size(scf_density_matrix_ao_beta_complex,1)) + END_PROVIDER BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_complex, (ao_num,ao_num) ]