mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-04 13:13:53 +01:00
basis set correction with cholesky works for hf
This commit is contained in:
eginer
parent
c6a6163944
commit
ce042fbd78
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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
|
||||
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) = 2.D0 * u_dot_v(mos_times_cholesky_r2(1,ipoint),mos_times_cholesky_r1(1,ipoint),cholesky_mo_num)
|
||||
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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user