mirror of
https://github.com/LCPQ/quantum_package
synced 2024-10-19 22:41:48 +02:00
added kohn_sham_range_separated
This commit is contained in:
parent
ffee456d62
commit
e2c380a1f3
1
src/kohn_sham_range_separated/NEED
Normal file
1
src/kohn_sham_range_separated/NEED
Normal file
@ -0,0 +1 @@
|
||||
DFT_Utils_one_body Integrals_erf SCF_Utils Integrals_Bielec
|
294
src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f
Normal file
294
src/kohn_sham_range_separated/fock_matrix_rs_ks.irp.f
Normal file
@ -0,0 +1,294 @@
|
||||
BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_alpha, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num, ao_num) ]
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha Fock matrix in AO basis set
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l,k1,r,s
|
||||
integer :: i0,j0,k0,l0
|
||||
integer*8 :: p,q
|
||||
double precision :: integral, c0, c1, c2
|
||||
double precision :: ao_bielec_integral, local_threshold
|
||||
double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:)
|
||||
double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_alpha_tmp
|
||||
|
||||
ao_bi_elec_integral_alpha = 0.d0
|
||||
ao_bi_elec_integral_beta = 0.d0
|
||||
if (do_direct_integrals) then
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, &
|
||||
!$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, &
|
||||
!$OMP local_threshold)&
|
||||
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
|
||||
!$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, &
|
||||
!$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta)
|
||||
|
||||
allocate(keys(1), values(1))
|
||||
allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), &
|
||||
ao_bi_elec_integral_beta_tmp(ao_num,ao_num))
|
||||
ao_bi_elec_integral_alpha_tmp = 0.d0
|
||||
ao_bi_elec_integral_beta_tmp = 0.d0
|
||||
|
||||
q = ao_num*ao_num*ao_num*ao_num
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do p=1_8,q
|
||||
call bielec_integrals_index_reverse(kk,ii,ll,jj,p)
|
||||
if ( (kk(1)>ao_num).or. &
|
||||
(ii(1)>ao_num).or. &
|
||||
(jj(1)>ao_num).or. &
|
||||
(ll(1)>ao_num) ) then
|
||||
cycle
|
||||
endif
|
||||
k = kk(1)
|
||||
i = ii(1)
|
||||
l = ll(1)
|
||||
j = jj(1)
|
||||
|
||||
if (ao_overlap_abs(k,l)*ao_overlap_abs(i,j) &
|
||||
< ao_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j)
|
||||
if (local_threshold < ao_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
i0 = i
|
||||
j0 = j
|
||||
k0 = k
|
||||
l0 = l
|
||||
values(1) = 0.d0
|
||||
local_threshold = ao_integrals_threshold/local_threshold
|
||||
do k2=1,8
|
||||
if (kk(k2)==0) then
|
||||
cycle
|
||||
endif
|
||||
i = ii(k2)
|
||||
j = jj(k2)
|
||||
k = kk(k2)
|
||||
l = ll(k2)
|
||||
c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)
|
||||
c1 = SCF_density_matrix_ao_alpha(k,i)
|
||||
c2 = SCF_density_matrix_ao_beta(k,i)
|
||||
if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then
|
||||
cycle
|
||||
endif
|
||||
if (values(1) == 0.d0) then
|
||||
values(1) = ao_bielec_integral(k0,l0,i0,j0)
|
||||
endif
|
||||
integral = c0 * values(1)
|
||||
ao_bi_elec_integral_alpha_tmp(i,j) += integral
|
||||
ao_bi_elec_integral_beta_tmp (i,j) += integral
|
||||
integral = values(1)
|
||||
ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral
|
||||
ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp
|
||||
!$OMP END CRITICAL
|
||||
deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)
|
||||
!$OMP END PARALLEL
|
||||
else
|
||||
PROVIDE ao_bielec_integrals_in_map
|
||||
PROVIDE ao_bielec_integrals_erf_in_map
|
||||
|
||||
integer(omp_lock_kind) :: lck(ao_num)
|
||||
integer*8 :: i8
|
||||
integer :: ii(8), jj(8), kk(8), ll(8), k2
|
||||
integer(cache_map_size_kind) :: n_elements_max, n_elements
|
||||
integer(key_kind), allocatable :: keys(:)
|
||||
double precision, allocatable :: values(:)
|
||||
integer(cache_map_size_kind) :: n_elements_max_erf, n_elements_erf
|
||||
integer(key_kind), allocatable :: keys_erf(:)
|
||||
double precision, allocatable :: values_erf(:)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, &
|
||||
!$OMP n_elements,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)&
|
||||
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
|
||||
!$OMP ao_integrals_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta)
|
||||
|
||||
call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
|
||||
allocate(keys(n_elements_max), values(n_elements_max))
|
||||
allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), &
|
||||
ao_bi_elec_integral_beta_tmp(ao_num,ao_num))
|
||||
ao_bi_elec_integral_alpha_tmp = 0.d0
|
||||
ao_bi_elec_integral_beta_tmp = 0.d0
|
||||
|
||||
!$OMP DO SCHEDULE(dynamic,64)
|
||||
!DIR$ NOVECTOR
|
||||
do i8=0_8,ao_integrals_map%map_size
|
||||
n_elements = n_elements_max
|
||||
call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
|
||||
do k1=1,n_elements
|
||||
call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
|
||||
|
||||
do k2=1,8
|
||||
if (kk(k2)==0) then
|
||||
cycle
|
||||
endif
|
||||
i = ii(k2)
|
||||
j = jj(k2)
|
||||
k = kk(k2)
|
||||
l = ll(k2)
|
||||
integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1)
|
||||
ao_bi_elec_integral_alpha_tmp(i,j) += integral
|
||||
ao_bi_elec_integral_beta_tmp (i,j) += integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_alpha += ao_bi_elec_integral_alpha_tmp
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_beta += ao_bi_elec_integral_beta_tmp
|
||||
!$OMP END CRITICAL
|
||||
deallocate(keys,values,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i,j,l,k1,k,integral_erf,ii,jj,kk,ll,i8,keys_erf,values_erf,n_elements_max_erf, &
|
||||
!$OMP n_elements_erf,ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)&
|
||||
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
|
||||
!$OMP ao_integrals_erf_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta)
|
||||
|
||||
|
||||
call get_cache_map_n_elements_max(ao_integrals_erf_map,n_elements_max_erf)
|
||||
allocate(ao_bi_elec_integral_alpha_tmp(ao_num,ao_num), &
|
||||
ao_bi_elec_integral_beta_tmp(ao_num,ao_num))
|
||||
allocate(keys_Erf(n_elements_max_erf), values_erf(n_elements_max_erf))
|
||||
|
||||
ao_bi_elec_integral_alpha_tmp = 0.d0
|
||||
ao_bi_elec_integral_beta_tmp = 0.d0
|
||||
!$OMP DO SCHEDULE(dynamic,64)
|
||||
!DIR$ NOVECTOR
|
||||
do i8=0_8,ao_integrals_erf_map%map_size
|
||||
n_elements_erf = n_elements_max_erf
|
||||
call get_cache_map(ao_integrals_erf_map,i8,keys_erf,values_erf,n_elements_erf)
|
||||
do k1=1,n_elements_erf
|
||||
call bielec_integrals_index_reverse(kk,ii,ll,jj,keys_erf(k1))
|
||||
|
||||
do k2=1,8
|
||||
if (kk(k2)==0) then
|
||||
cycle
|
||||
endif
|
||||
i = ii(k2)
|
||||
j = jj(k2)
|
||||
k = kk(k2)
|
||||
l = ll(k2)
|
||||
double precision :: integral_erf
|
||||
integral_erf = values_erf(k1)
|
||||
ao_bi_elec_integral_alpha_tmp(l,j) -= (SCF_density_matrix_ao_alpha(k,i) * integral_erf)
|
||||
ao_bi_elec_integral_beta_tmp (l,j) -= (SCF_density_matrix_ao_beta (k,i) * integral_erf)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_alpha = ao_bi_elec_integral_alpha + ao_bi_elec_integral_alpha_tmp
|
||||
!$OMP END CRITICAL
|
||||
!$OMP CRITICAL
|
||||
ao_bi_elec_integral_beta = ao_bi_elec_integral_beta + ao_bi_elec_integral_beta_tmp
|
||||
!$OMP END CRITICAL
|
||||
deallocate(ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)
|
||||
deallocate(keys_erf,values_erf)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha Fock matrix in AO basis set
|
||||
END_DOC
|
||||
|
||||
integer :: i,j
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
Fock_matrix_ao_alpha(i,j) = Fock_matrix_alpha_no_xc_ao(i,j) + ao_potential_alpha_xc(i,j)
|
||||
Fock_matrix_ao_beta (i,j) = Fock_matrix_beta_no_xc_ao(i,j) + ao_potential_beta_xc(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_no_xc_ao, (ao_num, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_no_xc_ao, (ao_num, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Mono electronic an Coulomb matrix in AO basis set
|
||||
END_DOC
|
||||
|
||||
integer :: i,j
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
Fock_matrix_alpha_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j)
|
||||
Fock_matrix_beta_no_xc_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, RS_KS_energy ]
|
||||
!BEGIN_PROVIDER [ double precision, SCF_energy ]
|
||||
&BEGIN_PROVIDER [ double precision, two_electron_energy]
|
||||
&BEGIN_PROVIDER [ double precision, one_electron_energy]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_energy]
|
||||
&BEGIN_PROVIDER [ double precision, trace_potential_xc ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Range-separated Kohn-Sham energy
|
||||
END_DOC
|
||||
RS_KS_energy = nuclear_repulsion
|
||||
|
||||
integer :: i,j
|
||||
double precision :: accu_mono,accu_fock
|
||||
one_electron_energy = 0.d0
|
||||
two_electron_energy = 0.d0
|
||||
Fock_matrix_energy = 0.d0
|
||||
trace_potential_xc = 0.d0
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
Fock_matrix_energy += Fock_matrix_ao_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) + &
|
||||
Fock_matrix_ao_beta(i,j) * SCF_density_matrix_ao_beta(i,j)
|
||||
two_electron_energy += 0.5d0 * ( ao_bi_elec_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
|
||||
+ao_bi_elec_integral_beta(i,j) * SCF_density_matrix_ao_beta(i,j) )
|
||||
one_electron_energy += ao_mono_elec_integral(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
|
||||
! possible bug fix for open-shell
|
||||
! trace_potential_xc += (ao_potential_alpha_xc(i,j) + ao_potential_beta_xc(i,j) ) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
|
||||
trace_potential_xc += ao_potential_alpha_xc(i,j) * SCF_density_matrix_ao_alpha(i,j) + ao_potential_beta_xc(i,j) * SCF_density_matrix_ao_beta (i,j)
|
||||
enddo
|
||||
enddo
|
||||
RS_KS_energy += e_exchange_dft + e_correlation_dft + one_electron_energy + two_electron_energy
|
||||
!SCF_energy = RS_KS_energy
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, extra_energy_contrib_from_density]
|
||||
implicit none
|
||||
! possible bug fix for open-shell:
|
||||
! extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.25d0 * trace_potential_xc
|
||||
extra_energy_contrib_from_density = e_exchange_dft + e_correlation_dft - 0.5d0 * trace_potential_xc
|
||||
END_PROVIDER
|
||||
|
||||
!BEGIN_PROVIDER [ double precision, SCF_energy ]
|
||||
! implicit none
|
||||
! SCF_energy = RS_KS_energy
|
||||
!END_PROVIDER
|
25
src/kohn_sham_range_separated/potential_functional.irp.f
Normal file
25
src/kohn_sham_range_separated/potential_functional.irp.f
Normal file
@ -0,0 +1,25 @@
|
||||
BEGIN_PROVIDER [double precision, ao_potential_alpha_xc, (ao_num, ao_num)]
|
||||
&BEGIN_PROVIDER [double precision, ao_potential_beta_xc, (ao_num, ao_num)]
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
ao_potential_alpha_xc = 0.d0
|
||||
ao_potential_beta_xc = 0.d0
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
ao_potential_alpha_xc(i,j) = potential_c_alpha_ao(i,j,1) + potential_x_alpha_ao(i,j,1)
|
||||
ao_potential_beta_xc(i,j) = potential_c_beta_ao(i,j,1) + potential_x_beta_ao(i,j,1)
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, e_exchange_dft]
|
||||
implicit none
|
||||
e_exchange_dft = energy_x(1)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, e_correlation_dft]
|
||||
implicit none
|
||||
e_correlation_dft = energy_c(1)
|
||||
|
||||
END_PROVIDER
|
103
src/kohn_sham_range_separated/rs_ks_scf.irp.f
Normal file
103
src/kohn_sham_range_separated/rs_ks_scf.irp.f
Normal file
@ -0,0 +1,103 @@
|
||||
program scf
|
||||
BEGIN_DOC
|
||||
! Produce `Kohn_Sham` MO orbital
|
||||
! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ
|
||||
! output: kohn_sham.energy
|
||||
! optional: mo_basis.mo_coef
|
||||
END_DOC
|
||||
read_wf = .False.
|
||||
density_for_dft ="WFT"
|
||||
touch density_for_dft
|
||||
touch read_wf
|
||||
print*, '**************************'
|
||||
print*, 'mu_erf = ',mu_erf
|
||||
print*, '**************************'
|
||||
call check_coherence_functional
|
||||
call create_guess
|
||||
call orthonormalize_mos
|
||||
call run
|
||||
end
|
||||
|
||||
subroutine check_coherence_functional
|
||||
implicit none
|
||||
integer :: ifound_x,ifound_c
|
||||
if(exchange_functional.eq."None")then
|
||||
ifound_x = 1
|
||||
else
|
||||
ifound_x = index(exchange_functional,"short_range")
|
||||
endif
|
||||
|
||||
if(correlation_functional.eq."None")then
|
||||
ifound_c = 1
|
||||
else
|
||||
ifound_c = index(correlation_functional,"short_range")
|
||||
endif
|
||||
print*,ifound_x,ifound_c
|
||||
if(ifound_x .eq.0 .or. ifound_c .eq. 0)then
|
||||
print*,'YOU ARE USING THE RANGE SEPARATED KS PROGRAM BUT YOUR INPUT KEYWORD FOR '
|
||||
print*,'exchange_functional is ',exchange_functional
|
||||
print*,'correlation_functional is ',correlation_functional
|
||||
print*,'CHANGE THE exchange_functional and correlation_functional keywords to range separated functionals'
|
||||
print*,'or switch to the KS_SCF program that uses regular functionals'
|
||||
stop
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine create_guess
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Create a MO guess if no MOs are present in the EZFIO directory
|
||||
END_DOC
|
||||
logical :: exists
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_mo_basis_mo_coef(exists)
|
||||
if (.not.exists) then
|
||||
print*,'Creating a guess for the MOs'
|
||||
print*,'mo_guess_type = ',mo_guess_type
|
||||
if (mo_guess_type == "HCore") then
|
||||
mo_coef = ao_ortho_lowdin_coef
|
||||
TOUCH mo_coef
|
||||
mo_label = 'Guess'
|
||||
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
|
||||
SOFT_TOUCH mo_coef mo_label
|
||||
else if (mo_guess_type == "Huckel") then
|
||||
call huckel_guess
|
||||
else
|
||||
print *, 'Unrecognized MO guess type : '//mo_guess_type
|
||||
stop 1
|
||||
endif
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine run
|
||||
|
||||
BEGIN_DOC
|
||||
! Run SCF calculation
|
||||
END_DOC
|
||||
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
double precision :: EHF
|
||||
|
||||
EHF = RS_KS_energy
|
||||
|
||||
mo_label = "Canonical"
|
||||
|
||||
! Choose SCF algorithm
|
||||
|
||||
! call damping_SCF ! Deprecated routine
|
||||
call Roothaan_Hall_SCF
|
||||
|
||||
write(*, '(A22,X,F16.10)') 'one_electron_energy = ',one_electron_energy
|
||||
write(*, '(A22,X,F16.10)') 'two_electron_energy = ',two_electron_energy
|
||||
write(*, '(A22,X,F16.10)') 'e_exchange_dft = ',e_exchange_dft
|
||||
write(*, '(A22,X,F16.10)') 'e_correlation_dft = ',e_correlation_dft
|
||||
write(*, '(A22,X,F16.10)') 'Fock_matrix_energy = ',Fock_matrix_energy
|
||||
|
||||
|
||||
end
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user