From 1df5dced1e3e6dc968cdd2e098cbdc9735371493 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 14 Apr 2017 17:33:35 +0200 Subject: [PATCH] KS LDA is okay --- plugins/DFT_Utils/functional.irp.f | 43 +- plugins/Kohn_Sham/EZFIO.cfg | 54 +++ plugins/Kohn_Sham/Fock_matrix.irp.f | 468 +++++++++++++++++++ plugins/Kohn_Sham/HF_density_matrix_ao.irp.f | 41 ++ plugins/Kohn_Sham/KS_SCF.irp.f | 54 +++ plugins/Kohn_Sham/NEEDED_CHILDREN_MODULES | 1 + plugins/Kohn_Sham/damping_SCF.irp.f | 132 ++++++ plugins/Kohn_Sham/diagonalize_fock.irp.f | 119 +++++ plugins/Kohn_Sham/potential_functional.irp.f | 31 ++ src/AO_Basis/aos_value.irp.f | 1 + src/Utils/constants.include.F | 1 + 11 files changed, 938 insertions(+), 7 deletions(-) create mode 100644 plugins/Kohn_Sham/EZFIO.cfg create mode 100644 plugins/Kohn_Sham/Fock_matrix.irp.f create mode 100644 plugins/Kohn_Sham/HF_density_matrix_ao.irp.f create mode 100644 plugins/Kohn_Sham/KS_SCF.irp.f create mode 100644 plugins/Kohn_Sham/NEEDED_CHILDREN_MODULES create mode 100644 plugins/Kohn_Sham/damping_SCF.irp.f create mode 100644 plugins/Kohn_Sham/diagonalize_fock.irp.f create mode 100644 plugins/Kohn_Sham/potential_functional.irp.f diff --git a/plugins/DFT_Utils/functional.irp.f b/plugins/DFT_Utils/functional.irp.f index 8bcb7d80..e034a244 100644 --- a/plugins/DFT_Utils/functional.irp.f +++ b/plugins/DFT_Utils/functional.irp.f @@ -1,25 +1,54 @@ -double precision function ex_lda(rho) +subroutine ex_lda(rho_a,rho_b,ex,vx_a,vx_b) include 'constants.include.F' implicit none - double precision, intent(in) :: rho - ex_lda = cst_lda * rho**(c_4_3) + double precision, intent(in) :: rho_a,rho_b + double precision, intent(out) :: ex,vx_a,vx_b + double precision :: tmp_a,tmp_b + tmp_a = rho_a**(c_1_3) + tmp_b = rho_b**(c_1_3) + ex = cst_lda * (tmp_a*tmp_a*tmp_a*tmp_a + tmp_b*tmp_b*tmp_b*tmp_b) + vx_a = cst_lda * c_4_3 * tmp_a + vx_b = cst_lda * c_4_3 * tmp_b end -BEGIN_PROVIDER [double precision, lda_exchange, (N_states)] + BEGIN_PROVIDER [double precision, lda_exchange, (N_states)] +&BEGIN_PROVIDER [double precision, lda_ex_potential_alpha_ao,(ao_num_align,ao_num,N_states)] +&BEGIN_PROVIDER [double precision, lda_ex_potential_beta_ao,(ao_num_align,ao_num,N_states)] + implicit none integer :: i,j,k,l - double precision :: ex_lda + integer :: m,n + double precision :: aos_array(ao_num) + double precision :: r(3) + lda_ex_potential_alpha_ao = 0.d0 + lda_ex_potential_beta_ao = 0.d0 do l = 1, N_states lda_exchange(l) = 0.d0 do j = 1, nucl_num do i = 1, n_points_radial_grid do k = 1, n_points_integration_angular - lda_exchange(l) += final_weight_functions_at_grid_points(k,i,j) * & - (ex_lda(one_body_dm_mo_alpha_at_grid_points(k,i,j,l)) + ex_lda(one_body_dm_mo_beta_at_grid_points(k,i,j,l))) + double precision :: rho_a,rho_b,ex + double precision :: vx_a,vx_b + rho_a = one_body_dm_mo_alpha_at_grid_points(k,i,j,l) + rho_b = one_body_dm_mo_beta_at_grid_points(k,i,j,l) + call ex_lda(rho_a,rho_b,ex,vx_a,vx_b) + lda_exchange(l) += final_weight_functions_at_grid_points(k,i,j) * ex + r(1) = grid_points_per_atom(1,k,i,j) + r(2) = grid_points_per_atom(2,k,i,j) + r(3) = grid_points_per_atom(3,k,i,j) + call give_all_aos_at_r(r,aos_array) + do m = 1, ao_num +! lda_ex_potential_ao(m,m,l) += (vx_a + vx_b) * aos_array(m)*aos_array(m) + do n = 1, ao_num + lda_ex_potential_alpha_ao(m,n,l) += (vx_a ) * aos_array(m)*aos_array(n) * final_weight_functions_at_grid_points(k,i,j) + lda_ex_potential_beta_ao(m,n,l) += (vx_b) * aos_array(m)*aos_array(n) * final_weight_functions_at_grid_points(k,i,j) + enddo + enddo enddo enddo enddo enddo END_PROVIDER + diff --git a/plugins/Kohn_Sham/EZFIO.cfg b/plugins/Kohn_Sham/EZFIO.cfg new file mode 100644 index 00000000..33d3a793 --- /dev/null +++ b/plugins/Kohn_Sham/EZFIO.cfg @@ -0,0 +1,54 @@ +[thresh_scf] +type: Threshold +doc: Threshold on the convergence of the Hartree Fock energy +interface: ezfio,provider,ocaml +default: 1.e-10 + +[exchange_functional] +type: character*(256) +doc: name of the exchange functional +interface: ezfio, provider, ocaml +default: "LDA" + + +[correlation_functional] +type: character*(256) +doc: name of the correlation functional +interface: ezfio, provider, ocaml +default: "LDA" + +[HF_exchange] +type: double precision +doc: Percentage of HF exchange in the DFT model +interface: ezfio,provider,ocaml +default: 0. + +[n_it_scf_max] +type: Strictly_positive_int +doc: Maximum number of SCF iterations +interface: ezfio,provider,ocaml +default: 200 + +[level_shift] +type: Positive_float +doc: Energy shift on the virtual MOs to improve SCF convergence +interface: ezfio,provider,ocaml +default: 0.5 + +[mo_guess_type] +type: MO_guess +doc: Initial MO guess. Can be [ Huckel | HCore ] +interface: ezfio,provider,ocaml +default: Huckel + +[energy] +type: double precision +doc: Calculated HF energy +interface: ezfio + +[no_oa_or_av_opt] +type: logical +doc: If true, skip the (inactive+core) --> (active) and the (active) --> (virtual) orbital rotations within the SCF procedure +interface: ezfio,provider,ocaml +default: False + diff --git a/plugins/Kohn_Sham/Fock_matrix.irp.f b/plugins/Kohn_Sham/Fock_matrix.irp.f new file mode 100644 index 00000000..9c91ddc9 --- /dev/null +++ b/plugins/Kohn_Sham/Fock_matrix.irp.f @@ -0,0 +1,468 @@ + BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + integer :: i,j,n + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_mo = Fock_matrix_alpha_mo + else + + do j=1,elec_beta_num + ! F-K + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))& + - (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + ! F+K/2 + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))& + + 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + ! F + do i=elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) + enddo + enddo + + do j=elec_beta_num+1,elec_alpha_num + ! F+K/2 + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))& + + 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + ! F + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) + enddo + ! F-K/2 + do i=elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))& + - 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + enddo + + do j=elec_alpha_num+1, mo_tot_num + ! F + do i=1,elec_beta_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) + enddo + ! F-K/2 + do i=elec_beta_num+1,elec_alpha_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))& + - 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + ! F+K + do i=elec_alpha_num+1,mo_tot_num + Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) & + + (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j)) + enddo + enddo + + endif + + do i = 1, mo_tot_num + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + enddo +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Alpha Fock matrix in AO basis set + END_DOC + + integer :: i,j + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,ao_num + Fock_matrix_alpha_ao(i,j) = Fock_matrix_alpha_no_xc_ao(i,j) + ao_potential_alpha_xc(i,j) + Fock_matrix_beta_ao (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_align, ao_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_no_xc_ao, (ao_num_align, 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 + !DIR$ VECTOR ALIGNED + 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, ao_bi_elec_integral_alpha, (ao_num_align, ao_num) ] +&BEGIN_PROVIDER [ double precision, ao_bi_elec_integral_beta , (ao_num_align, 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,ao_num_align,HF_density_matrix_ao_alpha,HF_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_align,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num_align,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 = HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l) + c1 = HF_density_matrix_ao_alpha(k,i) + c2 = HF_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 + + 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(:) + +! !$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,ao_num_align,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,& +! !$OMP ao_integrals_map, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta,HF_exchange) + + 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_align,ao_num), & + ao_bi_elec_integral_beta_tmp(ao_num_align,ao_num)) + ao_bi_elec_integral_alpha_tmp = 0.d0 + ao_bi_elec_integral_beta_tmp = 0.d0 + +! !OMP DO SCHEDULE(dynamic) +! !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 = (HF_density_matrix_ao_alpha(k,l)+HF_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 + integral = values(k1) + ao_bi_elec_integral_alpha_tmp(l,j) -= HF_exchange * (HF_density_matrix_ao_alpha(k,i) * integral) + ao_bi_elec_integral_beta_tmp (l,j) -= HF_exchange * (HF_density_matrix_ao_beta (k,i) * 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 + + endif + +END_PROVIDER + + + + + + +BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + double precision, allocatable :: T(:,:) + allocate ( T(ao_num_align,mo_tot_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, Fock_matrix_alpha_ao,size(Fock_matrix_alpha_ao,1), & + mo_coef, size(mo_coef,1), & + 0.d0, T, ao_num_align) + call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & + 1.d0, mo_coef,size(mo_coef,1), & + T, size(T,1), & + 0.d0, Fock_matrix_alpha_mo, mo_tot_num_align) + deallocate(T) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + double precision, allocatable :: T(:,:) + allocate ( T(ao_num_align,mo_tot_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + call dgemm('N','N', ao_num, mo_tot_num, ao_num, & + 1.d0, Fock_matrix_beta_ao,size(Fock_matrix_beta_ao,1), & + mo_coef, size(mo_coef,1), & + 0.d0, T, ao_num_align) + call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, & + 1.d0, mo_coef,size(mo_coef,1), & + T, size(T,1), & + 0.d0, Fock_matrix_beta_mo, mo_tot_num_align) + deallocate(T) +END_PROVIDER + + BEGIN_PROVIDER [ double precision, HF_energy ] +&BEGIN_PROVIDER [ double precision, two_electron_energy] +&BEGIN_PROVIDER [ double precision, one_electron_energy] + implicit none + BEGIN_DOC + ! Hartree-Fock energy + END_DOC + HF_energy = nuclear_repulsion + + integer :: i,j + double precision :: accu_mono,accu_fock + one_electron_energy = 0.d0 + two_electron_energy = 0.d0 + do j=1,ao_num + do i=1,ao_num + two_electron_energy += 0.5d0 * ( ao_bi_elec_integral_alpha(i,j) * HF_density_matrix_ao_alpha(i,j) & + +ao_bi_elec_integral_beta(i,j) * HF_density_matrix_ao_beta(i,j) ) + one_electron_energy += ao_mono_elec_integral(i,j) * (HF_density_matrix_ao_alpha(i,j) + HF_density_matrix_ao_beta (i,j) ) + enddo + enddo + print*, 'one_electron_energy = ',one_electron_energy + print*, 'two_electron_energy = ',two_electron_energy + print*, 'e_exchange_dft = ',(1.d0 - HF_exchange) * e_exchange_dft +!print*, 'accu_cor = ',e_correlation_dft + HF_energy += (1.d0 - HF_exchange) * e_exchange_dft + e_correlation_dft + one_electron_energy + two_electron_energy +!print*, 'HF_energy ' + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + 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 + enddo + else + double precision, allocatable :: T(:,:), M(:,:) + integer :: ierr + ! F_ao = S C F_mo C^t S + allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr) + if (ierr /=0 ) then + print *, irp_here, ' : allocation failed' + endif + +! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) +! -> M(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + +! M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num) +! -> T(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + Fock_matrix_mo, size(Fock_matrix_mo,1), & + 0.d0, & + T, size(T,1)) + +! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) +! -> M(ao_num,ao_num) + call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + +! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) +! -> Fock_matrix_ao(ao_num,ao_num) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + Fock_matrix_ao, size(Fock_matrix_ao,1)) + + + deallocate(T) + endif +END_PROVIDER + +subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + integer :: ierr + ! F_ao = S C F_mo C^t S + allocate (T(ao_num_align,ao_num),M(ao_num_align,ao_num),stat=ierr) + if (ierr /=0 ) then + print *, irp_here, ' : allocation failed' + endif + +! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num) +! -> M(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + +! M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num) +! -> T(ao_num,mo_tot_num) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + +! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num) +! -> M(ao_num,ao_num) + call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + +! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num) +! -> Fock_matrix_ao(ao_num,ao_num) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) +end + diff --git a/plugins/Kohn_Sham/HF_density_matrix_ao.irp.f b/plugins/Kohn_Sham/HF_density_matrix_ao.irp.f new file mode 100644 index 00000000..e8585f59 --- /dev/null +++ b/plugins/Kohn_Sham/HF_density_matrix_ao.irp.f @@ -0,0 +1,41 @@ +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! S^-1 x Alpha density matrix in the AO basis x S^-1 + END_DOC + + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! S^-1 Beta density matrix in the AO basis x S^-1 + END_DOC + + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! S^-1 Density matrix in the AO basis S^-1 + END_DOC + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) + if (elec_alpha_num== elec_beta_num) then + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha + else + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1)) + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta + endif + +END_PROVIDER + diff --git a/plugins/Kohn_Sham/KS_SCF.irp.f b/plugins/Kohn_Sham/KS_SCF.irp.f new file mode 100644 index 00000000..dead61ee --- /dev/null +++ b/plugins/Kohn_Sham/KS_SCF.irp.f @@ -0,0 +1,54 @@ +program scf + BEGIN_DOC +! Produce `Hartree_Fock` 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: hartree_fock.energy +! optional: mo_basis.mo_coef + END_DOC + call create_guess + call orthonormalize_mos + call run +end + +subroutine create_guess + implicit none + BEGIN_DOC +! Create an 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 + 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) + 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 + + use bitmasks + implicit none + BEGIN_DOC +! Run SCF calculation + END_DOC + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem + double precision :: E0 + integer :: i_it, i, j, k + + E0 = HF_energy + + mo_label = "Canonical" + call damping_SCF + +end diff --git a/plugins/Kohn_Sham/NEEDED_CHILDREN_MODULES b/plugins/Kohn_Sham/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..d8c28b56 --- /dev/null +++ b/plugins/Kohn_Sham/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Integrals_Bielec MOGuess Bitmask DFT_Utils diff --git a/plugins/Kohn_Sham/damping_SCF.irp.f b/plugins/Kohn_Sham/damping_SCF.irp.f new file mode 100644 index 00000000..aa6f02b0 --- /dev/null +++ b/plugins/Kohn_Sham/damping_SCF.irp.f @@ -0,0 +1,132 @@ +subroutine damping_SCF + implicit none + double precision :: E + double precision, allocatable :: D_alpha(:,:), D_beta(:,:) + double precision :: E_new + double precision, allocatable :: D_new_alpha(:,:), D_new_beta(:,:), F_new(:,:) + double precision, allocatable :: delta_alpha(:,:), delta_beta(:,:) + double precision :: lambda, E_half, a, b, delta_D, delta_E, E_min + + integer :: i,j,k + logical :: saving + character :: save_char + + allocate( & + D_alpha( ao_num_align, ao_num ), & + D_beta( ao_num_align, ao_num ), & + F_new( ao_num_align, ao_num ), & + D_new_alpha( ao_num_align, ao_num ), & + D_new_beta( ao_num_align, ao_num ), & + delta_alpha( ao_num_align, ao_num ), & + delta_beta( ao_num_align, ao_num )) + + do j=1,ao_num + do i=1,ao_num + D_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) + D_beta (i,j) = HF_density_matrix_ao_beta (i,j) + enddo + enddo + + + call write_time(output_hartree_fock) + + write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + '====','================','================','================', '====' + write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + ' N ', 'Energy ', 'Energy diff ', 'Density diff ', 'Save' + write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') & + '====','================','================','================', '====' + + E = HF_energy + 1.d0 + E_min = HF_energy + delta_D = 0.d0 + do k=1,n_it_scf_max + + delta_E = HF_energy - E + E = HF_energy + + if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then + exit + endif + + saving = E < E_min + if (saving) then + call save_mos + save_char = 'X' + E_min = E + else + save_char = ' ' + endif + + write(output_hartree_fock,'(I4,1X,F16.10, 1X, F16.10, 1X, F16.10, 3X, A )') & + k, E, delta_E, delta_D, save_char + + D_alpha = HF_density_matrix_ao_alpha + D_beta = HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + + D_new_alpha = HF_density_matrix_ao_alpha + D_new_beta = HF_density_matrix_ao_beta + F_new = Fock_matrix_ao + E_new = HF_energy + + delta_alpha = D_new_alpha - D_alpha + delta_beta = D_new_beta - D_beta + + lambda = .5d0 + E_half = 0.d0 + do while (E_half > E) + HF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha + HF_density_matrix_ao_beta = D_beta + lambda * delta_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + E_half = HF_energy + if ((E_half > E).and.(E_new < E)) then + lambda = 1.d0 + exit + else if ((E_half > E).and.(lambda > 5.d-4)) then + lambda = 0.5d0 * lambda + E_new = E_half + else + exit + endif + enddo + + a = (E_new + E - 2.d0*E_half)*2.d0 + b = -E_new - 3.d0*E + 4.d0*E_half + lambda = -lambda*b/(a+1.d-16) + D_alpha = (1.d0-lambda) * D_alpha + lambda * D_new_alpha + D_beta = (1.d0-lambda) * D_beta + lambda * D_new_beta + delta_E = HF_energy - E + do j=1,ao_num + do i=1,ao_num + delta_D = delta_D + & + (D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j))*(D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j)) + & + (D_beta (i,j) - HF_density_matrix_ao_beta (i,j))*(D_beta (i,j) - HF_density_matrix_ao_beta (i,j)) + enddo + enddo + delta_D = dsqrt(delta_D/dble(ao_num)**2) + HF_density_matrix_ao_alpha = D_alpha + HF_density_matrix_ao_beta = D_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + + + enddo + write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') '====','================','================','================', '====' + write(output_hartree_fock,*) + + if(.not.no_oa_or_av_opt)then + call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1) + endif + + call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy') + call ezfio_set_hartree_fock_energy(E_min) + + call write_time(output_hartree_fock) + + deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta) +end diff --git a/plugins/Kohn_Sham/diagonalize_fock.irp.f b/plugins/Kohn_Sham/diagonalize_fock.irp.f new file mode 100644 index 00000000..c80077b3 --- /dev/null +++ b/plugins/Kohn_Sham/diagonalize_fock.irp.f @@ -0,0 +1,119 @@ + BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (ao_num) ] +&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Diagonal Fock matrix in the MO basis + END_DOC + + integer :: i,j + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), S(:,:) + + + allocate( F(mo_tot_num_align,mo_tot_num) ) + do j=1,mo_tot_num + do i=1,mo_tot_num + F(i,j) = Fock_matrix_mo(i,j) + enddo + enddo + if(no_oa_or_av_opt)then + integer :: iorb,jorb + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + F(iorb,jorb) = 0.d0 + F(jorb,iorb) = 0.d0 + enddo + enddo + endif + + + + + ! Insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F(i,i) += 0.5d0*level_shift + enddo + + do i = elec_alpha_num+1, mo_tot_num + F(i,i) += level_shift + enddo + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork), iwork(liwork) ) + + lwork = -1 + liwork = -1 + + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) + + call dsyevd( 'V', 'U', mo_tot_num, F, & + size(F,1), diagonal_Fock_matrix_mo, & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed : ', info + stop 1 + endif + + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), F, size(F,1), & + 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) + deallocate(work, iwork, F) + + +! endif + +END_PROVIDER + +BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] + implicit none + BEGIN_DOC + ! diagonal element of the fock matrix calculated as the sum over all the interactions + ! with all the electrons in the RHF determinant + ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij + END_DOC + integer :: i,j + double precision :: accu + do j = 1,elec_alpha_num + accu = 0.d0 + do i = 1, elec_alpha_num + 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(j) = accu + mo_mono_elec_integral(j,j) + enddo + do j = elec_alpha_num+1,mo_tot_num + accu = 0.d0 + do i = 1, elec_alpha_num + 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(j) = accu + mo_mono_elec_integral(j,j) + enddo + +END_PROVIDER diff --git a/plugins/Kohn_Sham/potential_functional.irp.f b/plugins/Kohn_Sham/potential_functional.irp.f new file mode 100644 index 00000000..3502581b --- /dev/null +++ b/plugins/Kohn_Sham/potential_functional.irp.f @@ -0,0 +1,31 @@ + BEGIN_PROVIDER [double precision, ao_potential_alpha_xc, (ao_num_align, ao_num)] +&BEGIN_PROVIDER [double precision, ao_potential_beta_xc, (ao_num_align, ao_num)] + implicit none + integer :: i,j,k,l + ao_potential_alpha_xc = 0.d0 + ao_potential_beta_xc = 0.d0 +!if(exchange_functional == "LDA")then + do i = 1, ao_num + do j = 1, ao_num + ao_potential_alpha_xc(i,j) = (1.d0 - HF_exchange) * lda_ex_potential_alpha_ao(i,j,1) + ao_potential_beta_xc(i,j) = (1.d0 - HF_exchange) * lda_ex_potential_beta_ao(i,j,1) + enddo + enddo +!endif +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_exchange_dft] + implicit none +!if(exchange_functional == "LDA")then + e_exchange_dft = lda_exchange(1) +!endif + +END_PROVIDER + +BEGIN_PROVIDER [double precision, e_correlation_dft] + implicit none +!if(correlation_functional == "LDA")then + e_correlation_dft = 0.d0 +!endif + +END_PROVIDER diff --git a/src/AO_Basis/aos_value.irp.f b/src/AO_Basis/aos_value.irp.f index a531ce50..4876844c 100644 --- a/src/AO_Basis/aos_value.irp.f +++ b/src/AO_Basis/aos_value.irp.f @@ -26,6 +26,7 @@ double precision function ao_value(i,r) do m=1,ao_prim_num(i) beta = ao_expo_ordered_transp(m,i) accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) +! accu += ao_coef_transp(m,i) * dexp(-beta*r2) enddo ao_value = accu * dx * dy * dz diff --git a/src/Utils/constants.include.F b/src/Utils/constants.include.F index c077eb53..4655a4fc 100644 --- a/src/Utils/constants.include.F +++ b/src/Utils/constants.include.F @@ -14,3 +14,4 @@ double precision, parameter :: cx_lda = -0.73855876638202234d0 double precision, parameter :: c_2_4_3 = 2.5198420997897464d0 double precision, parameter :: cst_lda = -0.93052573634909996d0 double precision, parameter :: c_4_3 = 1.3333333333333333d0 +double precision, parameter :: c_1_3 = 0.3333333333333333d0