From 380cbdcbb5a997c7f884748dc56b6579bb13e631 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 18 Mar 2020 15:55:53 -0500 Subject: [PATCH] working on scf kpts --- src/bitmask/bitmasks.irp.f | 249 ++++++++ src/bitmask/core_inact_act_virt.irp.f | 546 ++++++++++-------- src/mo_basis/mos_cplx.irp.f | 2 +- src/mo_one_e_ints/ao_to_mo_cplx.irp.f | 78 +++ src/scf_utils/diagonalize_fock_cplx.irp.f | 4 +- src/scf_utils/fock_matrix_cplx.irp.f | 229 +++++--- .../scf_density_matrix_ao_cplx.irp.f | 51 ++ 7 files changed, 844 insertions(+), 315 deletions(-) diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index 03127b1c..75421967 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -254,3 +254,252 @@ BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] closed_shell_ref_bitmask(i,2) = ior(ref_bitmask(i,2),act_bitmask(i,2)) enddo END_PROVIDER + +!============================================! +! ! +! kpts ! +! ! +!============================================! + +!BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int) ] +! implicit none +! BEGIN_DOC +! ! Bitmask to include all possible MOs +! END_DOC +! +! integer :: i,j,k +! k=0 +! do j=1,N_int +! full_ijkl_bitmask(j) = 0_bit_kind +! do i=0,bit_kind_size-1 +! k=k+1 +! if (mo_class(k) /= 'Deleted') then +! full_ijkl_bitmask(j) = ibset(full_ijkl_bitmask(j),i) +! endif +! if (k == mo_num) exit +! enddo +! enddo +!END_PROVIDER +! +!BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask_4, (N_int,4) ] +! implicit none +! integer :: i +! do i=1,N_int +! full_ijkl_bitmask_4(i,1) = full_ijkl_bitmask(i) +! full_ijkl_bitmask_4(i,2) = full_ijkl_bitmask(i) +! full_ijkl_bitmask_4(i,3) = full_ijkl_bitmask(i) +! full_ijkl_bitmask_4(i,4) = full_ijkl_bitmask(i) +! enddo +!END_PROVIDER +! +!BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ] +! implicit none +! integer :: i +! do i=1,N_int +! core_inact_act_bitmask_4(i,1) = reunion_of_core_inact_act_bitmask(i,1) +! core_inact_act_bitmask_4(i,2) = reunion_of_core_inact_act_bitmask(i,1) +! core_inact_act_bitmask_4(i,3) = reunion_of_core_inact_act_bitmask(i,1) +! core_inact_act_bitmask_4(i,4) = reunion_of_core_inact_act_bitmask(i,1) +! enddo +!END_PROVIDER +! +!BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask_4, (N_int,4) ] +! implicit none +! integer :: i +! do i=1,N_int +! virt_bitmask_4(i,1) = virt_bitmask(i,1) +! virt_bitmask_4(i,2) = virt_bitmask(i,1) +! virt_bitmask_4(i,3) = virt_bitmask(i,1) +! virt_bitmask_4(i,4) = virt_bitmask(i,1) +! enddo +!END_PROVIDER +! +! +! +! +BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Hartree Fock bit mask + END_DOC + integer :: i,k + + hf_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + hf_bitmask_kpts(i,1,k) = iand(hf_bitmask(i,1),kpts_bitmask(i,k)) + hf_bitmask_kpts(i,2,k) = iand(hf_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), ref_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask + END_DOC + ref_bitmask_kpts = HF_bitmask_kpts +END_PROVIDER + + + +!BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6) ] +! implicit none +! BEGIN_DOC +! ! Bitmasks for generator determinants. +! ! (N_int, alpha/beta, hole/particle, generator). +! ! +! ! 3rd index is : +! ! +! ! * 1 : hole for single exc +! ! +! ! * 2 : particle for single exc +! ! +! ! * 3 : hole for 1st exc of double +! ! +! ! * 4 : particle for 1st exc of double +! ! +! ! * 5 : hole for 2nd exc of double +! ! +! ! * 6 : particle for 2nd exc of double +! ! +! END_DOC +! logical :: exists +! PROVIDE ezfio_filename full_ijkl_bitmask +! +! integer :: ispin, i +! do ispin=1,2 +! do i=1,N_int +! generators_bitmask(i,ispin,s_hole ) = reunion_of_inact_act_bitmask(i,ispin) +! generators_bitmask(i,ispin,s_part ) = reunion_of_act_virt_bitmask(i,ispin) +! generators_bitmask(i,ispin,d_hole1) = reunion_of_inact_act_bitmask(i,ispin) +! generators_bitmask(i,ispin,d_part1) = reunion_of_act_virt_bitmask(i,ispin) +! generators_bitmask(i,ispin,d_hole2) = reunion_of_inact_act_bitmask(i,ispin) +! generators_bitmask(i,ispin,d_part2) = reunion_of_act_virt_bitmask(i,ispin) +! enddo +! enddo +! +!END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the core and inactive and virtual bitmasks + END_DOC + integer :: i,k + do k=1,kpt_num + do i = 1, N_int + reunion_of_core_inact_bitmask_kpts(i,1,k) = ior(core_bitmask_kpts(i,1,k),inact_bitmask_kpts(i,1,k)) + reunion_of_core_inact_bitmask_kpts(i,2,k) = ior(core_bitmask_kpts(i,2,k),inact_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [integer(bit_kind), reunion_of_inact_act_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the inactive and active bitmasks + END_DOC + integer :: i,k + + do k=1,kpt_num + do i = 1, N_int + reunion_of_inact_act_bitmask_kpts(i,1,k) = ior(inact_bitmask_kpts(i,1,k),act_bitmask_kpts(i,1,k)) + reunion_of_inact_act_bitmask_kpts(i,2,k) = ior(inact_bitmask_kpts(i,2,k),act_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer(bit_kind), reunion_of_act_virt_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the inactive and active bitmasks + END_DOC + integer :: i,k + + do k=1,kpt_num + do i = 1, N_int + reunion_of_act_virt_bitmask_kpts(i,1,k) = ior(virt_bitmask_kpts(i,1,k),act_bitmask_kpts(i,1,k)) + reunion_of_act_virt_bitmask_kpts(i,2,k) = ior(virt_bitmask_kpts(i,2,k),act_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the core, inactive and active bitmasks + END_DOC + integer :: i,k + + do k=1,kpt_num + do i = 1, N_int + reunion_of_core_inact_act_bitmask_kpts(i,1,k) = ior(reunion_of_core_inact_bitmask_kpts(i,1,k),act_bitmask_kpts(i,1,k)) + reunion_of_core_inact_act_bitmask_kpts(i,2,k) = ior(reunion_of_core_inact_bitmask_kpts(i,2,k),act_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the inactive, active and virtual bitmasks + END_DOC + integer :: i,k + do k=1,kpt_num + do i = 1, N_int + reunion_of_bitmask_kpts(i,1,k) = ior(ior(act_bitmask_kpts(i,1,k),inact_bitmask_kpts(i,1,k)),virt_bitmask_kpts(i,1,k)) + reunion_of_bitmask_kpts(i,2,k) = ior(ior(act_bitmask_kpts(i,2,k),inact_bitmask_kpts(i,2,k)),virt_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), inact_virt_bitmask_kpts, (N_int,2,kpt_num)] +&BEGIN_PROVIDER [ integer(bit_kind), core_inact_virt_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + BEGIN_DOC + ! Reunion of the inactive and virtual bitmasks + END_DOC + integer :: i,k + do k=1,kpt_num + do i = 1, N_int + inact_virt_bitmask_kpts(i,1,k) = ior(inact_bitmask_kpts(i,1,k),virt_bitmask_kpts(i,1,k)) + inact_virt_bitmask_kpts(i,2,k) = ior(inact_bitmask_kpts(i,2,k),virt_bitmask_kpts(i,2,k)) + core_inact_virt_bitmask_kpts(i,1,k) = ior(core_bitmask_kpts(i,1,k),inact_virt_bitmask_kpts(i,1,k)) + core_inact_virt_bitmask_kpts(i,2,k) = ior(core_bitmask_kpts(i,2,k),inact_virt_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer(bit_kind), unpaired_alpha_electrons_kpts, (N_int,kpt_num)] + implicit none + BEGIN_DOC + ! Bitmask reprenting the unpaired alpha electrons in the HF_bitmask + END_DOC + integer :: i,k + unpaired_alpha_electrons_kpts = 0_bit_kind + do k = 1, kpt_num + do i = 1, N_int + unpaired_alpha_electrons_kpts(i,k) = xor(HF_bitmask_kpts(i,1,k),HF_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + integer :: i,k + + closed_shell_ref_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i = 1, N_int + closed_shell_ref_bitmask_kpts(i,1,k) = ior(ref_bitmask_kpts(i,1,k),act_bitmask_kpts(i,1,k)) + closed_shell_ref_bitmask_kpts(i,2,k) = ior(ref_bitmask_kpts(i,2,k),act_bitmask_kpts(i,2,k)) + enddo + enddo +END_PROVIDER + diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index c0484057..337e275b 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -418,6 +418,23 @@ END_PROVIDER ! kpts ! ! ! !============================================! +BEGIN_PROVIDER [ integer(bit_kind), kpts_bitmask , (N_int,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying each kpt + END_DOC + integer :: k,i,di + integer :: tmp_mo_list(mo_num_per_kpt) + kpts_bitmask = 0_bit_kind + do k=1,kpt_num + di=(k-1)*mo_num_per_kpt + do i=1,mo_num_per_kpt + tmp_mo_list(i) = i+di + enddo + call list_to_bitstring( kpts_bitmask(1,k), tmp_mo_list, mo_num_per_kpt, N_int) + enddo +END_PROVIDER + BEGIN_PROVIDER [ integer, n_core_orb_kpts, (kpt_num)] implicit none BEGIN_DOC @@ -524,7 +541,7 @@ BEGIN_PROVIDER [ integer, n_del_orb_kpts, (kpt_num)] END_PROVIDER BEGIN_PROVIDER [ integer, n_core_inact_orb_kpts, (kpt_num) ] - !todo: finish implementation for kpts (will need kpt_mask) + !todo: finish implementation for kpts (will need kpts_bitmask) implicit none BEGIN_DOC ! n_core + n_inact @@ -533,7 +550,7 @@ BEGIN_PROVIDER [ integer, n_core_inact_orb_kpts, (kpt_num) ] 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))) + n_core_inact_orb_kpts(k) += popcnt(iand(kpts_bitmask(i,k),reunion_of_core_inact_bitmask(i,1))) enddo enddo END_PROVIDER @@ -603,6 +620,24 @@ BEGIN_PROVIDER [integer, dim_list_del_orb_kpts] dim_list_del_orb_kpts = max(maxval(n_del_orb_kpts),1) END_PROVIDER +BEGIN_PROVIDER [integer, dim_list_core_inact_act_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_core_inact_act. + ! it is at least 1 + END_DOC + dim_list_core_inact_act_orb_kpts = max(maxval(n_core_inact_act_orb_kpts),1) +END_PROVIDER + +BEGIN_PROVIDER [integer, dim_list_inact_act_orb_kpts] + implicit none + BEGIN_DOC + ! dimensions for the allocation of list_inact_act. + ! it is at least 1 + END_DOC + dim_list_inact_act_orb_kpts = max(maxval(n_inact_act_orb_kpts),1) +END_PROVIDER + BEGIN_PROVIDER [integer, n_core_inact_act_orb_kpts, (kpt_num) ] implicit none BEGIN_DOC @@ -615,242 +650,273 @@ BEGIN_PROVIDER [integer, n_core_inact_act_orb_kpts, (kpt_num) ] END_PROVIDER + + +BEGIN_PROVIDER [ integer(bit_kind), core_bitmask_kpts , (N_int,2,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the core MOs + END_DOC + integer :: k,i + core_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + core_bitmask_kpts(i,1,k) = iand(core_bitmask(i,1),kpts_bitmask(i,k)) + core_bitmask_kpts(i,2,k) = iand(core_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), inact_bitmask_kpts , (N_int,2,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the inactive MOs + END_DOC + integer :: k,i + inact_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + inact_bitmask_kpts(i,1,k) = iand(inact_bitmask(i,1),kpts_bitmask(i,k)) + inact_bitmask_kpts(i,2,k) = iand(inact_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), act_bitmask_kpts , (N_int,2,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the active MOs + END_DOC + integer :: k,i + act_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + act_bitmask_kpts(i,1,k) = iand(act_bitmask(i,1),kpts_bitmask(i,k)) + act_bitmask_kpts(i,2,k) = iand(act_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask_kpts , (N_int,2,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the virtual MOs + END_DOC + integer :: k,i + virt_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + virt_bitmask_kpts(i,1,k) = iand(virt_bitmask(i,1),kpts_bitmask(i,k)) + virt_bitmask_kpts(i,2,k) = iand(virt_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), del_bitmask_kpts , (N_int,2,kpt_num) ] + implicit none + BEGIN_DOC + ! Bitmask identifying the deleted MOs + END_DOC + integer :: k,i + del_bitmask_kpts = 0_bit_kind + do k=1,kpt_num + do i=1,N_int + del_bitmask_kpts(i,1,k) = iand(del_bitmask(i,1),kpts_bitmask(i,k)) + del_bitmask_kpts(i,2,k) = iand(del_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_core_kpts , (dim_list_core_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_core_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are in the core. + END_DOC + integer :: i, n,k,di + list_core_kpts = 0 + list_core_kpts_reverse = 0 + + do k=1,kpt_num + n=0 + di = (k-1)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+di) == 'Core')then + n += 1 + list_core_kpts(n,k) = i + list_core_kpts_reverse(i,k) = n + endif + enddo + print *, 'Core MOs: ',k + print *, list_core_kpts(1:n_core_orb_kpts(k),k) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_inact_kpts , (dim_list_inact_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_inact_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are inactive. + END_DOC + integer :: i, n,k,di + list_inact_kpts = 0 + list_inact_kpts_reverse = 0 + + do k=1,kpt_num + n=0 + di = (k-1)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+di) == 'Inactive')then + n += 1 + list_inact_kpts(n,k) = i + list_inact_kpts_reverse(i,k) = n + endif + enddo + print *, 'Inactive MOs: ',k + print *, list_inact_kpts(1:n_inact_orb_kpts(k),k) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_virt_kpts , (dim_list_virt_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_virt_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are virtual. + END_DOC + integer :: i, n,k,di + list_virt_kpts = 0 + list_virt_kpts_reverse = 0 + + do k=1,kpt_num + n=0 + di = (k-1)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+di) == 'Virtual')then + n += 1 + list_virt_kpts(n,k) = i + list_virt_kpts_reverse(i,k) = n + endif + enddo + print *, 'Virtual MOs: ',k + print *, list_virt_kpts(1:n_virt_orb_kpts(k),k) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_del_kpts , (dim_list_del_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_del_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are deleted. + END_DOC + integer :: i, n,k,di + list_del_kpts = 0 + list_del_kpts_reverse = 0 + + do k=1,kpt_num + n=0 + di = (k-1)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+di) == 'Deleted')then + n += 1 + list_del_kpts(n,k) = i + list_del_kpts_reverse(i,k) = n + endif + enddo + print *, 'Deleted MOs: ',k + print *, list_del_kpts(1:n_del_orb_kpts(k),k) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ integer, list_act_kpts , (dim_list_act_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_act_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of MO indices which are active. + END_DOC + integer :: i, n,k,di + list_act_kpts = 0 + list_act_kpts_reverse = 0 + + do k=1,kpt_num + n=0 + di = (k-1)*mo_num_per_kpt + do i = 1, mo_num_per_kpt + if(mo_class(i+di) == 'Active')then + n += 1 + list_act_kpts(n,k) = i + list_act_kpts_reverse(i,k) = n + endif + enddo + print *, 'Active MOs: ',k + print *, list_act_kpts(1:n_act_orb_kpts(k),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 -! + + BEGIN_PROVIDER [ integer, list_core_inact_kpts , (dim_list_core_inact_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_core_inact_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of indices of the core and inactive MOs + END_DOC + integer :: i,itmp,k + list_core_inact_kpts_reverse = 0 + do k=1,kpt_num + !call bitstring_to_list(reunion_of_core_inact_bitmask(1,1), list_core_inact, itmp, N_int) + call bitstring_to_list(reunion_of_core_inact_bitmask_kpts(1,1,k), list_core_inact_kpts(1,k), itmp, N_int) + ASSERT (itmp == n_core_inact_orb_kpts(k)) + do i = 1, n_core_inact_orb_kpts(k) + list_core_inact_kpts_reverse(list_core_inact_kpts(i,k),k) = i + enddo + print *, 'Core and Inactive MOs: ',k + print *, list_core_inact_kpts(1:n_core_inact_orb_kpts(k),k) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer, list_core_inact_act_kpts , (dim_list_core_inact_act_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_core_inact_act_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of indices of the core inactive and active MOs + END_DOC + integer :: i,itmp,k + list_core_inact_act_kpts_reverse = 0 + do k=1,kpt_num + !call bitstring_to_list(reunion_of_core_inact_act_bitmask(1,1), list_core_inact_act, itmp, N_int) + call bitstring_to_list(reunion_of_core_inact_act_bitmask_kpts(1,1,k), list_core_inact_act_kpts(1,k), itmp, N_int) + ASSERT (itmp == n_core_inact_act_orb_kpts(k)) + do i = 1, n_core_inact_act_orb_kpts(k) + list_core_inact_act_kpts_reverse(list_core_inact_act_kpts(i,k),k) = i + enddo + print *, 'Core, Inactive and Active MOs: ',k + print *, list_core_inact_act_kpts(1:n_core_inact_act_orb_kpts(k),k) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer, list_inact_act_kpts , (dim_list_inact_act_orb_kpts,kpt_num) ] +&BEGIN_PROVIDER [ integer, list_inact_act_kpts_reverse, (mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! List of indices of the inactive and active MOs + END_DOC + integer :: i,itmp,k + list_inact_act_kpts_reverse = 0 + do k=1,kpt_num + call bitstring_to_list(reunion_of_inact_act_bitmask_kpts(1,1,k), list_inact_act_kpts(1,k), itmp, N_int) + ASSERT (itmp == n_inact_act_orb_kpts(k)) + do i = 1, n_inact_act_orb_kpts(k) + list_inact_act_kpts_reverse(list_inact_act_kpts(i,k),k) = i + enddo + print *, 'Inactive and Active MOs: ',k + print *, list_inact_act_kpts(1:n_inact_act_orb_kpts(k),k) + enddo +END_PROVIDER + diff --git a/src/mo_basis/mos_cplx.irp.f b/src/mo_basis/mos_cplx.irp.f index 6b7e14c7..e25e7717 100644 --- a/src/mo_basis/mos_cplx.irp.f +++ b/src/mo_basis/mos_cplx.irp.f @@ -317,7 +317,7 @@ subroutine ao_to_mo_kpts(A_ao,LDA_ao,A_mo,LDA_mo) 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, & + (1.d0,0.d0), A_ao(:,:,k),LDA_ao, & mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), & (0.d0,0.d0), T, size(T,1)) diff --git a/src/mo_one_e_ints/ao_to_mo_cplx.irp.f b/src/mo_one_e_ints/ao_to_mo_cplx.irp.f index 2530caf0..875d84a9 100644 --- a/src/mo_one_e_ints/ao_to_mo_cplx.irp.f +++ b/src/mo_one_e_ints/ao_to_mo_cplx.irp.f @@ -66,3 +66,81 @@ BEGIN_PROVIDER [ complex*16, S_mo_coef_complex, (ao_num, mo_num) ] END_PROVIDER +!============================================! +! ! +! kpts ! +! ! +!============================================! + +subroutine mo_to_ao_kpts(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis + ! + ! (S.C).A_mo.(S.C)t + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_mo(LDA_mo,mo_num_per_kpt,kpt_num) + complex*16, intent(out) :: A_ao(LDA_ao,ao_num_per_kpt,kpt_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(mo_num_per_kpt,ao_num_per_kpt) ) + integer :: k + do k=1,kpt_num + call zgemm('N','C', mo_num_per_kpt, ao_num_per_kpt, mo_num_per_kpt, & + (1.d0,0.d0), A_mo(:,:,k),size(A_mo,1), & + S_mo_coef_kpts(:,:,k), size(S_mo_coef_kpts,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num_per_kpt, ao_num_per_kpt, mo_num_per_kpt, & + (1.d0,0.d0), S_mo_coef_kpts(:,:,k), size(S_mo_coef_kpts,1), & + T, size(T,1), & + (0.d0,0.d0), A_ao(:,:,k), size(A_ao,1)) + enddo + deallocate(T) +end + +subroutine mo_to_ao_no_overlap_kpts(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the S^-1 AO basis + ! Useful for density matrix + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_mo(LDA_mo,mo_num_per_kpt,kpt_num) + complex*16, intent(out) :: A_ao(LDA_ao,ao_num_per_kpt,kpt_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(mo_num_per_kpt,ao_num_per_kpt) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + integer :: k + do k=1,kpt_num + call zgemm('N','C', mo_num_per_kpt, ao_num_per_kpt, mo_num_per_kpt, & + (1.d0,0.d0), A_mo(:,:,k),size(A_mo,1), & + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num_per_kpt, ao_num_per_kpt, mo_num_per_kpt, & + (1.d0,0.d0), mo_coef_kpts(:,:,k),size(mo_coef_kpts,1), & + T, size(T,1), & + (0.d0,0.d0), A_ao(:,:,k), size(A_ao,1)) + enddo + deallocate(T) +end + +BEGIN_PROVIDER [ complex*16, S_mo_coef_kpts, (ao_num_per_kpt, mo_num_per_kpt, kpt_num) ] + implicit none + BEGIN_DOC + ! Product S.C where S is the overlap matrix in the AO basis and C the mo_coef matrix. + 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_overlap_kpts(:,:,k), size(ao_overlap_kpts,1), & + mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), & + (0.d0,0.d0), & + S_mo_coef_kpts(:,:,k), size(S_mo_coef_kpts,1)) + enddo +END_PROVIDER + diff --git a/src/scf_utils/diagonalize_fock_cplx.irp.f b/src/scf_utils/diagonalize_fock_cplx.irp.f index de38c767..83d4b00f 100644 --- a/src/scf_utils/diagonalize_fock_cplx.irp.f +++ b/src/scf_utils/diagonalize_fock_cplx.irp.f @@ -94,11 +94,11 @@ BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts, (ao_num_per_kpt,m ! Insert level shift here !todo: elec per kpt - do i = elec_beta_num_per_kpt(k)+1, elec_alpha_num_per_kpt(k) + do i = elec_beta_num_kpts(k)+1, elec_alpha_num_kpts(k) F(i,i) += 0.5d0*level_shift enddo - do i = elec_alpha_num_per_kpt(k)+1, mo_num_per_kpt + do i = elec_alpha_num_kpts(k)+1, mo_num_per_kpt F(i,i) += level_shift enddo diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index 94508570..6b1fc808 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -475,54 +475,58 @@ END_PROVIDER -BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ] +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_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)) + call ao_to_mo_kpts(Fock_matrix_ao_alpha_kpts,size(Fock_matrix_ao_alpha_kpts,1), & + Fock_matrix_mo_alpha_kpts,size(Fock_matrix_mo_alpha_kpts,1)) END_PROVIDER -BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ] +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_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)) + call ao_to_mo_kpts(Fock_matrix_ao_beta_kpts,size(Fock_matrix_ao_beta_kpts,1), & + Fock_matrix_mo_beta_kpts,size(Fock_matrix_mo_beta_kpts,1)) END_PROVIDER -BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] +BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_kpts, (ao_num_per_kpt, ao_num_per_kpt,kpt_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)) + call mo_to_ao_kpts(Fock_matrix_mo_kpts,size(Fock_matrix_mo_kpts,1), & + Fock_matrix_ao_kpts,size(Fock_matrix_ao_kpts,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) + integer :: k + do k=1,kpt_num + if ( (elec_alpha_num_kpts(k) == elec_beta_num_kpts(k)).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num_per_kpt + do i=1,ao_num_per_kpt + Fock_matrix_ao_kpts(i,j,k) = Fock_matrix_ao_alpha_kpts(i,j,k) + enddo 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 + else + !call mo_to_ao_complex(Fock_matrix_mo_kpts,size(Fock_matrix_mo_kpts,1), & + call mo_to_ao_kpts(Fock_matrix_mo_kpts,size(Fock_matrix_mo_kpts,1), & + Fock_matrix_ao_kpts,size(Fock_matrix_ao_kpts,1)) + endif + enddo 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) ] + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] use map_module implicit none BEGIN_DOC @@ -534,11 +538,11 @@ END_PROVIDER 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(:,:) + 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) + ao_two_e_integral_alpha_kpts = (0.d0,0.d0) + ao_two_e_integral_beta_kpts = (0.d0,0.d0) PROVIDE ao_two_e_integrals_in_map integer(omp_lock_kind) :: lck(ao_num) @@ -549,19 +553,21 @@ END_PROVIDER 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 + integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l !$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 kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & !$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) + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) 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)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) @@ -587,6 +593,14 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) + kpt_i = (i-1)/kpt_num +1 + kpt_j = (j-1)/kpt_num +1 + kpt_k = (k-1)/kpt_num +1 + kpt_l = (l-1)/kpt_num +1 + idx_i = mod(i,kpt_num) + idx_j = mod(j,kpt_num) + idx_k = mod(k,kpt_num) + idx_l = mod(l,kpt_num) integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -594,13 +608,24 @@ END_PROVIDER !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 + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i) += c0 + endif - 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 + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif enddo else ! real part do k2=1,4 @@ -611,23 +636,42 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) + kpt_i = (i-1)/kpt_num +1 + kpt_j = (j-1)/kpt_num +1 + kpt_k = (k-1)/kpt_num +1 + kpt_l = (l-1)/kpt_num +1 + idx_i = mod(i,kpt_num) + idx_j = mod(j,kpt_num) + idx_k = mod(k,kpt_num) + idx_l = mod(l,kpt_num) integral = values(k1) - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i) += c0 + endif - 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 + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif 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 + ao_two_e_integral_alpha_kpts += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts += 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 @@ -636,15 +680,16 @@ END_PROVIDER !$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 kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & !$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) + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) 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)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) @@ -669,6 +714,14 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) + kpt_i = (i-1)/kpt_num +1 + kpt_j = (j-1)/kpt_num +1 + kpt_k = (k-1)/kpt_num +1 + kpt_l = (l-1)/kpt_num +1 + idx_i = mod(i,kpt_num) + idx_j = mod(j,kpt_num) + idx_k = mod(k,kpt_num) + idx_l = mod(l,kpt_num) integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -676,13 +729,24 @@ END_PROVIDER !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 + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i) += c0 + endif - 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 + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif enddo else ! real part do k2=1,4 @@ -693,23 +757,42 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) + kpt_i = (i-1)/kpt_num +1 + kpt_j = (j-1)/kpt_num +1 + kpt_k = (k-1)/kpt_num +1 + kpt_l = (l-1)/kpt_num +1 + idx_i = mod(i,kpt_num) + idx_j = mod(j,kpt_num) + idx_k = mod(k,kpt_num) + idx_l = mod(l,kpt_num) integral = values(k1) - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i) += c0 + endif - 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 + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif 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 + ao_two_e_integral_alpha_kpts += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts += 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 @@ -717,18 +800,20 @@ END_PROVIDER 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) ] + BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] +&BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_beta_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_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) + integer :: i,j,k + do k=1,kpt_num + do j=1,ao_num_per_kpt + do i=1,ao_num_per_kpt + Fock_matrix_ao_alpha_kpts(i,j,k) = ao_one_e_integrals_kpts(i,j,k) + ao_two_e_integral_alpha_kpts(i,j,k) + Fock_matrix_ao_beta_kpts (i,j,k) = ao_one_e_integrals_kpts(i,j,k) + ao_two_e_integral_beta_kpts (i,j,k) + enddo enddo 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 a6a66863..9726690c 100644 --- a/src/scf_utils/scf_density_matrix_ao_cplx.irp.f +++ b/src/scf_utils/scf_density_matrix_ao_cplx.irp.f @@ -73,3 +73,54 @@ BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_complex, (ao_num,ao_num) ] END_PROVIDER +!============================================! +! ! +! kpts ! +! ! +!============================================! + +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_alpha_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! $C.C^t$ over $\alpha$ MOs + END_DOC + + integer :: k + do k=1,kpt_num + call zgemm('N','C',ao_num_per_kpt,ao_num_per_kpt,elec_alpha_num_kpts(k),(1.d0,0.d0), & + mo_coef_kpts(1,1,k), size(mo_coef_kpts,1), & + mo_coef_kpts(1,1,k), size(mo_coef_kpts,1), (0.d0,0.d0), & + scf_density_matrix_ao_alpha_kpts(1,1,k), size(scf_density_matrix_ao_alpha_kpts,1)) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_beta_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! $C.C^t$ over $\beta$ MOs + END_DOC + + integer :: k + do k=1,kpt_num + call zgemm('N','C',ao_num_per_kpt,ao_num_per_kpt,elec_beta_num_kpts(k),(1.d0,0.d0), & + mo_coef_kpts(1,1,k), size(mo_coef_kpts,1), & + mo_coef_kpts(1,1,k), size(mo_coef_kpts,1), (0.d0,0.d0), & + scf_density_matrix_ao_beta_kpts(1,1,k), size(scf_density_matrix_ao_beta_kpts,1)) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, scf_density_matrix_ao_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! Sum of $\alpha$ and $\beta$ density matrices + END_DOC + ASSERT (size(scf_density_matrix_ao_kpts,1) == size(scf_density_matrix_ao_alpha_kpts,1)) + if (elec_alpha_num== elec_beta_num) then + scf_density_matrix_ao_kpts = scf_density_matrix_ao_alpha_kpts + scf_density_matrix_ao_alpha_kpts + else + ASSERT (size(scf_density_matrix_ao_kpts,1) == size(scf_density_matrix_ao_beta_kpts ,1)) + scf_density_matrix_ao_kpts = scf_density_matrix_ao_alpha_kpts + scf_density_matrix_ao_beta_kpts + endif + +END_PROVIDER +