From 6a9ebb343ba40cfafd998a05c298ef992959047d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 25 Apr 2014 23:55:54 +0200 Subject: [PATCH] Corrected bug in ROHF --- src/BiInts/ao_bi_integrals.irp.f | 11 +- src/BiInts/bi_integrals.ezfio_config | 1 + src/BiInts/mo_bi_integrals.irp.f | 215 ++++++++++++--------- src/BiInts/options.irp.f | 18 ++ src/Hartree_Fock/Fock_matrix_mo.irp.f | 52 +++-- src/Hartree_Fock/README.rst | 5 - src/Hartree_Fock/hartree_fock.ezfio_config | 1 - src/Hartree_Fock/options.irp.f | 18 -- 8 files changed, 183 insertions(+), 138 deletions(-) diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index fbd6fb98..c9beb9c3 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -160,20 +160,17 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) thresh = ao_integrals_threshold integer :: n_centers, i - integer*1 :: center_count(nucl_num) PROVIDE gauleg_t2 ao_nucl all_utils if (ao_overlap_abs(j,l) < thresh) then - buffer_value = 0. + buffer_value = 0._integral_kind return endif - center_count = 0 - do i = 1, ao_num if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then - buffer_value(i) = 0. + buffer_value(i) = 0._integral_kind cycle endif !DIR$ FORCEINLINE @@ -211,7 +208,6 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer :: n_integrals, n_centers 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 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) !$OMP PARALLEL PRIVATE(i,j,k,l,kk, & !$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 SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & !$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_value(size_buffer)) n_integrals = 0 - center_count = 0 !$OMP DO SCHEDULE(dynamic) do kk=1,lmax diff --git a/src/BiInts/bi_integrals.ezfio_config b/src/BiInts/bi_integrals.ezfio_config index c0a3989a..e780dce8 100644 --- a/src/BiInts/bi_integrals.ezfio_config +++ b/src/BiInts/bi_integrals.ezfio_config @@ -5,4 +5,5 @@ bielec_integrals write_mo_integrals logical threshold_ao double precision threshold_mo double precision + direct logical diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f index 349886c3..27ff4959 100644 --- a/src/BiInts/mo_bi_integrals.irp.f +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -296,95 +296,132 @@ end 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 and - 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(:) + implicit none + BEGIN_DOC + ! Transform Bi-electronic integrals and + END_DOC + + integer :: i,j,p,q,r,s + double precision :: c + real(integral_kind) :: integral + integer :: n, pp + real(integral_kind), allocatable :: int_value(:) + integer, allocatable :: int_idx(:) + + double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) + + 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_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 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), & + iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& + iqsr(mo_tot_num_align,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=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 + + if (do_direct_integrals) then + double precision :: ao_bielec_integral + do r=1,ao_num + call compute_ao_bielec_integrals(q,r,s,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + 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 compute_ao_bielec_integrals(q,s,r,ao_num,int_value) + do p=1,ao_num + integral = int_value(p) + 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 - PROVIDE ao_bielec_integrals_in_map + else - 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 - - + 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 - + diff --git a/src/BiInts/options.irp.f b/src/BiInts/options.irp.f index c821661c..d4af85e1 100644 --- a/src/BiInts/options.irp.f +++ b/src/BiInts/options.irp.f @@ -107,3 +107,21 @@ BEGIN_PROVIDER [ double precision, mo_integrals_threshold ] 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 diff --git a/src/Hartree_Fock/Fock_matrix_mo.irp.f b/src/Hartree_Fock/Fock_matrix_mo.irp.f index b27d9799..9b5c193d 100644 --- a/src/Hartree_Fock/Fock_matrix_mo.irp.f +++ b/src/Hartree_Fock/Fock_matrix_mo.irp.f @@ -92,24 +92,40 @@ END_PROVIDER 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 + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_ints_idx, ao_ints_val + 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 i=1,ao_num + do i=1,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) 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) + 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) += 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 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 + !$OMP END DO NOWAIT + !$OMP END PARALLEL + else !$OMP PARALLEL DEFAULT(NONE) & !$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 HF_density_matrix_ao_beta) allocate(ao_ints_idx(ao_num_align),ao_ints_val(ao_num_align)) - !$OMP DO - do i=1,ao_num - do j=1,ao_num + !$OMP DO SCHEDULE(GUIDED) + 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_beta_ao (i,j) = ao_mono_elec_integral(i,j) 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) + !DIR$ VECTOR ALIGNED 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 - integral = -ao_ints_val(k1) - Fock_matrix_alpha_ao(i,l) += HF_density_matrix_ao_alpha(k,j) * integral - Fock_matrix_beta_ao (i,l) += HF_density_matrix_ao_beta (k,j) * integral + 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) += integral + Fock_matrix_beta_ao (i,j) += integral + integral = ao_ints_val(k1) + Fock_matrix_alpha_ao(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral + Fock_matrix_beta_ao (l,j) -= HF_density_matrix_ao_beta (k,i) * integral enddo enddo enddo diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 410339a1..b8981dd0 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -1,8 +1,3 @@ -=================== -Hartree-Fock Module -=================== - - Needed Modules ============== diff --git a/src/Hartree_Fock/hartree_fock.ezfio_config b/src/Hartree_Fock/hartree_fock.ezfio_config index 134486f4..0a53b7a9 100644 --- a/src/Hartree_Fock/hartree_fock.ezfio_config +++ b/src/Hartree_Fock/hartree_fock.ezfio_config @@ -1,6 +1,5 @@ hartree_fock thresh_scf double precision n_it_scf_max integer - direct logical diis logical diff --git a/src/Hartree_Fock/options.irp.f b/src/Hartree_Fock/options.irp.f index 303ce89d..74819a63 100644 --- a/src/Hartree_Fock/options.irp.f +++ b/src/Hartree_Fock/options.irp.f @@ -35,24 +35,6 @@ BEGIN_PROVIDER [ integer, n_it_scf_max] 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 ] implicit none BEGIN_DOC