LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-03-13 20:32:26 +01:00
Code Issues Releases Wiki Activity

Hartree-Fock without MO bielec map

Browse Source
This commit is contained in:
Anthony Scemama 2014-04-25 22:18:26 +02:00
parent b9cfe96b17
commit 7db8b1eab0
3 changed files with 195 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

156
src/BiInts/mo_bi_integrals.irp.f
View File

@ -72,7 +72,6 @@ subroutine add_integrals_to_map(mask_ijkl)
PROVIDE N_int ao_bielec_integrals_in_map ao_integrals_map mo_coef mo_coef_transp
!Get list of MOs for i,j,k and l
!-------------------------------
@ -82,9 +81,6 @@ subroutine add_integrals_to_map(mask_ijkl)
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
l1_global = 0
size_buffer = min(ao_num*ao_num*ao_num,16000000)
write(output_BiInts,*) 'Providing the molecular integrals '
@ -294,63 +290,101 @@ subroutine add_integrals_to_map(mask_ijkl)
end
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j= 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i= 1, mo_tot_num
mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num
mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
enddo
! Padding
do i= mo_tot_num+1, mo_tot_num_align
mo_bielec_integral_jj_exchange(i,j) = 0.d0
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Bi-electronic integrals <ij|ij> - <ij|ji>
END_DOC
integer :: i,j
PROVIDE mo_bielec_integrals_in_map
do j = 1, mo_tot_num
!DIR$ VECTOR ALIGNED
do i = 1,mo_tot_num_align
mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Transform Bi-electronic integrals <ij|ij> and <ij|ji>
END_DOC
integer :: i,j,p,q,r,s
double precision :: integral, c
integer :: n, pp
double precision, allocatable :: int_value(:)
integer, allocatable :: int_idx(:)
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
PROVIDE ao_bielec_integrals_in_map
mo_bielec_integral_jj = 0.d0
mo_bielec_integral_jj_exchange = 0.d0
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, &
!$OMP iqrs, iqsr,iqri,iqis) &
!$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num) &
!$OMP REDUCTION(+:mo_bielec_integral_jj,mo_bielec_integral_jj_exchange)
allocate( int_value(ao_num), int_idx(ao_num), &
iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num), &
iqsr(mo_tot_num_align,ao_num) )
!$OMP DO
do q=1,ao_num
do s=1,ao_num
do j=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,j) = 0.d0
iqsr(i,j) = 0.d0
enddo
enddo
do r=1,ao_num
call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
enddo
call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n)
do pp=1,n
p = int_idx(pp)
integral = int_value(pp)
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
enddo
enddo
iqis = 0.d0
iqri = 0.d0
do r=1,ao_num
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqis(i) += mo_coef_transp(i,r) * iqrs(i,r)
iqri(i) += mo_coef_transp(i,r) * iqsr(i,r)
enddo
enddo
do i=1,mo_tot_num
!DIR$ VECTOR ALIGNED
do j=1,mo_tot_num
c = mo_coef_transp(j,q)*mo_coef_transp(j,s)
mo_bielec_integral_jj(j,i) += c * iqis(i)
mo_bielec_integral_jj_exchange(j,i) += c * iqri(i)
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
deallocate(iqrs,iqsr,int_value,int_idx)
!$OMP END PARALLEL
mo_bielec_integral_jj_anti = mo_bielec_integral_jj - mo_bielec_integral_jj_exchange
END_PROVIDER

131
src/Hartree_Fock/Fock_matrix_mo.irp.f
View File

@ -79,48 +79,115 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num,mo_tot_num) ]
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ]
&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao, (ao_num_align, ao_num) ]
implicit none
BEGIN_DOC
! Alpha Fock matrix in AO basis set
END_DOC
integer :: i,j,k,l,k1,kmax
double precision, allocatable :: ao_ints_val(:)
integer, allocatable :: ao_ints_idx(:)
double precision :: integral
double precision :: ao_bielec_integral
!DIR$ ATTRIBUTES ALIGN : 32 :: ao_ints_idx, ao_ints_val
if (do_direct_SCF) then
do j=1,ao_num
do i=1,ao_num
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j)
do l=1,ao_num
do k=1,ao_num
if ((abs(HF_density_matrix_ao_alpha(k,l)) > 1.d-9).or. &
(abs(HF_density_matrix_ao_beta (k,l)) > 1.d-9)) then
integral = 2.d0*ao_bielec_integral(k,l,i,j)-ao_bielec_integral(k,j,i,l)
Fock_matrix_alpha_ao(i,j) =Fock_matrix_alpha_ao(i,j) +( HF_density_matrix_ao_alpha(k,l) * integral)
Fock_matrix_beta_ao (i,j) =Fock_matrix_beta_ao (i,j) +( HF_density_matrix_ao_beta (k,l) * integral)
endif
enddo
enddo
enddo
enddo
else
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) &
!$OMP SHARED(ao_num,Fock_matrix_alpha_ao,ao_mono_elec_integral,&
!$OMP ao_num_align,Fock_matrix_beta_ao,HF_density_matrix_ao_alpha, &
!$OMP HF_density_matrix_ao_beta)
allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align))
!$OMP DO
do j=1,ao_num
do i=1,ao_num
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j)
do l=1,ao_num
call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax)
do k1=1,kmax
k = ao_ints_idx(k1)
integral = ao_ints_val(k1)+ao_ints_val(k1)
Fock_matrix_alpha_ao(i,j) += HF_density_matrix_ao_alpha(k,l) * integral
Fock_matrix_beta_ao (i,j) += HF_density_matrix_ao_beta (k,l) * integral
enddo
call get_ao_bielec_integrals_non_zero(i,j,l,ao_num,ao_ints_val,ao_ints_idx,kmax)
do k1=1,kmax
k = ao_ints_idx(k1)
integral = -ao_ints_val(k1)
Fock_matrix_alpha_ao(i,j) += HF_density_matrix_ao_alpha(k,l) * integral
Fock_matrix_beta_ao (i,j) += HF_density_matrix_ao_beta (k,l) * integral
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(ao_ints_val,ao_ints_idx)
!$OMP END PARALLEL
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
integer :: i,j,n
double precision :: get_mo_bielec_integral
do j=1,mo_tot_num
do i=1,mo_tot_num
Fock_matrix_alpha_mo(i,j) = mo_mono_elec_integral(i,j)
do n=1,elec_beta_num
Fock_matrix_alpha_mo(i,j) += 2.d0*get_mo_bielec_integral(i,n,j,n,mo_integrals_map) -&
get_mo_bielec_integral(i,n,n,j,mo_integrals_map)
enddo
do n=elec_beta_num+1,elec_alpha_num
Fock_matrix_alpha_mo(i,j) += get_mo_bielec_integral(i,n,j,n,mo_integrals_map) -&
get_mo_bielec_integral(i,n,n,j,mo_integrals_map)
enddo
enddo
enddo
double precision, allocatable :: T(:,:)
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_alpha_ao,size(Fock_matrix_alpha_ao,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, ao_num_align)
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, Fock_matrix_alpha_mo, mo_tot_num_align)
deallocate(T)
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num,mo_tot_num) ]
BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
integer :: i,j,n
double precision :: get_mo_bielec_integral
do j=1,mo_tot_num
do i=1,mo_tot_num
Fock_matrix_beta_mo(i,j) = mo_mono_elec_integral(i,j)
do n=1,elec_beta_num
Fock_matrix_beta_mo(i,j) += 2.d0*get_mo_bielec_integral(i,n,j,n,mo_integrals_map) -&
get_mo_bielec_integral(i,n,n,j,mo_integrals_map)
enddo
do n=elec_beta_num+1,elec_alpha_num
Fock_matrix_beta_mo(i,j) += get_mo_bielec_integral(i,n,j,n,mo_integrals_map)
enddo
enddo
enddo
double precision, allocatable :: T(:,:)
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_beta_ao,size(Fock_matrix_beta_ao,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, ao_num_align)
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, Fock_matrix_beta_mo, mo_tot_num_align)
deallocate(T)
END_PROVIDER

2
src/Hartree_Fock/mo_SCF_iterations.irp.f
View File

@ -28,6 +28,6 @@ program scf_iteration
endif
mo_label = "Canonical"
TOUCH mo_label mo_coef
call save_mos
! call save_mos
end