LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

working on scf kpts

Browse Source
This commit is contained in:
Kevin Gasperich 2020-03-18 15:55:53 -05:00
parent 84531d8021
commit 380cbdcbb5
7 changed files with 844 additions and 315 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

249
src/bitmask/bitmasks.irp.f
View File

@ -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

542
src/bitmask/core_inact_act_virt.irp.f
View File

@ -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
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)
! 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_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)
! 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
!
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

2
src/mo_basis/mos_cplx.irp.f
View File

@ -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))

78
src/mo_one_e_ints/ao_to_mo_cplx.irp.f
View File

@ -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

4
src/scf_utils/diagonalize_fock_cplx.irp.f
View File

@ -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

217
src/scf_utils/fock_matrix_cplx.irp.f
View File

@ -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. &
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
do i=1,ao_num
Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(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
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))
!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)*(<ij|kl>)
@ -594,13 +608,24 @@ END_PROVIDER
!G_a(i,l) -= D_a (k,j)*(<ij|kl>)
!G_b(i,l) -= D_b (k,j)*(<ij|kl>)
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)*(<ij|kl>)
@ -676,13 +729,24 @@ END_PROVIDER
!G_a(i,l) -= D_a (k,j)*(<ij|kl>)
!G_b(i,l) -= D_b (k,j)*(<ij|kl>)
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

51
src/scf_utils/scf_density_matrix_ao_cplx.irp.f
View File

@ -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