From f18f96e76ec81bbc75deb6ac66fcf69f8706a74e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Oct 2014 14:59:03 +0200 Subject: [PATCH] Integrals AO not recalcaulated if not necessary --- src/BiInts/map_integrals.irp.f | 1 + src/BiInts/mo_bi_integrals.irp.f | 51 +++++++++++++++----- src/Dets/NEEDED_MODULES | 2 +- src/Dets/determinants.irp.f | 2 +- src/{Hartree_Fock => Dets}/ref_bitmask.irp.f | 30 +++--------- src/Dets/save_fock_orb.irp.f | 45 ----------------- src/Hartree_Fock/Fock_matrix.irp.f | 22 ++++++--- src/Hartree_Fock/diagonalize_fock.irp.f | 4 +- 8 files changed, 68 insertions(+), 89 deletions(-) rename src/{Hartree_Fock => Dets}/ref_bitmask.irp.f (76%) delete mode 100644 src/Dets/save_fock_orb.irp.f diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index a9214347..098f585d 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -382,6 +382,7 @@ subroutine dump_$ao_integrals(filename) integer(cache_key_kind), pointer :: key(:) real(integral_kind), pointer :: val(:) integer*8 :: i,j, n + call ezfio_set_work_empty(.False.) open(unit=66,file=filename,FORM='unformatted') write(66) integral_kind, key_kind write(66) $ao_integrals_map%sorted, $ao_integrals_map%map_size, & diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f index 6e89aa25..eb9e650d 100644 --- a/src/BiInts/mo_bi_integrals.irp.f +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -319,14 +319,14 @@ 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) ] + BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_exchange_from_ao, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti_from_ao, (mo_tot_num_align,mo_tot_num) ] implicit none BEGIN_DOC - ! mo_bielec_integral_jj(i,j) = J_ij - ! mo_bielec_integral_jj_exchange(i,j) = J_ij - ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij + ! mo_bielec_integral_jj_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij + ! mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij END_DOC integer :: i,j,p,q,r,s @@ -342,8 +342,8 @@ end PROVIDE ao_bielec_integrals_in_map endif - mo_bielec_integral_jj = 0.d0 - mo_bielec_integral_jj_exchange = 0.d0 + mo_bielec_integral_jj_from_ao = 0.d0 + mo_bielec_integral_jj_exchange_from_ao = 0.d0 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr @@ -353,7 +353,7 @@ end !$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,abort_here) & - !$OMP REDUCTION(+:mo_bielec_integral_jj,mo_bielec_integral_jj_exchange) + !$OMP REDUCTION(+:mo_bielec_integral_jj_from_ao,mo_bielec_integral_jj_exchange_from_ao) allocate( int_value(ao_num), int_idx(ao_num), & iqrs(mo_tot_num_align,ao_num), iqis(mo_tot_num), iqri(mo_tot_num),& @@ -439,8 +439,8 @@ end !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) + mo_bielec_integral_jj_from_ao(j,i) += c * iqis(i) + mo_bielec_integral_jj_exchange_from_ao(j,i) += c * iqri(i) enddo enddo @@ -453,8 +453,35 @@ end stop 'Aborting in MO integrals calculation' endif - mo_bielec_integral_jj_anti = mo_bielec_integral_jj - mo_bielec_integral_jj_exchange + mo_bielec_integral_jj_anti_from_ao = mo_bielec_integral_jj_from_ao - mo_bielec_integral_jj_exchange_from_ao 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 + ! mo_bielec_integral_jj(i,j) = J_ij + ! mo_bielec_integral_jj_exchange(i,j) = K_ij + ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij + END_DOC + + integer :: i,j + double precision :: get_mo_bielec_integral + + PROVIDE mo_bielec_integrals_in_map + + mo_bielec_integral_jj = 0.d0 + mo_bielec_integral_jj_exchange = 0.d0 + do j=1,mo_tot_num + do i=1,mo_tot_num + mo_bielec_integral_jj(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map) + mo_bielec_integral_jj_exchange(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map) + mo_bielec_integral_jj_anti(i,j) = mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + enddo + +END_PROVIDER diff --git a/src/Dets/NEEDED_MODULES b/src/Dets/NEEDED_MODULES index 26eb5a9a..163bdf10 100644 --- a/src/Dets/NEEDED_MODULES +++ b/src/Dets/NEEDED_MODULES @@ -1 +1 @@ -AOs BiInts Bitmask Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils +AOs BiInts Bitmask Electrons Ezfio_files MonoInts MOs Nuclei Output Utils diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 4ac8eda4..073654ba 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -352,7 +352,7 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,psicoef) call ezfio_set_determinants_psi_det(psi_det_save) deallocate (psi_det_save) - progress_bar(7) = 7 + progress_bar(1) = 7 progress_value = dble(progress_bar(1)) allocate (psi_coef_save(ndet,nstates)) do k=1,nstates diff --git a/src/Hartree_Fock/ref_bitmask.irp.f b/src/Dets/ref_bitmask.irp.f similarity index 76% rename from src/Hartree_Fock/ref_bitmask.irp.f rename to src/Dets/ref_bitmask.irp.f index 2b0e927f..7f760562 100644 --- a/src/Hartree_Fock/ref_bitmask.irp.f +++ b/src/Dets/ref_bitmask.irp.f @@ -1,4 +1,5 @@ - BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ] + BEGIN_PROVIDER [ double precision, ref_bitmask_energy ] +&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ] &BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ] &BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ] &BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ] @@ -15,17 +16,20 @@ call bitstring_to_list(ref_bitmask(1,2), occ(1,2), i, N_int) + ref_bitmask_energy = 0.d0 mono_elec_ref_bitmask_energy = 0.d0 kinetic_ref_bitmask_energy = 0.d0 nucl_elec_ref_bitmask_energy = 0.d0 bi_elec_ref_bitmask_energy = 0.d0 do i = 1, elec_beta_num + ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) + mo_mono_elec_integral(occ(i,2),occ(i,2)) kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) + mo_kinetic_integral(occ(i,2),occ(i,2)) nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) + mo_nucl_elec_integral(occ(i,2),occ(i,2)) enddo do i = elec_beta_num+1,elec_alpha_num + ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) enddo @@ -33,39 +37,21 @@ do j= 1, elec_alpha_num do i = j+1, elec_alpha_num bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1)) + ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1)) enddo enddo do j= 1, elec_beta_num do i = j+1, elec_beta_num bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2)) + ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2)) enddo do i= 1, elec_alpha_num bi_elec_ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2)) + ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2)) enddo enddo mono_elec_ref_bitmask_energy = kinetic_ref_bitmask_energy + nucl_elec_ref_bitmask_energy END_PROVIDER - BEGIN_PROVIDER [ double precision, ref_bitmask_energy ] - use bitmasks - implicit none - BEGIN_DOC - ! Energy of the reference bitmask used in Slater rules - END_DOC - - integer :: i,j - - ref_bitmask_energy = 0.d0 - - do j=1,ao_num - do i=1,ao_num - ref_bitmask_energy += 0.5d0 * ( & - (ao_mono_elec_integral(i,j) + Fock_matrix_alpha_ao(i,j) ) * HF_density_matrix_ao_alpha(i,j) + & - (ao_mono_elec_integral(i,j) + Fock_matrix_beta_ao (i,j) ) * HF_density_matrix_ao_beta (i,j) ) - enddo - enddo - -END_PROVIDER - diff --git a/src/Dets/save_fock_orb.irp.f b/src/Dets/save_fock_orb.irp.f deleted file mode 100644 index 67c27db2..00000000 --- a/src/Dets/save_fock_orb.irp.f +++ /dev/null @@ -1,45 +0,0 @@ -program save_fock_orb - - double precision, allocatable :: one_body_dm_ao_alpha(:,:) - double precision, allocatable :: one_body_dm_ao_beta (:,:) - - read_wf = .True. - touch read_wf - print *, 'N_det = ',N_det - print *, HF_energy - allocate (one_body_dm_ao_alpha(ao_num_align,ao_num), & - one_body_dm_ao_beta (ao_num_align,ao_num)) - - call mo_to_ao_no_overlap(one_body_dm_mo_alpha,mo_tot_num_align, & - one_body_dm_ao_alpha,ao_num_align) - call mo_to_ao_no_overlap(one_body_dm_mo_beta ,mo_tot_num_align, & - one_body_dm_ao_beta ,ao_num_align) - - do i=1,mo_tot_num - do j=1,mo_tot_num - if (abs(fock_matrix_mo(i,j)) > 1.d-10) then - print *, i,j, fock_matrix_mo(i,j) - endif - enddo - enddo - hf_density_matrix_ao_alpha = one_body_dm_ao_alpha - hf_density_matrix_ao_beta = one_body_dm_ao_beta - touch hf_density_matrix_ao_alpha hf_density_matrix_ao_beta - - print *, '---' - do i=1,mo_tot_num - do j=1,mo_tot_num - if (abs(fock_matrix_mo(i,j)) > 1.d-10) then - print *, i,j, fock_matrix_mo(i,j) - endif - enddo - enddo - mo_coef = eigenvectors_fock_matrix_mo - mo_label = 'CASSCF' - TOUCH mo_coef mo_label - print *, HF_energy - call save_mos - - deallocate(one_body_dm_ao_alpha,one_body_dm_ao_beta) -end - diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 667264f6..7dd349b1 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -272,7 +272,17 @@ BEGIN_PROVIDER [ double precision, HF_energy ] BEGIN_DOC ! Hartree-Fock energy END_DOC - HF_energy = nuclear_repulsion + ref_bitmask_energy + HF_energy = nuclear_repulsion + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + HF_energy += 0.5d0 * ( & + (ao_mono_elec_integral(i,j) + Fock_matrix_alpha_ao(i,j) ) * HF_density_matrix_ao_alpha(i,j) +& + (ao_mono_elec_integral(i,j) + Fock_matrix_beta_ao (i,j) ) * HF_density_matrix_ao_beta (i,j) ) + enddo + enddo + END_PROVIDER @@ -283,12 +293,12 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] END_DOC if (elec_alpha_num == elec_beta_num) then - integer :: i,j + integer :: i,j do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j) - enddo + !DIR$ VECTOR ALIGNED + do i=1,ao_num_align + Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j) + enddo enddo else double precision, allocatable :: T(:,:), M(:,:) diff --git a/src/Hartree_Fock/diagonalize_fock.irp.f b/src/Hartree_Fock/diagonalize_fock.irp.f index 9a396055..3ee0d298 100644 --- a/src/Hartree_Fock/diagonalize_fock.irp.f +++ b/src/Hartree_Fock/diagonalize_fock.irp.f @@ -66,14 +66,14 @@ BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] do i = 1,elec_alpha_num accu = 0.d0 do j = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) enddo diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i) enddo do i = elec_alpha_num+1,mo_tot_num accu = 0.d0 do j = 1, elec_alpha_num - accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j) enddo diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i) enddo