From ce042fbd787a21a600830596fa3caa5f7aa2cdb1 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 21 May 2024 12:01:28 +0200 Subject: [PATCH] basis set correction with cholesky works for hf --- .../local/basis_correction/51.basis_c.bats | 8 -- .../{01.convert.bats => convert_bats_old} | 0 src/hartree_fock/10.hf.bats | 13 -- src/mu_of_r/basis_def.irp.f | 45 +++++++ .../{f_cholesky.irp.f => f_hf_cholesky.irp.f} | 121 +++++++++--------- 5 files changed, 104 insertions(+), 83 deletions(-) rename src/ezfio_files/{01.convert.bats => convert_bats_old} (100%) rename src/mu_of_r/{f_cholesky.irp.f => f_hf_cholesky.irp.f} (67%) diff --git a/plugins/local/basis_correction/51.basis_c.bats b/plugins/local/basis_correction/51.basis_c.bats index 914b482b..1e20bae3 100644 --- a/plugins/local/basis_correction/51.basis_c.bats +++ b/plugins/local/basis_correction/51.basis_c.bats @@ -37,14 +37,6 @@ function run_sd() { eq $energy1 $1 $thresh } -@test "O2 CAS" { - qp set_file o2_cas.gms.ezfio - qp set_mo_class -c "[1-2]" -a "[3-10]" -d "[11-46]" - run -149.72435425 3.e-4 10000 - qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]" - run_md -0.1160222327 1.e-6 -} - @test "LiF RHF" { qp set_file lif.ezfio diff --git a/src/ezfio_files/01.convert.bats b/src/ezfio_files/convert_bats_old similarity index 100% rename from src/ezfio_files/01.convert.bats rename to src/ezfio_files/convert_bats_old diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index b496a089..214dfa86 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -115,9 +115,6 @@ rm -rf $EZFIO run hco.ezfio -113.1841002944744 } -@test "HBO" { # 0.805600 1.4543s - run hbo.ezfio -100.018582259096 -} @test "H2S" { # 1.655600 4.21402s run h2s.ezfio -398.6944130421982 @@ -127,9 +124,6 @@ rm -rf $EZFIO run h3coh.ezfio -114.9865030596373 } -@test "H2O" { # 1.811100 1.84387s - run h2o.ezfio -0.760270218692179E+02 -} @test "H2O2" { # 2.217000 8.50267s run h2o2.ezfio -150.7806608469964 @@ -187,13 +181,6 @@ rm -rf $EZFIO run oh.ezfio -75.42025413469165 } -@test "[Cu(NH3)4]2+" { # 59.610100 4.18766m - [[ -n $TRAVIS ]] && skip - qp set_file cu_nh3_4_2plus.ezfio - qp set scf_utils thresh_scf 1.e-10 - run cu_nh3_4_2plus.ezfio -1862.97590358903 -} - @test "SO2" { # 71.894900 3.22567m [[ -n $TRAVIS ]] && skip run so2.ezfio -41.55800401346361 diff --git a/src/mu_of_r/basis_def.irp.f b/src/mu_of_r/basis_def.irp.f index fff9f581..e433f4d8 100644 --- a/src/mu_of_r/basis_def.irp.f +++ b/src/mu_of_r/basis_def.irp.f @@ -114,3 +114,48 @@ BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_fi enddo enddo END_PROVIDER + +! BEGIN_PROVIDER [integer, n_docc_val_orb_for_cas] +!&BEGIN_PROVIDER [integer, n_max_docc_val_orb_for_cas] +! implicit none +! BEGIN_DOC +! ! Number of DOUBLY OCCUPIED VALENCE ORBITALS for the CAS wave function +! ! +! ! This determines the size of the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 +! END_DOC +! integer :: i +! n_docc_val_orb_for_cas = 0 +! ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL +! do i = 1, elec_beta_num +! if( trim(mo_class(i))=="Inactive" & +! .or. trim(mo_class(i))=="Active" & +! .or. trim(mo_class(i))=="Virtual" )then +! n_docc_val_orb_for_cas +=1 +! endif +! enddo +! n_max_docc_val_orb_for_cas = maxval(n_docc_val_orb_for_cas) +! +!END_PROVIDER +! +!BEGIN_PROVIDER [integer, list_doc_valence_orb_for_cas, (n_max_docc_val_orb_for_cas)] +! implicit none +! BEGIN_DOC +! ! List of OCCUPIED valence orbitals for each spin to build the f_{HF}(r_1,r_2) function +! ! +! ! This corresponds to ALL OCCUPIED orbitals in the HF wave function, except those defined as "core" +! ! +! ! This determines the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 +! END_DOC +! j = 0 +! ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL +! do i = 1, elec_beta_num +! if( trim(mo_class(i))=="Inactive" & +! .or. trim(mo_class(i))=="Active" & +! .or. trim(mo_class(i))=="Virtual" )then +! j +=1 +! list_doc_valence_orb_for_cas(j) = i +! endif +! enddo +! +!END_PROVIDER + diff --git a/src/mu_of_r/f_cholesky.irp.f b/src/mu_of_r/f_hf_cholesky.irp.f similarity index 67% rename from src/mu_of_r/f_cholesky.irp.f rename to src/mu_of_r/f_hf_cholesky.irp.f index 1ad4ce36..84097f09 100644 --- a/src/mu_of_r/f_cholesky.irp.f +++ b/src/mu_of_r/f_hf_cholesky.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [integer, list_couple_orb_r1, (2,n_couple_orb_r1)] +BEGIN_PROVIDER [integer, list_couple_hf_orb_r1, (2,n_couple_orb_r1)] implicit none integer :: ii,i,mm,m,itmp itmp = 0 @@ -7,14 +7,14 @@ BEGIN_PROVIDER [integer, list_couple_orb_r1, (2,n_couple_orb_r1)] do mm = 1, n_basis_orb ! electron 1 m = list_basis(mm) itmp += 1 - list_couple_orb_r1(1,itmp) = i - list_couple_orb_r1(2,itmp) = m + list_couple_hf_orb_r1(1,itmp) = i + list_couple_hf_orb_r1(2,itmp) = m enddo enddo END_PROVIDER -BEGIN_PROVIDER [integer, list_couple_orb_r2, (2,n_couple_orb_r2)] +BEGIN_PROVIDER [integer, list_couple_hf_orb_r2, (2,n_couple_orb_r2)] implicit none integer :: ii,i,mm,m,itmp itmp = 0 @@ -23,8 +23,8 @@ BEGIN_PROVIDER [integer, list_couple_orb_r2, (2,n_couple_orb_r2)] do mm = 1, n_basis_orb ! electron 1 m = list_basis(mm) itmp += 1 - list_couple_orb_r2(1,itmp) = i - list_couple_orb_r2(2,itmp) = m + list_couple_hf_orb_r2(1,itmp) = i + list_couple_hf_orb_r2(2,itmp) = m enddo enddo END_PROVIDER @@ -87,31 +87,6 @@ BEGIN_PROVIDER [ double precision, mos_times_cholesky_r1, (cholesky_mo_num,n_poi enddo call get_AB_prod(mo_chol_r1,cholesky_mo_num,n_couple_orb_r1,mos_ib_r1,n_points_final_grid,mos_times_cholesky_r1) - allocate(test(cholesky_mo_num,n_points_final_grid)) - test = 0.d0 - do ipoint = 1, n_points_final_grid - do itmp = 1, n_couple_orb_r1 - i = list_couple_orb_r1(1,itmp) - m = list_couple_orb_r1(2,itmp) - mo_i_r1 = mos_in_r_array_omp(i,ipoint) - mo_b_r1 = mos_in_r_array_omp(m,ipoint) - do mm = 1, cholesky_mo_num - test(mm,ipoint) += mo_i_r1 * mo_b_r1 * mo_chol_r1(mm,itmp) - enddo - enddo - enddo - double precision :: accu - accu = 0.d0 - do ipoint = 1, n_points_final_grid - do mm = 1, cholesky_mo_num - accu += dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint) ) - if(dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then - print*,'problem ! ',dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)) & - , mos_times_cholesky_r1(mm,ipoint) , test(mm,ipoint) - endif - enddo - enddo - print*,'accu = ',accu END_PROVIDER @@ -157,53 +132,72 @@ BEGIN_PROVIDER [ double precision, mos_times_cholesky_r2, (cholesky_mo_num,n_poi enddo call get_AB_prod(mo_chol_r2,cholesky_mo_num,n_couple_orb_r2,mos_ib_r2,n_points_final_grid,mos_times_cholesky_r2) - allocate(test(cholesky_mo_num,n_points_final_grid)) - test = 0.d0 - do ipoint = 1, n_points_final_grid - do itmp = 1, n_couple_orb_r2 - i = list_couple_orb_r2(1,itmp) - m = list_couple_orb_r2(2,itmp) - mo_i_r2 = mos_in_r_array_omp(i,ipoint) - mo_b_r2 = mos_in_r_array_omp(m,ipoint) - do mm = 1, cholesky_mo_num - test(mm,ipoint) += mo_i_r2 * mo_b_r2 * mo_chol_r2(mm,itmp) - enddo - enddo - enddo - double precision :: accu - accu = 0.d0 - do ipoint = 1, n_points_final_grid - do mm = 1, cholesky_mo_num - accu += dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint) ) - if(dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then - print*,'problem ! ',dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)) & - , mos_times_cholesky_r2(mm,ipoint) , test(mm,ipoint) - endif - enddo - enddo - print*,'accu = ',accu END_PROVIDER BEGIN_PROVIDER [ double precision, f_hf_cholesky, (n_points_final_grid)] implicit none - integer :: ipoint + integer :: ipoint,m,k !!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ !! = \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ !! = \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R) !! = \sum_A V_AR G_AR !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI - double precision :: u_dot_v - do ipoint = 1, n_points_final_grid - f_hf_cholesky(ipoint) = 2.D0 * u_dot_v(mos_times_cholesky_r2(1,ipoint),mos_times_cholesky_r1(1,ipoint),cholesky_mo_num) - enddo + double precision :: u_dot_v,wall0,wall1 + if(elec_alpha_num == elec_beta_num)then + provide mos_times_cholesky_r1 + print*,'providing f_hf_cholesky ...' + call wall_time(wall0) + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,m) & + !$OMP ShARED (mos_times_cholesky_r1,cholesky_mo_num,f_hf_cholesky,n_points_final_grid) + do ipoint = 1, n_points_final_grid + f_hf_cholesky(ipoint) = 0.d0 + do m = 1, cholesky_mo_num + f_hf_cholesky(ipoint) = f_hf_cholesky(ipoint) + & + mos_times_cholesky_r1(m,ipoint) * mos_times_cholesky_r1(m,ipoint) + enddo + f_hf_cholesky(ipoint) *= 2.D0 + enddo + !$OMP END PARALLEL DO + + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky = ',wall1-wall0 + free mos_times_cholesky_r1 + else + provide mos_times_cholesky_r2 mos_times_cholesky_r1 + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,m) & + !$OMP ShARED (mos_times_cholesky_r2,mos_times_cholesky_r1,cholesky_mo_num,f_hf_cholesky,n_points_final_grid) + do ipoint = 1, n_points_final_grid + f_hf_cholesky(ipoint) = 0.D0 + do m = 1, cholesky_mo_num + f_hf_cholesky(ipoint) = f_hf_cholesky(ipoint) + & + mos_times_cholesky_r2(m,ipoint)*mos_times_cholesky_r1(m,ipoint) + enddo + f_hf_cholesky(ipoint) *= 2.D0 + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky = ',wall1-wall0 + free mos_times_cholesky_r2 mos_times_cholesky_r1 + endif END_PROVIDER BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] implicit none integer :: ipoint,i,ii - double precision :: dm_a, dm_b + double precision :: dm_a, dm_b,wall0,wall1 + print*,'providing on_top_hf_grid ...' + provide mos_in_r_array_omp + call wall_time(wall0) + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,dm_a,dm_b,ii,i) & + !$OMP ShARED (n_points_final_grid,n_occ_val_orb_for_hf,mos_in_r_array_omp,list_valence_orb_for_hf,on_top_hf_grid) do ipoint = 1, n_points_final_grid dm_a = 0.d0 do ii = 1, n_occ_val_orb_for_hf(1) @@ -217,5 +211,8 @@ BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] enddo on_top_hf_grid(ipoint) = 2.D0 * dm_a*dm_b enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide on_top_hf_grid = ',wall1-wall0 END_PROVIDER