LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-10-02 06:20:59 +02:00
Code Issues Releases Wiki Activity

Corrected bug in ROHF

Browse Source
This commit is contained in:
Anthony Scemama 2014-04-25 23:55:54 +02:00
parent 13d4bf8b0e
commit 6a9ebb343b
8 changed files with 183 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/BiInts/ao_bi_integrals.irp.f
View File

@ -160,20 +160,17 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
thresh = ao_integrals_threshold thresh = ao_integrals_threshold
integer :: n_centers, i integer :: n_centers, i
integer*1 :: center_count(nucl_num)
PROVIDE gauleg_t2 ao_nucl all_utils PROVIDE gauleg_t2 ao_nucl all_utils
if (ao_overlap_abs(j,l) < thresh) then if (ao_overlap_abs(j,l) < thresh) then
buffer_value = 0. buffer_value = 0._integral_kind
return return
endif endif
center_count = 0
do i = 1, ao_num do i = 1, ao_num
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
buffer_value(i) = 0. buffer_value(i) = 0._integral_kind
cycle cycle
endif endif
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
@ -211,7 +208,6 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
integer :: n_integrals, n_centers integer :: n_integrals, n_centers
integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
integer*1 :: center_count(nucl_num)
PROVIDE gauleg_t2 ao_integrals_map all_utils PROVIDE gauleg_t2 ao_integrals_map all_utils
integral = ao_bielec_integral(1,1,1,1) integral = ao_bielec_integral(1,1,1,1)
@ -243,7 +239,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
call cpu_time(cpu_1) call cpu_time(cpu_1)
!$OMP PARALLEL PRIVATE(i,j,k,l,kk, & !$OMP PARALLEL PRIVATE(i,j,k,l,kk, &
!$OMP integral,buffer_i,buffer_value,n_integrals, & !$OMP integral,buffer_i,buffer_value,n_integrals, &
!$OMP cpu_2,wall_2,i1,j1,center_count) & !$OMP cpu_2,wall_2,i1,j1) &
!$OMP DEFAULT(NONE) & !$OMP DEFAULT(NONE) &
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
@ -252,7 +248,6 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
allocate(buffer_i(size_buffer)) allocate(buffer_i(size_buffer))
allocate(buffer_value(size_buffer)) allocate(buffer_value(size_buffer))
n_integrals = 0 n_integrals = 0
center_count = 0
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do kk=1,lmax do kk=1,lmax

1
src/BiInts/bi_integrals.ezfio_config
View File

@ -5,4 +5,5 @@ bielec_integrals
write_mo_integrals logical write_mo_integrals logical
threshold_ao double precision threshold_ao double precision
threshold_mo double precision threshold_mo double precision
direct logical

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

@ -296,94 +296,131 @@ end
BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj, (mo_tot_num_align,mo_tot_num) ] 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_exchange, (mo_tot_num_align,mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (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 implicit none
BEGIN_DOC BEGIN_DOC
! Transform Bi-electronic integrals <ij|ij> and <ij|ji> ! Transform Bi-electronic integrals <ij|ij> and <ij|ji>
END_DOC END_DOC
integer :: i,j,p,q,r,s integer :: i,j,p,q,r,s
double precision :: integral, c double precision :: c
integer :: n, pp real(integral_kind) :: integral
double precision, allocatable :: int_value(:) integer :: n, pp
integer, allocatable :: int_idx(:) real(integral_kind), allocatable :: int_value(:)
integer, allocatable :: int_idx(:)
double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:)
PROVIDE ao_bielec_integrals_in_map PROVIDE ao_integrals_threshold
if (.not.do_direct_integrals) then
PROVIDE ao_bielec_integrals_in_map
endif
mo_bielec_integral_jj = 0.d0 mo_bielec_integral_jj = 0.d0
mo_bielec_integral_jj_exchange = 0.d0 mo_bielec_integral_jj_exchange = 0.d0
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & !$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, &
!$OMP iqrs, iqsr,iqri,iqis) & !$OMP iqrs, iqsr,iqri,iqis) &
!$OMP SHARED(mo_tot_num,mo_coef_transp,mo_tot_num_align,ao_num) & !$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) !$OMP ao_integrals_threshold,do_direct_integrals) &
!$OMP REDUCTION(+:mo_bielec_integral_jj,mo_bielec_integral_jj_exchange)
allocate( int_value(ao_num), int_idx(ao_num), & allocate( int_value(ao_num), int_idx(ao_num), &
iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num), & iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),&
iqsr(mo_tot_num_align,ao_num) ) iqsr(mo_tot_num_align,ao_num) )
!$OMP DO !$OMP DO SCHEDULE (guided)
do q=1,ao_num do s=1,ao_num
do s=1,ao_num do q=1,ao_num
do j=1,ao_num do j=1,ao_num
!DIR$ VECTOR ALIGNED !DIR$ VECTOR ALIGNED
do i=1,mo_tot_num do i=1,mo_tot_num
iqrs(i,j) = 0.d0 iqrs(i,j) = 0.d0
iqsr(i,j) = 0.d0 iqsr(i,j) = 0.d0
enddo enddo
enddo enddo
do r=1,ao_num if (do_direct_integrals) then
call get_ao_bielec_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) double precision :: ao_bielec_integral
do pp=1,n do r=1,ao_num
p = int_idx(pp) call compute_ao_bielec_integrals(q,r,s,ao_num,int_value)
integral = int_value(pp) do p=1,ao_num
!DIR$ VECTOR ALIGNED integral = int_value(p)
do i=1,mo_tot_num if (abs(integral) > ao_integrals_threshold) then
iqrs(i,r) += mo_coef_transp(i,p) * integral !DIR$ VECTOR ALIGNED
enddo do i=1,mo_tot_num
enddo iqrs(i,r) += mo_coef_transp(i,p) * integral
call get_ao_bielec_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) enddo
do pp=1,n endif
p = int_idx(pp) enddo
integral = int_value(pp) call compute_ao_bielec_integrals(q,s,r,ao_num,int_value)
!DIR$ VECTOR ALIGNED do p=1,ao_num
do i=1,mo_tot_num integral = int_value(p)
iqsr(i,r) += mo_coef_transp(i,p) * integral if (abs(integral) > ao_integrals_threshold) then
enddo !DIR$ VECTOR ALIGNED
enddo do i=1,mo_tot_num
enddo iqsr(i,r) += mo_coef_transp(i,p) * integral
iqis = 0.d0 enddo
iqri = 0.d0 endif
do r=1,ao_num enddo
!DIR$ VECTOR ALIGNED enddo
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 else
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 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)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqrs(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
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)
if (abs(integral) > ao_integrals_threshold) then
!DIR$ VECTOR ALIGNED
do i=1,mo_tot_num
iqsr(i,r) += mo_coef_transp(i,p) * integral
enddo
endif
enddo
enddo
endif
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 END_PROVIDER

18
src/BiInts/options.irp.f
View File

@ -107,3 +107,21 @@ BEGIN_PROVIDER [ double precision, mo_integrals_threshold ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ logical, do_direct_integrals ]
implicit none
BEGIN_DOC
! If True, compute integrals on the fly
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_bielec_integrals_direct(has)
if (has) then
call ezfio_get_bielec_integrals_direct(do_direct_integrals)
else
do_direct_integrals = .False.
call ezfio_set_bielec_integrals_direct(do_direct_integrals)
endif
END_PROVIDER

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

@ -92,24 +92,40 @@ END_PROVIDER
integer, allocatable :: ao_ints_idx(:) integer, allocatable :: ao_ints_idx(:)
double precision :: integral double precision :: integral
double precision :: ao_bielec_integral double precision :: ao_bielec_integral
!DIR$ ATTRIBUTES ALIGN : 32 :: ao_ints_idx, ao_ints_val !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_ints_idx, ao_ints_val
if (do_direct_SCF) then if (do_direct_integrals) then
PROVIDE all_utils
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral) &
!$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)
!$OMP DO SCHEDULE(GUIDED)
do j=1,ao_num do j=1,ao_num
do i=1,ao_num do i=1,j
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j) Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_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 l=1,ao_num
do k=1,ao_num do k=1,ao_num
if ((abs(HF_density_matrix_ao_alpha(k,l)) > 1.d-9).or. & 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 (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) integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta (k,l)) * ao_bielec_integral(k,l,i,j)
Fock_matrix_alpha_ao(i,j) =Fock_matrix_alpha_ao(i,j) +( HF_density_matrix_ao_alpha(k,l) * integral) Fock_matrix_alpha_ao(i,j) += integral
Fock_matrix_beta_ao (i,j) =Fock_matrix_beta_ao (i,j) +( HF_density_matrix_ao_beta (k,l) * integral) Fock_matrix_beta_ao (i,j) += integral
integral = ao_bielec_integral(k,j,i,l)
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
endif endif
enddo enddo
enddo enddo
Fock_matrix_alpha_ao(j,i) = Fock_matrix_alpha_ao(i,j)
Fock_matrix_beta_ao (j,i) = Fock_matrix_beta_ao (i,j)
enddo enddo
enddo enddo
!$OMP END DO NOWAIT
!$OMP END PARALLEL
else else
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) & !$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) &
@ -117,23 +133,25 @@ END_PROVIDER
!$OMP ao_num_align,Fock_matrix_beta_ao,HF_density_matrix_ao_alpha, & !$OMP ao_num_align,Fock_matrix_beta_ao,HF_density_matrix_ao_alpha, &
!$OMP HF_density_matrix_ao_beta) !$OMP HF_density_matrix_ao_beta)
allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align)) allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align))
!$OMP DO !$OMP DO SCHEDULE(GUIDED)
do i=1,ao_num do j=1,ao_num
do j=1,ao_num !DIR$ VECTOR ALIGNED
do i=1,ao_num
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j) Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j)
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j) Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j)
enddo enddo
do j=1,ao_num do l=1,ao_num
do l=1,ao_num do i=1,ao_num
call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax) call get_ao_bielec_integrals_non_zero(i,l,j,ao_num,ao_ints_val,ao_ints_idx,kmax)
!DIR$ VECTOR ALIGNED
do k1=1,kmax do k1=1,kmax
k = ao_ints_idx(k1) k = ao_ints_idx(k1)
integral = ao_ints_val(k1)+ao_ints_val(k1) integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * ao_ints_val(k1)
Fock_matrix_alpha_ao(i,j) += HF_density_matrix_ao_alpha(k,l) * integral Fock_matrix_alpha_ao(i,j) += integral
Fock_matrix_beta_ao (i,j) += HF_density_matrix_ao_beta (k,l) * integral Fock_matrix_beta_ao (i,j) += integral
integral = -ao_ints_val(k1) integral = ao_ints_val(k1)
Fock_matrix_alpha_ao(i,l) += HF_density_matrix_ao_alpha(k,j) * integral Fock_matrix_alpha_ao(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
Fock_matrix_beta_ao (i,l) += HF_density_matrix_ao_beta (k,j) * integral Fock_matrix_beta_ao (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
enddo enddo
enddo enddo
enddo enddo

5
src/Hartree_Fock/README.rst
View File

@ -1,8 +1,3 @@
===================
Hartree-Fock Module
===================
Needed Modules Needed Modules
============== ==============

1
src/Hartree_Fock/hartree_fock.ezfio_config
View File

@ -1,6 +1,5 @@
hartree_fock hartree_fock
thresh_scf double precision thresh_scf double precision
n_it_scf_max integer n_it_scf_max integer
direct logical
diis logical diis logical

18
src/Hartree_Fock/options.irp.f
View File

@ -35,24 +35,6 @@ BEGIN_PROVIDER [ integer, n_it_scf_max]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ logical, do_direct_SCF ]
implicit none
BEGIN_DOC
! If True, compute integrals on the fly
END_DOC
logical :: has
PROVIDE ezfio_filename
call ezfio_has_Hartree_Fock_direct(has)
if (has) then
call ezfio_get_Hartree_Fock_direct(do_direct_SCF)
else
do_direct_SCF = .False.
call ezfio_set_Hartree_Fock_direct(do_direct_SCF)
endif
END_PROVIDER
BEGIN_PROVIDER [ logical, do_DIIS ] BEGIN_PROVIDER [ logical, do_DIIS ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC