From 3b75503bc189527aaaef0f0d03934ec1ef544978 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 7 Aug 2020 11:47:51 +0200 Subject: [PATCH 01/26] minor changes --- src/determinants/psi_energy_mono_elec.irp.f | 30 +++++++++++++++++++++ src/utils/integration.irp.f | 10 ++++++- 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/src/determinants/psi_energy_mono_elec.irp.f b/src/determinants/psi_energy_mono_elec.irp.f index 74e69160..fc0ffba9 100644 --- a/src/determinants/psi_energy_mono_elec.irp.f +++ b/src/determinants/psi_energy_mono_elec.irp.f @@ -27,3 +27,33 @@ psi_energy_h_core(i) = psi_energy_h_core(i) * accu enddo END_PROVIDER + + BEGIN_PROVIDER [ double precision, v_ne_psi_energy, (N_states) ] + implicit none + integer :: i + integer :: j,k + double precision :: tmp(mo_num,mo_num),mono_ints(mo_num,mo_num) + BEGIN_DOC +! v_ne_psi_energy = $\langle \Psi | v_ne |\Psi \rangle$ +! +! computed using the :c:data:`one_e_dm_mo_alpha` + +! :c:data:`one_e_dm_mo_beta` and :c:data:`mo_one_e_integrals` + END_DOC + v_ne_psi_energy = 0.d0 + do i = 1, N_states + do j = 1, mo_num + do k = 1, mo_num + v_ne_psi_energy(i) += mo_integrals_n_e(k,j) * (one_e_dm_mo_alpha(k,j,i) + one_e_dm_mo_beta(k,j,i)) + enddo + enddo + enddo + double precision :: accu + do i = 1, N_states + accu = 0.d0 + do j = 1, mo_num + accu += one_e_dm_mo_alpha(j,j,i) + one_e_dm_mo_beta(j,j,i) + enddo + accu = (elec_alpha_num + elec_beta_num ) / accu + v_ne_psi_energy(i) = v_ne_psi_energy(i) * accu + enddo +END_PROVIDER diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index c907e425..da8120a1 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -30,7 +30,11 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph ab = alpha * beta d_AB = (A_center - B_center) * (A_center - B_center) P_center = (alpha * A_center + beta * B_center) * p_inv - fact_k = exp(-ab*p_inv * d_AB) + if(dabs(ab*p_inv * d_AB).lt.50.d0)then + fact_k = exp(-ab*p_inv * d_AB) + else + fact_k = 0.d0 + endif ! Recenter the polynomials P_a and P_b on x !DIR$ FORCEINLINE @@ -78,6 +82,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, !DIR$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) if (fact_k < thresh) then + P_center = 0.d0 + p = 1.d-10 + P_new = 0.d0 + iorder = -1 fact_k = 0.d0 return endif From 03445e1a6ec05e02b8c5268b2db67e3e3c6a447a Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 11 Sep 2020 12:53:06 +0200 Subject: [PATCH 02/26] added hcore_guess in tools --- src/tools/hcore_guess.irp.f | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 src/tools/hcore_guess.irp.f diff --git a/src/tools/hcore_guess.irp.f b/src/tools/hcore_guess.irp.f new file mode 100644 index 00000000..87d0cb7d --- /dev/null +++ b/src/tools/hcore_guess.irp.f @@ -0,0 +1,3 @@ +program hcore_guess_prog + call hcore_guess +end From ee267e27e9b6027ea76a05e2b43271cc8906fd55 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 21 Sep 2020 15:38:26 +0200 Subject: [PATCH 03/26] minor modifs --- config/gfortran_debug.cfg | 3 +- src/basis_correction/print_routine.irp.f | 1 + src/dft_utils_in_r/mo_in_r.irp.f | 3 +- src/mo_basis/mos_in_r.irp.f | 10 ++-- src/mo_guess/h_core_guess_routine.irp.f | 4 +- src/mu_of_r/mu_of_r_conditions.irp.f | 67 ++++++++++++++++++++++++ src/utils/util.irp.f | 12 +++++ 7 files changed, 91 insertions(+), 9 deletions(-) diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 342acae9..926255e0 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -51,7 +51,8 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [DEBUG] -FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan +#FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan +FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan # OpenMP flags ################# diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 05fbbf60..3d497552 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -75,6 +75,7 @@ subroutine print_basis_correction print*,'**************' do istate = 1, N_states write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate) + write(*, '(A29,X,I3,X,A3,X,F16.10)') 'mu_average_trans_corr, state ',istate,' = ',mu_average_trans_corr(istate) enddo end diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 81863c3a..fa3ec73c 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -79,7 +79,7 @@ END_DOC integer :: m integer :: i,j - mos_grad_in_r_array = 0.d0 + mos_grad_in_r_array_tranp = 0.d0 do i = 1, n_points_final_grid do j = 1, mo_num do m = 1, 3 @@ -126,4 +126,3 @@ enddo END_PROVIDER - diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index 049db8aa..ffcd5a9a 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -12,18 +12,18 @@ subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array) implicit none double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) - double precision, intent(out) :: mos_grad_array(mo_num,3) + double precision, intent(out) :: mos_grad_array(3,mo_num) integer :: i,j,k - double precision :: aos_array(ao_num),aos_grad_array(ao_num,3) + double precision :: aos_array(ao_num),aos_grad_array(3,ao_num) call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) mos_array=0d0 mos_grad_array=0d0 do j = 1, mo_num do k=1, ao_num mos_array(j) += mo_coef(k,j)*aos_array(k) - mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1) - mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2) - mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3) + mos_grad_array(1,j) += mo_coef(k,j)*aos_grad_array(1,k) + mos_grad_array(2,j) += mo_coef(k,j)*aos_grad_array(2,k) + mos_grad_array(3,j) += mo_coef(k,j)*aos_grad_array(3,k) enddo enddo end diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index 8fc3f6f2..9c5e76ff 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -9,5 +9,7 @@ subroutine hcore_guess size(mo_one_e_integrals,1), & size(mo_one_e_integrals,2),label,1,.false.) call save_mos - SOFT_TOUCH mo_coef mo_label +! SOFT_TOUCH mo_coef mo_label + TOUCH mo_coef mo_label + print*,'mo_one_e_integrals(1,1) = ',mo_one_e_integrals(1,1) end diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 05a2a4d7..b35a2a8a 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -28,6 +28,8 @@ mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + else if(mu_of_r_potential.EQ."transcorr")then + mu_of_r_prov(ipoint,istate) = mu_of_r_transcorr(ipoint,istate) else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential @@ -124,6 +126,47 @@ print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0 END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_transcorr, (n_points_final_grid,N_states) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a wave function developped in an active space + ! + ! corresponds to \sqrt(pi) * (W(0) + 1/4)/3 + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + ! + ! + END_DOC + integer :: ipoint,istate + double precision :: wall0,wall1,f_psi,on_top,w_psi,sqpi + print*,'providing mu_of_r_transcorr ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + + provide f_psi_cas_ab + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & + !$OMP SHARED (n_points_final_grid,mu_of_r_transcorr,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + f_psi = f_psi_cas_ab(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) + if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then + w_psi = 1.d+10 + else + w_psi = f_psi / on_top + endif + mu_of_r_transcorr(ipoint,istate) = (0.25d0 + w_psi) * sqpi / 3.d0 + enddo + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_transcorr = ',wall1-wall0 + END_PROVIDER + BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] implicit none @@ -148,3 +191,27 @@ mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) enddo END_PROVIDER + + BEGIN_PROVIDER [double precision, mu_average_trans_corr, (N_states)] + implicit none + BEGIN_DOC + ! average value of mu(r) weighted with the total one-e density and divised by the number of electrons + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: weight,density + mu_average_trans_corr = 0.d0 + do istate = 1, N_states + do ipoint = 1, n_points_final_grid + weight =final_weight_at_r_vector(ipoint) + density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) + if(mu_of_r_transcorr(ipoint,istate).gt.1.d+09)cycle + mu_average_trans_corr(istate) += mu_of_r_transcorr(ipoint,istate) * weight * density + enddo + mu_average_trans_corr(istate) = mu_average_trans_corr(istate) / elec_num_grid_becke(istate) + enddo + END_PROVIDER diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index 1b01a1ec..cfb42fd1 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -1,3 +1,15 @@ +double precision function derf_mu_x(mu,x) + implicit none + include 'utils/constants.include.F' + double precision, intent(in) :: mu,x + if(dabs(x).gt.1.d-6)then + derf_mu_x = derf(mu * x)/x + else + derf_mu_x = inv_sq_pi * 2.d0 * mu * (1.d0 - mu*mu*x*x/3.d0) + endif +end + + double precision function binom_func(i,j) implicit none BEGIN_DOC From 01b1ee327392efbab00af9ff0badef73f9f91ca9 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 7 Oct 2020 11:03:23 +0200 Subject: [PATCH 04/26] fixed bug in laplacians --- src/dft_utils_in_r/ao_in_r.irp.f | 39 -------------------------------- src/dft_utils_in_r/mo_in_r.irp.f | 21 +++++++++++++++-- src/mo_basis/mos_in_r.irp.f | 16 ++++++------- 3 files changed, 27 insertions(+), 49 deletions(-) diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 4b1526dd..2892d2e2 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -121,42 +121,3 @@ END_PROVIDER - BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)] - implicit none - ! - ! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point - ! - ! k = 1 : x, k= 2, y, k 3, z - integer :: i,j,m - do m = 1, 3 - do i = 1, n_points_final_grid - do j = 1, ao_num - aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_in_r_array(j,i,m) - enddo - enddo - enddo - END_PROVIDER - - - BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)] -&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)] - implicit none - BEGIN_DOC - ! aos_in_r_array_per_atom(i,j,k) = value of the ith ao on the jth grid point attached on the kth atom - END_DOC - integer :: i,j,k - double precision :: aos_array(ao_num), r(3) - do k = 1, nucl_num - do i = 1, n_pts_per_atom(k) - r(1) = final_grid_points_per_atom(1,i,k) - r(2) = final_grid_points_per_atom(2,i,k) - r(3) = final_grid_points_per_atom(3,i,k) - call give_all_aos_at_r(r,aos_array) - do j = 1, ao_num - aos_in_r_array_per_atom(j,i,k) = aos_array(j) - aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j) - enddo - enddo - enddo - END_PROVIDER - diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index fa3ec73c..d4afcf76 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -115,8 +115,6 @@ BEGIN_DOC ! mos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith mo on the jth grid point ! - ! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth mo on the ith grid point - ! ! k = 1 : x, k= 2, y, k 3, z END_DOC integer :: m @@ -126,3 +124,22 @@ enddo END_PROVIDER + + BEGIN_PROVIDER[double precision, mos_lapl_in_r_array_tranp,(3,mo_num,n_points_final_grid)] + implicit none + BEGIN_DOC + ! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplient of jth mo on the ith grid point + ! + ! k = 1 : x, k= 2, y, k 3, z + END_DOC + integer :: m + integer :: i,j + mos_lapl_in_r_array_tranp = 0.d0 + do i = 1, n_points_final_grid + do j = 1, mo_num + do m = 1, 3 + mos_lapl_in_r_array_tranp(m,j,i) = mos_lapl_in_r_array(j,i,m) + enddo + enddo + enddo + END_PROVIDER diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index ffcd5a9a..ee2795d0 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -33,9 +33,9 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la implicit none double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) - double precision, intent(out) :: mos_grad_array(mo_num,3),mos_lapl_array(mo_num,3) + double precision, intent(out) :: mos_grad_array(3,mo_num),mos_lapl_array(3,mo_num) integer :: i,j,k - double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3) + double precision :: aos_array(ao_num),aos_grad_array(3,ao_num),aos_lapl_array(3,ao_num) call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) mos_array = 0.d0 mos_grad_array = 0.d0 @@ -43,12 +43,12 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la do j = 1, mo_num do k=1, ao_num mos_array(j) += mo_coef(k,j) * aos_array(k) - mos_grad_array(j,1) += mo_coef(k,j) * aos_grad_array(k,1) - mos_grad_array(j,2) += mo_coef(k,j) * aos_grad_array(k,2) - mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3) - mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1) - mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2) - mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3) + mos_grad_array(1,j) += mo_coef(k,j) * aos_grad_array(1,k) + mos_grad_array(2,j) += mo_coef(k,j) * aos_grad_array(2,k) + mos_grad_array(3,j) += mo_coef(k,j) * aos_grad_array(3,k) + mos_lapl_array(1,j) += mo_coef(k,j) * aos_lapl_array(1,k) + mos_lapl_array(2,j) += mo_coef(k,j) * aos_lapl_array(2,k) + mos_lapl_array(3,j) += mo_coef(k,j) * aos_lapl_array(3,k) enddo enddo end From c545098af083a423e2b4f6ce7597c1d9ce566c38 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 19 Oct 2020 12:20:14 +0200 Subject: [PATCH 05/26] added tools/huckel_guess.irp.f --- src/tools/huckel_guess.irp.f | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 src/tools/huckel_guess.irp.f diff --git a/src/tools/huckel_guess.irp.f b/src/tools/huckel_guess.irp.f new file mode 100644 index 00000000..5ec37df4 --- /dev/null +++ b/src/tools/huckel_guess.irp.f @@ -0,0 +1,5 @@ +program pouet + implicit none + call huckel_guess + +end From 2c33bca408aaca9003d2af93c8ba63721cefcb88 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 28 Oct 2020 11:24:32 +0100 Subject: [PATCH 06/26] minor changes --- src/dft_utils_in_r/ao_in_r.irp.f | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 2892d2e2..6abda550 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -120,4 +120,18 @@ !$OMP END PARALLEL DO END_PROVIDER + BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)] + implicit none + integer :: i,j,m + double precision :: aos_array(ao_num), r(3) + double precision :: aos_grad_array(3,ao_num) + do m = 1, 3 + do j = 1, ao_num + do i = 1, n_points_final_grid + aos_grad_in_r_array_transp_bis(i,j,m) = aos_grad_in_r_array(j,i,m) + enddo + enddo + enddo + END_PROVIDER + From bd8fcdb03030af73d5f4bef74a6ccd0a8899702d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 2 Nov 2020 17:24:35 +0100 Subject: [PATCH 07/26] removed stupid bug in utils/one_e_integration.irp.f --- src/ao_one_e_ints/ao_overlap.irp.f | 7 +++++++ src/utils/one_e_integration.irp.f | 22 +++++++++++----------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index 11c95e42..d9061d67 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -54,6 +54,13 @@ call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) ao_overlap(i,j) += c * overlap + if(isnan(ao_overlap(i,j)))then + print*,'i,j',i,j + print*,'l,n',l,n + print*,'c,overlap',c,overlap + print*,overlap_x,overlap_y,overlap_z + stop + endif ao_overlap_x(i,j) += c * overlap_x ao_overlap_y(i,j) += c * overlap_y ao_overlap_z(i,j) += c * overlap_z diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index 97eef89d..e7b3e051 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_ call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,& beta,power_A,power_B,A_center,B_center,dim) -! if(fact_p.lt.0.000001d0)then -! overlap_gaussian_x = 0.d0 -! return -! endif + if(fact_p.lt.1.d-20)then + overlap_gaussian_x = 0.d0 + return + endif overlap_gaussian_x = 0.d0 integer :: i @@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& integer :: iorder_p(3) call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) -! if(fact_p.lt.1d-20)then -! overlap_x = 0.d0 -! overlap_y = 0.d0 -! overlap_z = 0.d0 -! overlap = 0.d0 -! return -! endif + if(fact_p.lt.1d-20)then + overlap_x = 0.d0 + overlap_y = 0.d0 + overlap_z = 0.d0 + overlap = 0.d0 + return + endif integer :: nmax double precision :: F_integral nmax = maxval(iorder_p) From bdbb08207956051bb427051f184449985628fd23 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sun, 8 Nov 2020 17:53:37 +0100 Subject: [PATCH 08/26] fixed a little bun in src/utils/one_e_integration.irp.f --- src/utils/one_e_integration.irp.f | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f index 97eef89d..2c597104 100644 --- a/src/utils/one_e_integration.irp.f +++ b/src/utils/one_e_integration.irp.f @@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& integer :: iorder_p(3) call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) -! if(fact_p.lt.1d-20)then -! overlap_x = 0.d0 -! overlap_y = 0.d0 -! overlap_z = 0.d0 -! overlap = 0.d0 -! return -! endif + if(fact_p.lt.1d-20)then + overlap_x = 1.d-10 + overlap_y = 1.d-10 + overlap_z = 1.d-10 + overlap = 1.d-10 + return + endif integer :: nmax double precision :: F_integral nmax = maxval(iorder_p) From 6652da23f57dbb84f490913346ad31c13e628f68 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sun, 8 Nov 2020 18:26:54 +0100 Subject: [PATCH 09/26] removed symmetry bug in a-b-2-RDM : warning it has no impacts on any quantity as long as all the elements of the two-rdm were used --- src/two_rdm_routines/update_rdm.irp.f | 20 +++++++++++++++++-- .../update_state_av_rdm.irp.f | 16 +++++++++++++-- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/src/two_rdm_routines/update_rdm.irp.f b/src/two_rdm_routines/update_rdm.irp.f index 4d74280e..662b3ee9 100644 --- a/src/two_rdm_routines/update_rdm.irp.f +++ b/src/two_rdm_routines/update_rdm.irp.f @@ -356,12 +356,20 @@ h2 = list_orb_reverse(h2) nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) * phase + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo else ! Mono beta @@ -377,12 +385,20 @@ h2 = list_orb_reverse(h2) nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) * phase + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = h2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo endif else if(spin_trace)then diff --git a/src/two_rdm_routines/update_state_av_rdm.irp.f b/src/two_rdm_routines/update_state_av_rdm.irp.f index 35024331..bf7b6148 100644 --- a/src/two_rdm_routines/update_state_av_rdm.irp.f +++ b/src/two_rdm_routines/update_state_av_rdm.irp.f @@ -327,11 +327,17 @@ if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 - values(nkeys) = c_1 * phase + values(nkeys) = 0.5d0 * c_1 * phase keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo else ! Mono beta @@ -346,11 +352,17 @@ if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 - values(nkeys) = c_1 * phase + values(nkeys) = 0.5d0 * c_1 * phase keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = h2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 enddo endif else if(spin_trace)then From aa0c44959c69d0355e35c09d2798b3c053122305 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 9 Nov 2020 18:30:10 +0100 Subject: [PATCH 10/26] removed stupid dependency in basis_set_correction --- src/basis_correction/print_routine.irp.f | 1 - 1 file changed, 1 deletion(-) diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 3d497552..05fbbf60 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -75,7 +75,6 @@ subroutine print_basis_correction print*,'**************' do istate = 1, N_states write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate) - write(*, '(A29,X,I3,X,A3,X,F16.10)') 'mu_average_trans_corr, state ',istate,' = ',mu_average_trans_corr(istate) enddo end From 5b8580fe2dc993a469f3376d4e01ee7b31b7a552 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 2 Jan 2021 15:40:03 +0100 Subject: [PATCH 11/26] some cleaning --- src/dft_utils_in_r/ao_in_r.irp.f | 4 ++ src/mo_guess/h_core_guess_routine.irp.f | 2 - src/mu_of_r/mu_of_r_conditions.irp.f | 66 ------------------------- 3 files changed, 4 insertions(+), 68 deletions(-) diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 6abda550..a98afd5d 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -122,6 +122,10 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)] implicit none + BEGIN_DOC +! Transposed gradients +! + END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) double precision :: aos_grad_array(3,ao_num) diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index 9c5e76ff..246dfef2 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -9,7 +9,5 @@ subroutine hcore_guess size(mo_one_e_integrals,1), & size(mo_one_e_integrals,2),label,1,.false.) call save_mos -! SOFT_TOUCH mo_coef mo_label TOUCH mo_coef mo_label - print*,'mo_one_e_integrals(1,1) = ',mo_one_e_integrals(1,1) end diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index b35a2a8a..148c65b3 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -28,8 +28,6 @@ mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) - else if(mu_of_r_potential.EQ."transcorr")then - mu_of_r_prov(ipoint,istate) = mu_of_r_transcorr(ipoint,istate) else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential @@ -126,47 +124,6 @@ print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0 END_PROVIDER - BEGIN_PROVIDER [double precision, mu_of_r_transcorr, (n_points_final_grid,N_states) ] - implicit none - BEGIN_DOC - ! mu(r) computed with a wave function developped in an active space - ! - ! corresponds to \sqrt(pi) * (W(0) + 1/4)/3 - ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals - ! - ! in the one- and two-body density matrix are excluded - ! - ! - END_DOC - integer :: ipoint,istate - double precision :: wall0,wall1,f_psi,on_top,w_psi,sqpi - print*,'providing mu_of_r_transcorr ...' - call wall_time(wall0) - sqpi = dsqrt(dacos(-1.d0)) - - provide f_psi_cas_ab - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & - !$OMP SHARED (n_points_final_grid,mu_of_r_transcorr,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) - do istate = 1, N_states - do ipoint = 1, n_points_final_grid - f_psi = f_psi_cas_ab(ipoint,istate) - on_top = on_top_cas_mu_r(ipoint,istate) - if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then - w_psi = 1.d+10 - else - w_psi = f_psi / on_top - endif - mu_of_r_transcorr(ipoint,istate) = (0.25d0 + w_psi) * sqpi / 3.d0 - enddo - enddo - !$OMP END PARALLEL DO - call wall_time(wall1) - print*,'Time to provide mu_of_r_transcorr = ',wall1-wall0 - END_PROVIDER - BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] implicit none @@ -192,26 +149,3 @@ enddo END_PROVIDER - BEGIN_PROVIDER [double precision, mu_average_trans_corr, (N_states)] - implicit none - BEGIN_DOC - ! average value of mu(r) weighted with the total one-e density and divised by the number of electrons - ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals - ! - ! in the one- and two-body density matrix are excluded - END_DOC - integer :: ipoint,istate - double precision :: weight,density - mu_average_trans_corr = 0.d0 - do istate = 1, N_states - do ipoint = 1, n_points_final_grid - weight =final_weight_at_r_vector(ipoint) - density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & - + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) - if(mu_of_r_transcorr(ipoint,istate).gt.1.d+09)cycle - mu_average_trans_corr(istate) += mu_of_r_transcorr(ipoint,istate) * weight * density - enddo - mu_average_trans_corr(istate) = mu_average_trans_corr(istate) / elec_num_grid_becke(istate) - enddo - END_PROVIDER From 970f846a4d4615d7cb73853f9c8b7d2b1e2802c9 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 7 Jan 2021 14:20:11 +0100 Subject: [PATCH 12/26] changed the read/write two_rdm --- src/two_body_rdm/act_2_rdm.irp.f | 24 ++++++++++++++++-------- src/two_body_rdm/io_two_rdm.irp.f | 29 +++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 8 deletions(-) create mode 100644 src/two_body_rdm/io_two_rdm.irp.f diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index e3265572..41b28aea 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -22,6 +22,8 @@ END_DOC integer :: ispin double precision :: wall_1, wall_2 + character*(128) :: name_file + name_file = 'act_2_rdm_ab_mo' ! condition for alpha/beta spin print*,'' print*,'Providing act_2_rdm_ab_mo ' @@ -31,13 +33,13 @@ call wall_time(wall_1) if(read_two_body_rdm_ab)then print*,'Reading act_2_rdm_ab_mo from disk ...' - call ezfio_get_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo) + call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_ab_mo,name_file) else call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_ab)then print*,'Writing act_2_rdm_ab_mo on disk ...' - call ezfio_set_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo) + call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_ab_mo,name_file) call ezfio_set_two_body_rdm_io_two_body_rdm_ab("Read") endif call wall_time(wall_2) @@ -63,18 +65,20 @@ ! condition for alpha/beta spin print*,'' print*,'Providing act_2_rdm_aa_mo ' + character*(128) :: name_file + name_file = 'act_2_rdm_aa_mo' ispin = 1 act_2_rdm_aa_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_aa)then print*,'Reading act_2_rdm_aa_mo from disk ...' - call ezfio_get_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo) + call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_aa_mo,name_file) else call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_aa)then print*,'Writing act_2_rdm_aa_mo on disk ...' - call ezfio_set_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo) + call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_aa_mo,name_file) call ezfio_set_two_body_rdm_io_two_body_rdm_aa("Read") endif @@ -101,18 +105,20 @@ ! condition for beta/beta spin print*,'' print*,'Providing act_2_rdm_bb_mo ' + character*(128) :: name_file + name_file = 'act_2_rdm_bb_mo' ispin = 2 act_2_rdm_bb_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_bb)then print*,'Reading act_2_rdm_bb_mo from disk ...' - call ezfio_get_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo) + call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_bb_mo,name_file) else call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_bb)then print*,'Writing act_2_rdm_bb_mo on disk ...' - call ezfio_set_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo) + call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_bb_mo,name_file) call ezfio_set_two_body_rdm_io_two_body_rdm_bb("Read") endif @@ -138,18 +144,20 @@ ! condition for beta/beta spin print*,'' print*,'Providing act_2_rdm_spin_trace_mo ' + character*(128) :: name_file + name_file = 'act_2_rdm_spin_trace_mo' ispin = 4 act_2_rdm_spin_trace_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_spin_trace)then print*,'Reading act_2_rdm_spin_trace_mo from disk ...' - call ezfio_get_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo) + call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_spin_trace_mo,name_file) else call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_spin_trace)then print*,'Writing act_2_rdm_spin_trace_mo on disk ...' - call ezfio_set_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo) + call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_spin_trace_mo,name_file) call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read") endif diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f new file mode 100644 index 00000000..f7008ca9 --- /dev/null +++ b/src/two_body_rdm/io_two_rdm.irp.f @@ -0,0 +1,29 @@ +subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates) + + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'W') + write(i_unit_output)array_tmp + close(unit=i_unit_output) +end + +subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file) + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + integer, intent(in) :: n_orb,nstates + character*(128), intent(in) :: name_file + double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states) + PROVIDE ezfio_filename + output=trim(ezfio_filename)//'/work/'//trim(name_file) + i_unit_output = getUnitAndOpen(output,'R') + read(i_unit_output)array_tmp + close(unit=i_unit_output) +end + From 7874e4f0373addc220e153d4fee362469629f850 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 7 Jan 2021 14:24:44 +0100 Subject: [PATCH 13/26] removed useless arrays for two_rdm --- src/two_body_rdm/EZFIO.cfg | 24 ------------------------ 1 file changed, 24 deletions(-) diff --git a/src/two_body_rdm/EZFIO.cfg b/src/two_body_rdm/EZFIO.cfg index 4ca39d73..d73f704d 100644 --- a/src/two_body_rdm/EZFIO.cfg +++ b/src/two_body_rdm/EZFIO.cfg @@ -1,45 +1,21 @@ -[two_rdm_ab_disk] -type: double precision -doc: active part of the two body rdm alpha/beta stored on disk -interface: ezfio -size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states) - [io_two_body_rdm_ab] type: Disk_access doc: Read/Write the active part of the two-body rdm for alpha/beta electrons from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[two_rdm_aa_disk] -type: double precision -doc: active part of the two body rdm alpha/alpha stored on disk -interface: ezfio -size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states) - [io_two_body_rdm_aa] type: Disk_access doc: Read/Write the active part of the two-body rdm for alpha/alpha electrons from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[two_rdm_bb_disk] -type: double precision -doc: active part of the two body rdm beta/beta stored on disk -interface: ezfio -size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states) - [io_two_body_rdm_bb] type: Disk_access doc: Read/Write the active part of the two-body rdm for beta/beta electrons from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[two_rdm_spin_trace_disk] -type: double precision -doc: active part of the two body rdm spin trace stored on disk -interface: ezfio -size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states) - [io_two_body_rdm_spin_trace] type: Disk_access doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ] From e5da77577fdaa29f2990256d3a245e035f2feb60 Mon Sep 17 00:00:00 2001 From: bpradines Date: Thu, 14 Jan 2021 11:05:03 +0100 Subject: [PATCH 14/26] removed stupid dependency on the AO integrals in basis_correction --- src/basis_correction/basis_correction.irp.f | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/src/basis_correction/basis_correction.irp.f b/src/basis_correction/basis_correction.irp.f index 8d3264bc..a7ea7244 100644 --- a/src/basis_correction/basis_correction.irp.f +++ b/src/basis_correction/basis_correction.irp.f @@ -7,22 +7,10 @@ program basis_correction touch read_wf no_core_density = .True. touch no_core_density - provide ao_two_e_integrals_in_map + if(io_mo_two_e_integrals .ne. "Read")then + provide ao_two_e_integrals_in_map + endif provide mo_two_e_integrals_in_map call print_basis_correction -! call print_e_b end -subroutine print_e_b - implicit none - print *, 'Hello world' - print*,'ecmd_lda_mu_of_r = ',ecmd_lda_mu_of_r - print*,'ecmd_pbe_ueg_mu_of_r = ',ecmd_pbe_ueg_mu_of_r - print*,'ecmd_pbe_ueg_eff_xi_mu_of_r = ',ecmd_pbe_ueg_eff_xi_mu_of_r - print*,'' - print*,'psi_energy + E^B_LDA = ',psi_energy + ecmd_lda_mu_of_r - print*,'psi_energy + E^B_PBE_UEG = ',psi_energy + ecmd_pbe_ueg_mu_of_r - print*,'psi_energy + E^B_PBE_UEG_Xi = ',psi_energy + ecmd_pbe_ueg_eff_xi_mu_of_r - print*,'' - print*,'mu_average_prov = ',mu_average_prov -end From acd9192aa4393f41d1c81cea20248d8b6d22083b Mon Sep 17 00:00:00 2001 From: bpradines Date: Tue, 19 Jan 2021 19:01:52 +0100 Subject: [PATCH 15/26] fixed bug in weird dependence of mu(r) --- src/mu_of_r/f_hf_utils.irp.f | 1 + src/mu_of_r/f_psi_i_a_v_utils.irp.f | 3 +++ 2 files changed, 4 insertions(+) diff --git a/src/mu_of_r/f_hf_utils.irp.f b/src/mu_of_r/f_hf_utils.irp.f index b89dda18..8480a288 100644 --- a/src/mu_of_r/f_hf_utils.irp.f +++ b/src/mu_of_r/f_hf_utils.irp.f @@ -9,6 +9,7 @@ BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: get_two_e_integral + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 m = list_valence_orb_for_hf(orb_m,1) do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f index aed054ae..427da199 100644 --- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -235,6 +235,7 @@ BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: integrals_array(mo_num,mo_num),get_two_e_integral + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals do orb_m = 1, n_act_orb ! electron 1 m = list_act(orb_m) do orb_n = 1, n_act_orb ! electron 2 @@ -264,6 +265,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_ina END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: integrals_array(mo_num,mo_num),get_two_e_integral + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals do orb_m = 1, n_act_orb ! electron 1 m = list_act(orb_m) do orb_n = 1, n_inact_orb ! electron 2 @@ -293,6 +295,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_ina END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: get_two_e_integral,integrals_array(mo_num,mo_num) + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals do orb_m = 1, n_inact_orb ! electron 1 m = list_inact(orb_m) do orb_n = 1, n_inact_orb ! electron 2 From d6d8b0586c60cd3184e9e419b476ddaf5f7219d1 Mon Sep 17 00:00:00 2001 From: bpradines Date: Fri, 26 Feb 2021 10:18:42 +0100 Subject: [PATCH 16/26] added basis_correction/print_su_pbe_ot.irp.f --- src/basis_correction/print_su_pbe_ot.irp.f | 28 ++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 src/basis_correction/print_su_pbe_ot.irp.f diff --git a/src/basis_correction/print_su_pbe_ot.irp.f b/src/basis_correction/print_su_pbe_ot.irp.f new file mode 100644 index 00000000..49f90ade --- /dev/null +++ b/src/basis_correction/print_su_pbe_ot.irp.f @@ -0,0 +1,28 @@ +program basis_corr_su_pbe_ot + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + read_wf = .True. + touch read_wf + no_core_density = .True. + touch no_core_density + if(io_mo_two_e_integrals .ne. "Read")then + provide ao_two_e_integrals_in_map + endif + provide mo_two_e_integrals_in_map + call print_su_pbe_ot + +end + +subroutine print_su_pbe_ot + implicit none + integer :: istate + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) + enddo + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate) + enddo + +end From 3e551e153742c6effd0680ef97620816a382b5f6 Mon Sep 17 00:00:00 2001 From: bpradines Date: Fri, 26 Feb 2021 10:48:26 +0100 Subject: [PATCH 17/26] added dipole moments in determinants --- src/determinants/dipole_moments.irp.f | 65 +++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 src/determinants/dipole_moments.irp.f diff --git a/src/determinants/dipole_moments.irp.f b/src/determinants/dipole_moments.irp.f new file mode 100644 index 00000000..8a5f1a2d --- /dev/null +++ b/src/determinants/dipole_moments.irp.f @@ -0,0 +1,65 @@ + BEGIN_PROVIDER [double precision, z_dipole_moment, (N_states)] +&BEGIN_PROVIDER [double precision, y_dipole_moment, (N_states)] +&BEGIN_PROVIDER [double precision, x_dipole_moment, (N_states)] + implicit none + BEGIN_DOC + ! blablabla + END_DOC + integer :: ipoint,istate,i,j + double precision :: weight, r(3) + double precision :: cpu0,cpu1,nuclei_part_z,nuclei_part_y,nuclei_part_x + + call cpu_time(cpu0) + z_dipole_moment = 0.d0 + y_dipole_moment = 0.d0 + x_dipole_moment = 0.d0 + do istate = 1, N_states + do i = 1, mo_num + z_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_z(i,i) + y_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_y(i,i) + x_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_x(i,i) + do j = i+1, mo_num + z_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_z(j,i) + y_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_y(j,i) + x_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_x(j,i) + enddo + enddo + enddo + + print*,'electron part for z_dipole = ',z_dipole_moment + print*,'electron part for y_dipole = ',y_dipole_moment + print*,'electron part for x_dipole = ',x_dipole_moment + + nuclei_part_z = 0.d0 + nuclei_part_y = 0.d0 + nuclei_part_x = 0.d0 + do i = 1,nucl_num + nuclei_part_z += nucl_charge(i) * nucl_coord(i,3) + nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) + nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) + enddo + print*,'nuclei part for z_dipole = ',nuclei_part_z + print*,'nuclei part for y_dipole = ',nuclei_part_y + print*,'nuclei part for x_dipole = ',nuclei_part_x + + do istate = 1, N_states + z_dipole_moment(istate) += nuclei_part_z + y_dipole_moment(istate) += nuclei_part_y + x_dipole_moment(istate) += nuclei_part_x + enddo + + call cpu_time(cpu1) + print*,'Time to provide the dipole moment :',cpu1-cpu0 +END_PROVIDER + + + + + subroutine print_z_dipole_moment_only + implicit none + print*, '' + print*, '' + print*, '****************************************' + print*, 'z_dipole_moment = ',z_dipole_moment + print*, '****************************************' + end From c1669181b2e66b2a4f2c24090001ad06d53a2152 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 15 Mar 2021 17:40:37 +0100 Subject: [PATCH 18/26] minor changes --- src/mo_two_e_erf_ints/core_quantities_erf.irp.f | 2 +- src/tools/print_dipole.irp.f | 5 +++++ src/tools/print_var_energy.irp.f | 11 +++++++++++ 3 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 src/tools/print_dipole.irp.f create mode 100644 src/tools/print_var_energy.irp.f diff --git a/src/mo_two_e_erf_ints/core_quantities_erf.irp.f b/src/mo_two_e_erf_ints/core_quantities_erf.irp.f index 3cd68205..43d77904 100644 --- a/src/mo_two_e_erf_ints/core_quantities_erf.irp.f +++ b/src/mo_two_e_erf_ints/core_quantities_erf.irp.f @@ -7,7 +7,7 @@ BEGIN_PROVIDER [double precision, core_energy_erf] core_energy_erf = 0.d0 do i = 1, n_core_orb j = list_core(i) - core_energy_erf += 2.d0 * mo_one_e_integrals(j,j) + mo_two_e_int_erf_jj(j,j) + core_energy_erf += mo_two_e_int_erf_jj(j,j) do k = i+1, n_core_orb l = list_core(k) core_energy_erf += 2.d0 * (2.d0 * mo_two_e_int_erf_jj(j,l) - mo_two_e_int_erf_jj_exchange(j,l)) diff --git a/src/tools/print_dipole.irp.f b/src/tools/print_dipole.irp.f new file mode 100644 index 00000000..8351308e --- /dev/null +++ b/src/tools/print_dipole.irp.f @@ -0,0 +1,5 @@ +program print_dipole + implicit none + call print_z_dipole_moment_only + +end diff --git a/src/tools/print_var_energy.irp.f b/src/tools/print_var_energy.irp.f new file mode 100644 index 00000000..bab57d8c --- /dev/null +++ b/src/tools/print_var_energy.irp.f @@ -0,0 +1,11 @@ +program print_var_energy + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + print*,'psi_energy = ',psi_energy +end From 1935cc845e9a6ad11c8181e5afb45a821a2a76f4 Mon Sep 17 00:00:00 2001 From: bpradines Date: Tue, 30 Mar 2021 14:58:21 +0200 Subject: [PATCH 19/26] removed nan in density_for_dft when no beta electrons --- src/density_for_dft/density_for_dft.irp.f | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/density_for_dft/density_for_dft.irp.f b/src/density_for_dft/density_for_dft.irp.f index ee70cd38..364cedd8 100644 --- a/src/density_for_dft/density_for_dft.irp.f +++ b/src/density_for_dft/density_for_dft.irp.f @@ -41,7 +41,11 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s elec_alpha_valence(istate) += one_e_dm_mo_alpha_for_dft(i,i,istate) enddo elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate) - one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate) + if( dabs(elec_alpha_valence(istate)) .lt.1.d-12)then + one_e_dm_mo_alpha_for_dft = 0.d0 + else + one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate) + endif enddo endif @@ -55,6 +59,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st ! density matrix for beta electrons in the MO basis used for all DFT calculations based on the density END_DOC double precision :: delta_beta(mo_num,mo_num,N_states) + one_e_dm_mo_beta_for_dft = 0.d0 if(density_for_dft .EQ. "damping_rs_dft")then delta_beta = one_e_dm_mo_beta - data_one_e_dm_beta_mo one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta @@ -73,6 +78,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:) endif + if(no_core_density)then integer :: ii,i,j do ii = 1, n_core_orb @@ -82,17 +88,21 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st one_e_dm_mo_beta_for_dft(i,j,:) = 0.d0 enddo enddo + double precision :: elec_beta_valence(N_states),elec_beta_frozen_num + integer :: istate if(normalize_dm)then - double precision :: elec_beta_valence(N_states),elec_beta_frozen_num elec_beta_frozen_num = elec_beta_num - n_core_orb elec_beta_valence = 0.d0 - integer :: istate do istate = 1, N_states do i = 1, mo_num elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate) enddo - elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate) - one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate) + if(dabs(elec_beta_valence(istate)).lt.1.d-12)then + one_e_dm_mo_beta_for_dft = 0.d0 + else + elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate) + one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate) + endif enddo endif endif From 4a6f7a3a928fa634fc1a5384eb43494c59f3d2bf Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 8 Apr 2021 20:37:17 +0200 Subject: [PATCH 20/26] added some stuffs for foboscf --- src/bitmask/modify_bitmasks.irp.f | 8 ++-- src/cas_based_on_top/example.irp.f | 16 ++++--- src/dft_utils_in_r/dm_in_r_routines.irp.f | 56 +++++++++++++++++++++++ src/scf_utils/EZFIO.cfg | 7 +++ src/scf_utils/diagonalize_fock.irp.f | 21 +++++++++ src/scf_utils/fock_matrix.irp.f | 21 +++++++++ 6 files changed, 118 insertions(+), 11 deletions(-) diff --git a/src/bitmask/modify_bitmasks.irp.f b/src/bitmask/modify_bitmasks.irp.f index 834be6c8..ed291adf 100644 --- a/src/bitmask/modify_bitmasks.irp.f +++ b/src/bitmask/modify_bitmasks.irp.f @@ -143,10 +143,10 @@ subroutine print_generators_bitmasks_holes key_tmp(j,1) = generators_bitmask(j,1,i) key_tmp(j,2) = generators_bitmask(j,2,i) enddo - print*,'' - print*,'index hole = ',i - call print_det(key_tmp,N_int) - print*,'' +! print*,'' +! print*,'index hole = ',i +! call print_det(key_tmp,N_int) +! print*,'' enddo deallocate(key_tmp) diff --git a/src/cas_based_on_top/example.irp.f b/src/cas_based_on_top/example.irp.f index f00956e3..2f709495 100644 --- a/src/cas_based_on_top/example.irp.f +++ b/src/cas_based_on_top/example.irp.f @@ -5,37 +5,39 @@ subroutine write_on_top_in_real_space ! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid END_DOC double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist + double precision :: core_dens, inact_dens,act_dens(2,1) integer :: nz,i,istate character*(128) :: output integer :: i_unit_output,getUnitAndOpen PROVIDE ezfio_filename output=trim(ezfio_filename)//'.on_top' - print*,'output = ',trim(output) + print*,'output = ',trim(output) i_unit_output = getUnitAndOpen(output,'w') - zmax = 2.0d0 + zmax = 5.0d0 print*,'nucl_coord(1,3) = ',nucl_coord(1,3) print*,'nucl_coord(2,3) = ',nucl_coord(2,3) dist = dabs(nucl_coord(1,3) - nucl_coord(2,3)) - zmax += dist + zmax += dist zcenter = (nucl_coord(1,3) + nucl_coord(2,3))*0.5d0 print*,'zcenter = ',zcenter print*,'zmax = ',zmax nz = 1000 dz = zmax / dble(nz) - r(:) = 0.d0 - r(3) = zcenter -zmax * 0.5d0 + r(:) = 0.d0 + r(3) = zcenter -zmax * 0.5d0 print*,'r(3) = ',r(3) istate = 1 + + write(i_unit_output,*)" z, on-top(z), n(z) " do i = 1, nz call give_on_top_in_r_one_state(r,istate,on_top_in_r) - call give_cas_density_in_r(r,total_density) + call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_density,r) write(i_unit_output,*)r(3),on_top_in_r,total_density r(3) += dz enddo - end diff --git a/src/dft_utils_in_r/dm_in_r_routines.irp.f b/src/dft_utils_in_r/dm_in_r_routines.irp.f index 6fa99e22..9991289c 100644 --- a/src/dft_utils_in_r/dm_in_r_routines.irp.f +++ b/src/dft_utils_in_r/dm_in_r_routines.irp.f @@ -110,6 +110,62 @@ end grad_dm_b *= 2.d0 end + subroutine density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) + implicit none + BEGIN_DOC +! input: +! +! * r(1) ==> r(1) = x, r(2) = y, r(3) = z +! +! output: +! +! * dm_a = alpha density evaluated at r +! * dm_b = beta density evaluated at r +! * grad_dm_a(1) = X gradient of the alpha density evaluated in r +! * grad_dm_a(1) = X gradient of the beta density evaluated in r +! + END_DOC + double precision, intent(in) :: r(3) + double precision, intent(out) :: dm_a(N_states),dm_b(N_states) + double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states) + double precision :: grad_aos_array(3,ao_num) + integer :: i,j,istate + double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v + double precision :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3) + + call give_all_aos_and_grad_at_r(r,aos_array,grad_aos_array) + do i = 1, ao_num + do j = 1, 3 + aos_grad_array(i,j) = grad_aos_array(j,i) + enddo + enddo + + do istate = 1, N_states + ! alpha density + ! aos_array_bis = \rho_ao * aos_array + call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1) + dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num) + + ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i) + grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num) + grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num) + grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num) + ! aos_grad_array_bis = \rho_ao * aos_grad_array + + ! beta density + call dsymv('U',ao_num,1.d0,one_e_dm_beta_ao_for_dft(1,1,istate),size(one_e_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1) + dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num) + + ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i) + grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num) + grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num) + grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num) + ! aos_grad_array_bis = \rho_ao * aos_grad_array + enddo + grad_dm_a *= 2.d0 + grad_dm_b *= 2.d0 + end + subroutine density_and_grad_lapl_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, lapl_dm_a, lapl_dm_b, aos_array, grad_aos_array, lapl_aos_array) diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg index 4a56a35b..694590ec 100644 --- a/src/scf_utils/EZFIO.cfg +++ b/src/scf_utils/EZFIO.cfg @@ -51,3 +51,10 @@ doc: If true, leave untouched all the orbitals defined as core and optimize all interface: ezfio,provider,ocaml default: False + +[no_oa_or_av_opt] +type: logical +doc: If true, you set to zero all Fock elements between the orbital set to active and all the other orbitals +interface: ezfio,provider,ocaml +default: False + diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f index d501278f..5188581a 100644 --- a/src/scf_utils/diagonalize_fock.irp.f +++ b/src/scf_utils/diagonalize_fock.irp.f @@ -31,6 +31,27 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num) enddo enddo endif + if(no_oa_or_av_opt)then + 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 diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index fc9eaadd..61633d3b 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -92,6 +92,27 @@ enddo endif + if(no_oa_or_av_opt)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + Fock_matrix_mo(iorb,jorb) = 0.d0 + Fock_matrix_mo(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + Fock_matrix_mo(iorb,jorb) = 0.d0 + Fock_matrix_mo(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + Fock_matrix_mo(iorb,jorb) = 0.d0 + Fock_matrix_mo(jorb,iorb) = 0.d0 + enddo + enddo + endif + END_PROVIDER From 6ecb8d1be6c5a7a7183e5306d7363e79466c5d8c Mon Sep 17 00:00:00 2001 From: Gilles Frison Date: Thu, 15 Apr 2021 15:46:41 +0200 Subject: [PATCH 21/26] added opam_installer_no_usr_bin.sh --- scripts/opam_installer_no_usr_bin.sh | 365 +++++++++++++++++++++++++++ 1 file changed, 365 insertions(+) create mode 100755 scripts/opam_installer_no_usr_bin.sh diff --git a/scripts/opam_installer_no_usr_bin.sh b/scripts/opam_installer_no_usr_bin.sh new file mode 100755 index 00000000..0f7219d1 --- /dev/null +++ b/scripts/opam_installer_no_usr_bin.sh @@ -0,0 +1,365 @@ +#!/bin/sh + +set -ue + +# (c) Copyright Fabrice Le Fessant INRIA/OCamlPro 2013 +# (c) Copyright Louis Gesbert OCamlPro 2014-2017 + +VERSION='2.0.7' +DEV_VERSION='2.1.0~beta2' +DEFAULT_BINDIR=/usr/local/bin + +bin_sha512() { + case "$OPAM_BIN" in + opam-2.0.6-arm64-linux) echo "d2b3d92fd5fae7f053702b53ddbc7c224fcfbfc9b232247ba4e40cbf1cda28f160d8c14fde87aebeebfd2545e13265c0ee4a47e292f035767fb944b1b8ff5c90";; + opam-2.0.6-armhf-linux) echo "a42a7ad8c1afdb20ac5746934306576e6364f5453b176ccd42a3e5a116a5db05c2758cec31800ffab11411290cf671f9eee3f299df48c7ceca8e4d7e33dfedc8";; + opam-2.0.6-i686-linux) echo "6c0d965f89a2026ead3120e217d12b2df7426740d54bc94e2c46faaeff5893081e68aac162621bfa694ab597a18be28165f10cdda1217a4d73653789a9928b64";; + opam-2.0.6-x86_64-linux) echo "2b9d4a99aa28a193c88c7c6f6265203bd3cfeef98929d6f5cfce4b52cd9ddbd7be7eddc1d3d9c440f81d65074dd7851b8d29cd397fb06d2cfccffb54d3cdcc6a";; + opam-2.0.6-x86_64-macos) echo "cf02546b22ca91b1d97a3657b970b34d4acf4dc745696b7200ff185d25ebb5914ea8b6a94b503eb8c999634de6fdb944998a970105cd6a4c6df538c262b48b7f";; + opam-2.0.6-x86_64-openbsd) echo "2f58b3d4902d4c3fb823d251a50e034f9101b0c5a3827725876bb3bcb6c013c4f54138054d82abba0a9e917675275e26f05b98630cf7116c465d2110756f1309";; + + opam-2.0.7-arm64-linux) echo "0dd4d80496545f684af39dc5b4b28867bc19a74186577c38bd2a8934d871c2cbcdb9891bfd41c080b5f12d5a3c8801e203df8a76d55e1e22fe80d31447402e46";; + opam-2.0.7-armhf-linux) echo "ea691bc9565acc1207dea3dfb89192b1865b5b5809efe804a329f39878640fb19771edcb05c5699f8e914e88e3155f31132b845c54b0095bedd3952d336bae0b";; + opam-2.0.7-i686-linux) echo "5fa8fb9664d36ead5760e7e1c337f6ae7b0fd4be5089ddfb50ae74028deec30893b1f4dee040402bc3f15da197ba89a45c7d626ecf6e5be80d176f43526c4bad";; + opam-2.0.7-x86_64-linux) echo "da75b0cb5ad50f95d31857a7d72f1836132a1fa1cdbfdedf684342b798e7107b4add4c74c05d5ce44881309fa1e57707538dbcda874e7f74b269b1bb204f3ae3";; + opam-2.0.7-x86_64-macos) echo "de1194c8e97e53956e5e47502c28881bbf26d1beaac4f33a43a922b8ca7ce97725533cfaf65a33fc0e183eab5a95e9ecd2e20f72faeaec333dc3850b79b5fe8a";; + opam-2.0.7-x86_64-openbsd) echo "b253809c4388847e1a33b5c4f1f5d72bef79a2f0c43b19ef65b40d0c10341aa0bee4a4b1f3a9ab70eb026e4cc220a63cfc56a18c035b6b0297c92f2bdb7f9a78";; + + opam-2.1.0-alpha-arm64-linux) echo "1bf0acfa64aa01c3244e65eed60eef1caaa6de53aa8b32dd0d2446f91905a1e41591f53cd350e85b2b9f5edba9b137d723c32949115623e9753e77b707bb25b0";; + opam-2.1.0-alpha-armhf-linux) echo "87c12a422bd14a0d10a94ddaaa46de23700e3b89810a0c06232eff8d96b37c2fd43dcb5a8da5a2004aa8040d1b93293209f1ff1aab865ffd150364e24c87c716";; + opam-2.1.0-alpha-i686-linux) echo "b8369da6d4795a461ff1b49e687b027325d4e90bc8f19517e52a94ee3be167c4faaaf33bd0b3536be552d2add54865d0e33933acaa674f2e1a17249b022738af";; + opam-2.1.0-alpha-x86_64-linux) echo "2e22747829fb0bada3a74a23f5e0ff2228520d647fc4fe08a1ce76f3cb357cc7240f7b45e422c5f4b8eafe832ae3a8973ecbd4814ae0e8ce1096bcff39482020";; + opam-2.1.0-alpha-x86_64-macos) echo "c440e8ae1970fa7533e6e1b96ba3e3dd65b04432d41bc57ce4c768ed9b4229954546d59ec06f3d4ee49cbe00bb4bfd0b3f509d6d9a27de2db17725e097a61c86";; + opam-2.1.0-alpha-x86_64-openbsd) echo "d87afe99fee541a1c6fae30b72653db7a5ea2abdec3fa3b2b480daddf3fcd8d4096e2a40458310755faec3722119f29ed981ffbfa65142e618f99b70572f892f";; + + opam-2.1.0-alpha2-arm64-linux) echo "b67520bb2a6c59f800da100278d74e58f2bbf66924f94643023dc46b97b16f17a30de95e439c6f9b032bd555c062ddba325f3e5169cac186615b959a8c434788";; + opam-2.1.0-alpha2-armhf-linux) echo "9a6312eb54d6c9c2036ca90f7816789c27c23f1b1d325cd69d27a910cdd8760b82f19c9e9b61b5b6214818f1f40f8b4d2ef081acb43f0dad68c976986a7c6a45";; + opam-2.1.0-alpha2-i686-linux) echo "0dc07f236405777ad74d58fcc6cb6c3247e7dfc31408df4a199599077d5cb41ec86895f1d0c5eaa2a9c70842a2a998226674f986ba0044c82896c073ac90b209";; + opam-2.1.0-alpha2-x86_64-linux) echo "21509e8abd8463f4e18a55398f690700772e25f0ddb9f3fd7644e2f9a9a89ebbf5c09efbeceafe4a0ab5015d0d03b2f29506be514aae813a2f3dac7dd01261f3";; + opam-2.1.0-alpha2-x86_64-macos) echo "1c1bd26621eebb5bf3783dec80d5555aa5ff02dcbf43eb44398798e6162c1964bc1964e3980391ea115e5c068c1bb66960f8ebdd91bc4f0bac844f3a61433f1e";; + opam-2.1.0-alpha2-x86_64-openbsd) echo "941f3e306bc36e8e44e4245ca5e635b04e0a54f33439d55d41875ced47384cad8c222b649027d3c4eacc3c2c569cf5006c872763b19c490d9b289c9cfe4f491a";; + + opam-2.1.0-alpha3-arm64-linux) echo "ad906bb2ab764a92fabdf0b906310c5034bf5daf0ebfb2529e9b87661ddbf8fd14f51dee5ce75b4fd4bb5789e29c7be71063f1ebcc92e92333be12aa62efdff9";; + opam-2.1.0-alpha3-armhf-linux) echo "2a7022c1f5dbc855a0d067f29677b13253dccbc9792b8170fa72a743802bbcd6e41ce7512c4845091af0f73b8ba7573038ec53ea9aaf74be04367ac1767e7220";; + opam-2.1.0-alpha3-i686-linux) echo "6f2fce0c45ae700e7a1b32d0a24988645c9aed3afc45998c8fbe70e97a65e3ba5d824069914a892bb3f9b1336383cfd492c28678ff16db5cada863da924b07d8";; + opam-2.1.0-alpha3-x86_64-linux) echo "1d219dbf670e1550bf71c28e586d14f1d8af2605f0e13bea2f11ad52a7f176bd9a89637e44a91a024f0088db1b2aba8dc3207bc81fa930580e54f4031255c178";; + opam-2.1.0-alpha3-x86_64-macos) echo "93edb6c1151f8f5bd017f230ffd9277f6ad943e3f5032ea000c37f012738fb3ab4b4add172e1f624c37e6564963fef0716b876b0113c8e43f5943d77bbbc173c";; + opam-2.1.0-alpha3-x86_64-openbsd) echo "0e3b3761e877c57f5b333aacb70c86bf60f50eecdca6e9e1a552e3d666cea034d8873f3a87e585a5970b1aef7e540adb18c71e0e8fd8794843dd5d1d421a87ec";; + + opam-2.1.0-beta-arm64-linux) echo "954670c74ea8244b440756e4f7755bd2b5548ab67428ce577c4c507fc33c8d00eb73c4d7b59ccb0ef800f4465b5c704573c63486b78a23e9568f3751bf9aef78";; + opam-2.1.0-beta-armhf-linux) echo "cc666f2c6b1ac07d1bc8a035c6b3a9455794b51a827c54bb92786ae1a75c6c55839d3f48b378508f42a66ac887fdc68f7628a67e2826813cb6df048c906755ca";; + opam-2.1.0-beta-i686-linux) echo "66ac48b298741f753ca868be362851ccd9bf84fd8772d18f3307e99cf72c8c68ac9fa17bf2d610d7f3b5dc6209eb8371bf0e10b363e963fc6c31d70e5938017f";; + opam-2.1.0-beta-x86_64-linux) echo "e316f1b5f1c668affba6c2819f692c28776e131a17fb64b2c0e23f8a3b7d456575a8109fcdcb9babfad13bc33c17fa619cbb4a48ca6198765f86296b7e611f24";; + opam-2.1.0-beta-x86_64-macos) echo "acb29b7c64df314c6629e14f6d8f079504d39b7fd3104867fd22df3395ccfea9f1014a3a87dff9c12bf03ca451e9ee2918b9d9d8f17ce1a6d7de0c0649452fa9";; + opam-2.1.0-beta-x86_64-openbsd) echo "ff9fa1ee0ae7e54b4e18999cf5ea9b899c0b4039b744a950e96221e3e86c21eaa50904bdbc836ff8103f7713506d0de3d32ec77b169561e0cd694bfeea812cae";; + + opam-2.1.0-beta2-arm64-linux) echo "a58ba3ebb4431d3cabfe96b806c9897205153e8a546ebe74f0229982758d140b4fcbcea421db70589b1eb3080dc86534522a3cba0330ce82e0898a60048d51ba";; + opam-2.1.0-beta2-armhf-linux) echo "fc4e6b753ce6368f75a0d3005f4b21ce9606599d21607a67015db55a38b6ef473b4205f5b128c5808189feed8ae58f93bd79348988be7c5007ae1b39307a5cd0";; + opam-2.1.0-beta2-i686-linux) echo "a376a6e0e1e2b08ea4d0a5c1c38487e67984bef2e89f978536dd08283f945f74dd31ee287bc68d91690603ba0fa657e91ff0d30bea217743f79ed99d2390eba5";; + opam-2.1.0-beta2-x86_64-linux) echo "12c5e2b0087ed389fa12fdb0e1f6f7dc0b3df3f95c59e8bc576279b7780921d47bbc4ebcba6caddde30f4fb1cc9e4a873cc8a6aef80fcc48a878aba69be7af44";; + opam-2.1.0-beta2-x86_64-macos) echo "4acc12672a2e3ad7e78540634edcae2e7e84860057b86a56b1cdf7eacf8d97957aaa864f571d6fb8f61ee8280f8a4ed73df7881d91a22c9d8c2d73e8a558f61d";; + opam-2.1.0-beta2-x86_64-openbsd) echo "84d7d409220c72e3ed7e6acdd7cce3b5a208f2966d232648a57a48641ab8ce4fa58e94e40b7176201455d82260e6c501a6ba4a30b1426a552f8d09cfd027ddde";; + + *) echo "no sha";; + esac +} + +usage() { + echo "opam binary installer v.$VERSION" + echo "Downloads and installs a pre-compiled binary of opam $VERSION to the system." + echo "This can also be used to switch between opam versions" + echo + echo "Options:" + echo " --dev Install the latest alpha or beta instead: $DEV_VERSION" + echo " --no-backup Don't attempt to backup the current opam root" + echo " --backup Force the backup the current opam root (even if it" + echo " is from the 2.0 branch already)" + echo " --fresh Create the opam $VERSION root from scratch" + echo " --restore VERSION Restore a backed up opam binary and root" + echo + echo "The default is to backup if the current version of opam is 1.*, or when" + echo "using '--fresh' or '--dev'" +} + +RESTORE= +NOBACKUP= +FRESH= +DOWNLOAD_ONLY= + +while [ $# -gt 0 ]; do + case "$1" in + --dev) + VERSION=$DEV_VERSION + if [ -z "$NOBACKUP" ]; then NOBACKUP=0; fi;; + --restore) + if [ $# -lt 2 ]; then echo "Option $1 requires an argument"; exit 2; fi + shift; + RESTORE=$1;; + --no-backup) + NOBACKUP=1;; + --backup) + NOBACKUP=0;; + --fresh) + FRESH=1;; + --download-only) + DOWNLOAD_ONLY=1;; + --help|-h) + usage; exit 0;; + *) + usage; exit 2;; + esac + shift +done + + +TMP=${TMPDIR:-/tmp} + +ARCH=$(uname -m || echo unknown) +case "$ARCH" in + x86|i?86) ARCH="i686";; + x86_64|amd64) ARCH="x86_64";; + ppc|powerpc|ppcle) ARCH="ppc";; + aarch64_be|aarch64|armv8b|armv8l) ARCH="arm64";; + armv5*|armv6*|earmv6*|armv7*|earmv7*) ARCH="armhf";; + *) ARCH=$(echo "$ARCH" | awk '{print tolower($0)}') +esac + +OS=$( (uname -s || echo unknown) | awk '{print tolower($0)}') + +if [ "$OS" = "darwin" ] ; then + OS=macos +fi + +TAG=$(echo "$VERSION" | tr '~' '-') + +OPAM_BIN_URL_BASE='https://github.com/ocaml/opam/releases/download/' +OPAM_BIN="opam-${TAG}-${ARCH}-${OS}" +OPAM_BIN_URL="${OPAM_BIN_URL_BASE}${TAG}/${OPAM_BIN}" + +download() { + if command -v wget >/dev/null; then wget -q -O "$@" + else curl -s -L -o "$@" + fi +} + +check_sha512() { + OPAM_BIN_LOC="$1" + if command -v openssl > /dev/null; then + sha512_devnull="cf83e1357eefb8bdf1542850d66d8007d620e4050b5715dc83f4a921d36ce9ce47d0d13c5d85f2b0ff8318d2877eec2f63b931bd47417a81a538327af927da3e" + sha512_check=`openssl sha512 2>&1 < /dev/null | cut -f 2 -d ' '` + if [ "x$sha512_devnull" = "x$sha512_check" ]; then + sha512=`openssl sha512 "$OPAM_BIN_LOC" 2> /dev/null | cut -f 2 -d ' '` + check=`bin_sha512` + test "x$sha512" = "x$check" + else + echo "openssl 512 option not handled, binary integrity check can't be performed." + return 0 + fi + else + echo "openssl not found, binary integrity check can't be performed." + return 0 + fi +} + +download_and_check() { + OPAM_BIN_LOC="$1" + echo "## Downloading opam $VERSION for $OS on $ARCH..." + + if ! download "$OPAM_BIN_LOC" "$OPAM_BIN_URL"; then + echo "There may not yet be a binary release for your architecture or OS, sorry." + echo "See https://github.com/ocaml/opam/releases/tag/$TAG for pre-compiled binaries," + echo "or run 'make cold' from https://github.com/ocaml/opam/archive/$TAG.tar.gz" + echo "to build from scratch" + exit 10 + else + if check_sha512 "$OPAM_BIN_LOC"; then + echo "## Downloaded." + else + echo "Checksum mismatch, a problem occurred during download." + exit 10 + fi + fi +} + +DOWNLOAD_ONLY=${DOWNLOAD_ONLY:-0} + +if [ $DOWNLOAD_ONLY -eq 1 ]; then + OPAM_BIN_LOC="$PWD/$OPAM_BIN" + if [ -e "$OPAM_BIN_LOC" ]; then + echo "Found opam binary in $OPAM_BIN_LOC ..." + if check_sha512 "$OPAM_BIN_LOC" ; then + echo "... with matching sha512" + exit 0; + else + echo "... with mismatching sha512, download the good one." + fi + fi + download_and_check "$OPAM_BIN_LOC" + exit 0; +fi + +EXISTING_OPAM=$(command -v opam || echo) +EXISTING_OPAMV= +if [ -n "$EXISTING_OPAM" ]; then + EXISTING_OPAMV=$("$EXISTING_OPAM" --version || echo "unknown") +fi + +FRESH=${FRESH:-0} + +OPAMROOT=${OPAMROOT:-$HOME/.opam} + +if [ ! -d "$OPAMROOT" ]; then FRESH=1; fi + +if [ -z "$NOBACKUP" ] && [ ! "$FRESH" = 1 ] && [ -z "$RESTORE" ]; then + case "$EXISTING_OPAMV" in + 2.*) NOBACKUP=1;; + *) NOBACKUP=0;; + esac +fi + +xsudo() { + local CMD=$1; shift + local DST + for DST in "$@"; do : ; done + + local DSTDIR=$(dirname "$DST") + if [ ! -w "$DSTDIR" ]; then + echo "Write access to $DSTDIR required, using 'sudo'." + echo "Command: $CMD $@" + if [ "$CMD" = "install" ]; then + sudo "$CMD" -g 0 -o root "$@" + else + sudo "$CMD" "$@" + fi + else + "$CMD" "$@" + fi +} + +if [ -n "$RESTORE" ]; then + OPAM=$(command -v opam) + OPAMV=$("$OPAM" --version) + OPAM_BAK="$OPAM.$RESTORE" + OPAMROOT_BAK="$OPAMROOT.$RESTORE" + if [ ! -e "$OPAM_BAK" ] || [ ! -d "$OPAMROOT_BAK" ]; then + echo "No backup of opam $RESTORE was found" + exit 1 + fi + if [ "$NOBACKUP" = 1 ]; then + printf "## This will clear $OPAM and $OPAMROOT. Continue ? [Y/n] " + read R + case "$R" in + ""|"y"|"Y"|"yes") + xsudo rm -f "$OPAM" + rm -rf "$OPAMROOT";; + *) exit 1 + esac + else + xsudo mv "$OPAM" "$OPAM.$OPAMV" + mv "$OPAMROOT" "$OPAMROOT.$OPAMV" + fi + xsudo mv "$OPAM_BAK" "$OPAM" + mv "$OPAMROOT_BAK" "$OPAMROOT" + printf "## Opam $RESTORE and its root were restored." + if [ "$NOBACKUP" = 1 ]; then echo + else echo " Opam $OPAMV was backed up." + fi + exit 0 +fi + +#if [ -e "$TMP/$OPAM_BIN" ] && ! check_sha512 "$TMP/$OPAM_BIN" || [ ! -e "$TMP/$OPAM_BIN" ]; then + download_and_check "$TMP/$OPAM_BIN" +#else +# echo "## Using already downloaded \"$TMP/$OPAM_BIN\"" +#fi + +if [ -n "$EXISTING_OPAM" ]; then + DEFAULT_BINDIR=$(dirname "$EXISTING_OPAM") +fi + +while true; do + printf "## Where should it be installed ? [$DEFAULT_BINDIR] " + read BINDIR + if [ -z "$BINDIR" ]; then BINDIR="$DEFAULT_BINDIR"; fi + + if [ -d "$BINDIR" ]; then break + else + printf "## $BINDIR does not exist. Create ? [Y/n] " + read R + case "$R" in + ""|"y"|"Y"|"yes") + xsudo mkdir -p $BINDIR + break;; + esac + fi +done + +#if [ -e "$EXISTING_OPAM" ]; then +# if [ "$NOBACKUP" = 1 ]; then +# xsudo rm -f "$EXISTING_OPAM" +# else +# xsudo mv "$EXISTING_OPAM" "$EXISTING_OPAM.$EXISTING_OPAMV" +# echo "## $EXISTING_OPAM backed up as $(basename $EXISTING_OPAM).$EXISTING_OPAMV" +# fi +#fi + +if [ -d "$OPAMROOT" ]; then + if [ "$FRESH" = 1 ]; then + if [ "$NOBACKUP" = 1 ]; then + printf "## This will clear $OPAMROOT. Continue ? [Y/n] " + read R + case "$R" in + ""|"y"|"Y"|"yes") + rm -rf "$OPAMROOT";; + *) exit 1 + esac + else + mv "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV" + echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV" + fi + echo "## opam $VERSION installed. Please run 'opam init' to get started" + elif [ ! "$NOBACKUP" = 1 ]; then + echo "## Backing up $OPAMROOT to $(basename $OPAMROOT).$EXISTING_OPAMV (this may take a while)" + if [ -e "$OPAMROOT.$EXISTING_OPAMV" ]; then + echo "ERROR: there is already a backup at $OPAMROOT.$EXISTING_OPAMV" + echo "Please move it away or run with --no-backup" + fi + FREE=$(df -k "$OPAMROOT" | awk 'NR>1 {print $4}') + NEEDED=$(du -sk "$OPAMROOT" | awk '{print $1}') + if ! [ $NEEDED -lt $FREE ]; then + echo "Error: not enough free space to backup. You can retry with --no-backup," + echo "--fresh, or remove '$OPAMROOT'" + exit 1 + fi + cp -a "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV" + echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV" + fi + rm -f "$OPAMROOT"/repo/*/*.tar.gz* +fi + +xsudo install -m 755 "$TMP/$OPAM_BIN" "$BINDIR/opam" +echo "## opam $VERSION installed to $BINDIR" + +if [ ! "$FRESH" = 1 ]; then + echo "## Converting the opam root format & updating" + "$BINDIR/opam" init --reinit -ni +fi + +WHICH=$(command -v opam || echo notfound) + +case "$WHICH" in + "$BINDIR/opam") ;; + notfound) echo "## Remember to add $BINDIR to your PATH";; + *) + echo "## WARNING: 'opam' command found in PATH does not match the installed one:" + echo " - Installed: '$BINDIR/opam'" + echo " - Found: '$WHICH'" + echo "Make sure to remove the second or fix your PATH to use the new opam" + echo +esac + +if [ ! "$NOBACKUP" = 1 ]; then + echo "## Run this script again with '--restore $EXISTING_OPAMV' to revert." +fi + +rm -f $TMP/$OPAM_BIN From 06e97356d3bdec7ce9c3788dacc1257e6380e1fe Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 19 May 2021 15:23:06 +0200 Subject: [PATCH 22/26] added additional grid --- src/becke_numerical_grid/EZFIO.cfg | 31 +++ .../angular_extra_grid.irp.f | 104 ++++++++++ src/becke_numerical_grid/extra_grid.irp.f | 178 ++++++++++++++++++ 3 files changed, 313 insertions(+) create mode 100644 src/becke_numerical_grid/angular_extra_grid.irp.f create mode 100644 src/becke_numerical_grid/extra_grid.irp.f diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index 17e8fd90..53b11fbd 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -33,3 +33,34 @@ doc: Number of angular grid points given from input. Warning, this number cannot interface: ezfio,provider,ocaml default: 1202 + +[extra_grid_type_sgn] +type: integer +doc: Type of extra_grid used for the Becke's numerical extra_grid. Can be, by increasing accuracy: [ 0 | 1 | 2 | 3 ] +interface: ezfio,provider,ocaml +default: 2 + +[thresh_extra_grid] +type: double precision +doc: threshold on the weight of a given extra_grid point +interface: ezfio,provider,ocaml +default: 1.e-20 + +[my_extra_grid_becke] +type: logical +doc: if True, the number of angular and radial extra_grid points are read from EZFIO +interface: ezfio,provider,ocaml +default: False + +[my_n_pt_r_extra_grid] +type: integer +doc: Number of radial extra_grid points given from input +interface: ezfio,provider,ocaml +default: 300 + +[my_n_pt_a_extra_grid] +type: integer +doc: Number of angular extra_grid points given from input. Warning, this number cannot be any integer. See file list_angular_extra_grid +interface: ezfio,provider,ocaml +default: 1202 + diff --git a/src/becke_numerical_grid/angular_extra_grid.irp.f b/src/becke_numerical_grid/angular_extra_grid.irp.f new file mode 100644 index 00000000..bbf869e2 --- /dev/null +++ b/src/becke_numerical_grid/angular_extra_grid.irp.f @@ -0,0 +1,104 @@ + + BEGIN_PROVIDER [double precision, angular_quadrature_points_extra, (n_points_extra_integration_angular,3) ] +&BEGIN_PROVIDER [double precision, weights_angular_points_extra, (n_points_extra_integration_angular)] + implicit none + BEGIN_DOC + ! weights and grid points_extra for the integration on the angular variables on + ! the unit sphere centered on (0,0,0) + ! According to the LEBEDEV scheme + END_DOC + + include 'constants.include.F' + integer :: i + double precision :: accu + double precision :: degre_rad + double precision :: x(n_points_extra_integration_angular) + double precision :: y(n_points_extra_integration_angular) + double precision :: z(n_points_extra_integration_angular) + double precision :: w(n_points_extra_integration_angular) + + degre_rad = pi/180.d0 + accu = 0.d0 + + select case (n_points_extra_integration_angular) + + + case (0006) + call LD0006(X,Y,Z,W,n_points_extra_integration_angular) + case (0014) + call LD0014(X,Y,Z,W,n_points_extra_integration_angular) + case (0026) + call LD0026(X,Y,Z,W,n_points_extra_integration_angular) + case (0038) + call LD0038(X,Y,Z,W,n_points_extra_integration_angular) + case (0050) + call LD0050(X,Y,Z,W,n_points_extra_integration_angular) + case (0074) + call LD0074(X,Y,Z,W,n_points_extra_integration_angular) + case (0086) + call LD0086(X,Y,Z,W,n_points_extra_integration_angular) + case (0110) + call LD0110(X,Y,Z,W,n_points_extra_integration_angular) + case (0146) + call LD0146(X,Y,Z,W,n_points_extra_integration_angular) + case (0170) + call LD0170(X,Y,Z,W,n_points_extra_integration_angular) + case (0194) + call LD0194(X,Y,Z,W,n_points_extra_integration_angular) + case (0230) + call LD0230(X,Y,Z,W,n_points_extra_integration_angular) + case (0266) + call LD0266(X,Y,Z,W,n_points_extra_integration_angular) + case (0302) + call LD0302(X,Y,Z,W,n_points_extra_integration_angular) + case (0350) + call LD0350(X,Y,Z,W,n_points_extra_integration_angular) + case (0434) + call LD0434(X,Y,Z,W,n_points_extra_integration_angular) + case (0590) + call LD0590(X,Y,Z,W,n_points_extra_integration_angular) + case (0770) + call LD0770(X,Y,Z,W,n_points_extra_integration_angular) + case (0974) + call LD0974(X,Y,Z,W,n_points_extra_integration_angular) + case (1202) + call LD1202(X,Y,Z,W,n_points_extra_integration_angular) + case (1454) + call LD1454(X,Y,Z,W,n_points_extra_integration_angular) + case (1730) + call LD1730(X,Y,Z,W,n_points_extra_integration_angular) + case (2030) + call LD2030(X,Y,Z,W,n_points_extra_integration_angular) + case (2354) + call LD2354(X,Y,Z,W,n_points_extra_integration_angular) + case (2702) + call LD2702(X,Y,Z,W,n_points_extra_integration_angular) + case (3074) + call LD3074(X,Y,Z,W,n_points_extra_integration_angular) + case (3470) + call LD3470(X,Y,Z,W,n_points_extra_integration_angular) + case (3890) + call LD3890(X,Y,Z,W,n_points_extra_integration_angular) + case (4334) + call LD4334(X,Y,Z,W,n_points_extra_integration_angular) + case (4802) + call LD4802(X,Y,Z,W,n_points_extra_integration_angular) + case (5294) + call LD5294(X,Y,Z,W,n_points_extra_integration_angular) + case (5810) + call LD5810(X,Y,Z,W,n_points_extra_integration_angular) + case default + print *, irp_here//': wrong n_points_extra_integration_angular. See in ${QP_ROOT}/src/becke_numerical_grid/list_angular_grid to see the possible angular grid points_extra. Ex: ' + print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]' + stop -1 + end select + + do i = 1, n_points_extra_integration_angular + angular_quadrature_points_extra(i,1) = x(i) + angular_quadrature_points_extra(i,2) = y(i) + angular_quadrature_points_extra(i,3) = z(i) + weights_angular_points_extra(i) = w(i) * 4.d0 * pi + accu += w(i) + enddo + +END_PROVIDER diff --git a/src/becke_numerical_grid/extra_grid.irp.f b/src/becke_numerical_grid/extra_grid.irp.f new file mode 100644 index 00000000..db691e55 --- /dev/null +++ b/src/becke_numerical_grid/extra_grid.irp.f @@ -0,0 +1,178 @@ + + BEGIN_PROVIDER [integer, n_points_extra_radial_grid] +&BEGIN_PROVIDER [integer, n_points_extra_integration_angular] + implicit none + BEGIN_DOC + ! n_points_extra_radial_grid = number of radial grid points_extra per atom + ! + ! n_points_extra_integration_angular = number of angular grid points_extra per atom + ! + ! These numbers are automatically set by setting the grid_type_sgn parameter + END_DOC +if(.not.my_extra_grid_becke)then + select case (extra_grid_type_sgn) + case(0) + n_points_extra_radial_grid = 23 + n_points_extra_integration_angular = 170 + case(1) + n_points_extra_radial_grid = 50 + n_points_extra_integration_angular = 194 + case(2) + n_points_extra_radial_grid = 75 + n_points_extra_integration_angular = 302 + case(3) + n_points_extra_radial_grid = 99 + n_points_extra_integration_angular = 590 + case default + write(*,*) '!!! Quadrature grid not available !!!' + stop + end select +else + n_points_extra_radial_grid = my_n_pt_r_extra_grid + n_points_extra_integration_angular = my_n_pt_a_extra_grid +endif +END_PROVIDER + +BEGIN_PROVIDER [integer, n_points_extra_grid_per_atom] + implicit none + BEGIN_DOC + ! Number of grid points_extra per atom + END_DOC + n_points_extra_grid_per_atom = n_points_extra_integration_angular * n_points_extra_radial_grid + +END_PROVIDER + + BEGIN_PROVIDER [double precision, grid_points_extra_radial, (n_points_extra_radial_grid)] +&BEGIN_PROVIDER [double precision, dr_radial_extra_integral] + + implicit none + BEGIN_DOC + ! points_extra in [0,1] to map the radial integral [0,\infty] + END_DOC + dr_radial_extra_integral = 1.d0/dble(n_points_extra_radial_grid-1) + integer :: i + do i = 1, n_points_extra_radial_grid + grid_points_extra_radial(i) = dble(i-1) * dr_radial_extra_integral + enddo + +END_PROVIDER + +BEGIN_PROVIDER [double precision, grid_points_extra_per_atom, (3,n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)] + BEGIN_DOC + ! x,y,z coordinates of grid points_extra used for integration in 3d space + END_DOC + implicit none + integer :: i,j,k + double precision :: dr,x_ref,y_ref,z_ref + double precision :: knowles_function + do i = 1, nucl_num + x_ref = nucl_coord(i,1) + y_ref = nucl_coord(i,2) + z_ref = nucl_coord(i,3) + do j = 1, n_points_extra_radial_grid-1 + double precision :: x,r + ! x value for the mapping of the [0, +\infty] to [0,1] + x = grid_points_extra_radial(j) + + ! value of the radial coordinate for the integration + r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x) + + ! explicit values of the grid points_extra centered around each atom + do k = 1, n_points_extra_integration_angular + grid_points_extra_per_atom(1,k,j,i) = & + x_ref + angular_quadrature_points_extra(k,1) * r + grid_points_extra_per_atom(2,k,j,i) = & + y_ref + angular_quadrature_points_extra(k,2) * r + grid_points_extra_per_atom(3,k,j,i) = & + z_ref + angular_quadrature_points_extra(k,3) * r + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] + BEGIN_DOC + ! Weight function at grid points_extra : w_n(r) according to the equation (22) + ! of Becke original paper (JCP, 88, 1988) + ! + ! The "n" discrete variable represents the nucleis which in this array is + ! represented by the last dimension and the points_extra are labelled by the + ! other dimensions. + END_DOC + implicit none + integer :: i,j,k,l,m + double precision :: r(3) + double precision :: accu,cell_function_becke + double precision :: tmp_array(nucl_num) + ! run over all points_extra in space + ! that are referred to each atom + do j = 1, nucl_num + !for each radial grid attached to the "jth" atom + do k = 1, n_points_extra_radial_grid -1 + ! for each angular point attached to the "jth" atom + do l = 1, n_points_extra_integration_angular + r(1) = grid_points_extra_per_atom(1,l,k,j) + r(2) = grid_points_extra_per_atom(2,l,k,j) + r(3) = grid_points_extra_per_atom(3,l,k,j) + accu = 0.d0 + ! For each of these points_extra in space, ou need to evaluate the P_n(r) + do i = 1, nucl_num + ! function defined for each atom "i" by equation (13) and (21) with k == 3 + tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + ! Then you compute the summ the P_n(r) function for each of the "r" points_extra + accu += tmp_array(i) + enddo + accu = 1.d0/accu + weight_at_r_extra(l,k,j) = tmp_array(j) * accu + if(isnan(weight_at_r_extra(l,k,j)))then + print*,'isnan(weight_at_r_extra(l,k,j))' + print*,l,k,j + accu = 0.d0 + do i = 1, nucl_num + ! function defined for each atom "i" by equation (13) and (21) with k == 3 + tmp_array(i) = cell_function_becke(r,i) ! P_n(r) + print*,i,tmp_array(i) + ! Then you compute the summ the P_n(r) function for each of the "r" points_extra + accu += tmp_array(i) + enddo + write(*,'(100(F16.10,X))')tmp_array(j) , accu + stop + endif + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, final_weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] + BEGIN_DOC + ! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights. + END_DOC + implicit none + integer :: i,j,k,l,m + double precision :: r(3) + double precision :: accu,cell_function_becke + double precision :: tmp_array(nucl_num) + double precision :: contrib_integration,x + double precision :: derivative_knowles_function,knowles_function + ! run over all points_extra in space + do j = 1, nucl_num ! that are referred to each atom + do i = 1, n_points_extra_radial_grid -1 !for each radial grid attached to the "jth" atom + x = grid_points_extra_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] + do k = 1, n_points_extra_integration_angular ! for each angular point attached to the "jth" atom + contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)& + *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2 + final_weight_at_r_extra(k,i,j) = weights_angular_points_extra(k) * weight_at_r_extra(k,i,j) * contrib_integration * dr_radial_extra_integral + if(isnan(final_weight_at_r_extra(k,i,j)))then + print*,'isnan(final_weight_at_r_extra(k,i,j))' + print*,k,i,j + write(*,'(100(F16.10,X))')weights_angular_points_extra(k) , weight_at_r_extra(k,i,j) , contrib_integration , dr_radial_extra_integral + stop + endif + enddo + enddo + enddo + +END_PROVIDER + From 95d470ea52ff4e3a43ad036e9c23aec9ceda23ef Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 19 May 2021 16:12:41 +0200 Subject: [PATCH 23/26] added extra becke grid --- src/becke_numerical_grid/EZFIO.cfg | 2 +- .../extra_grid_vector.irp.f | 60 +++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 src/becke_numerical_grid/extra_grid_vector.irp.f diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index 53b11fbd..4083e0e7 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -38,7 +38,7 @@ default: 1202 type: integer doc: Type of extra_grid used for the Becke's numerical extra_grid. Can be, by increasing accuracy: [ 0 | 1 | 2 | 3 ] interface: ezfio,provider,ocaml -default: 2 +default: 0 [thresh_extra_grid] type: double precision diff --git a/src/becke_numerical_grid/extra_grid_vector.irp.f b/src/becke_numerical_grid/extra_grid_vector.irp.f new file mode 100644 index 00000000..3a5e6d3c --- /dev/null +++ b/src/becke_numerical_grid/extra_grid_vector.irp.f @@ -0,0 +1,60 @@ + +BEGIN_PROVIDER [integer, n_points_extra_final_grid] + implicit none + BEGIN_DOC + ! Number of points_extra which are non zero + END_DOC + integer :: i,j,k,l + n_points_extra_final_grid = 0 + do j = 1, nucl_num + do i = 1, n_points_extra_radial_grid -1 + do k = 1, n_points_extra_integration_angular + if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then + cycle + endif + n_points_extra_final_grid += 1 + enddo + enddo + enddo + print*,'n_points_extra_final_grid = ',n_points_extra_final_grid + print*,'n max point = ',n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1) +! call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid) +END_PROVIDER + + BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)] +&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid) ] +&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid) ] +&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ] + implicit none + BEGIN_DOC +! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point +! +! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions +! +! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point +! +! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices + END_DOC + integer :: i,j,k,l,i_count + double precision :: r(3) + i_count = 0 + do j = 1, nucl_num + do i = 1, n_points_extra_radial_grid -1 + do k = 1, n_points_extra_integration_angular + if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then + cycle + endif + i_count += 1 + final_grid_points_extra(1,i_count) = grid_points_extra_per_atom(1,k,i,j) + final_grid_points_extra(2,i_count) = grid_points_extra_per_atom(2,k,i,j) + final_grid_points_extra(3,i_count) = grid_points_extra_per_atom(3,k,i,j) + final_weight_at_r_vector_extra(i_count) = final_weight_at_r_extra(k,i,j) + index_final_points_extra(1,i_count) = k + index_final_points_extra(2,i_count) = i + index_final_points_extra(3,i_count) = j + index_final_points_extra_reverse(k,i,j) = i_count + enddo + enddo + enddo + +END_PROVIDER From ca2b58b495eaa783b9c6857857bced0f811f235d Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 3 Jun 2021 14:57:45 +0200 Subject: [PATCH 24/26] fixed two rdms --- src/dft_utils_in_r/mo_in_r.irp.f | 18 +++++ src/two_body_rdm/example.irp.f | 72 +++++++++++++++++-- src/two_body_rdm/state_av_act_2rdm.irp.f | 6 +- src/two_rdm_routines/update_rdm.irp.f | 30 ++++++-- .../update_state_av_rdm.irp.f | 25 +++++-- src/utils/constants.include.F | 1 + src/utils/linear_algebra.irp.f | 6 +- 7 files changed, 138 insertions(+), 20 deletions(-) diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index d4afcf76..b901ea50 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -89,6 +89,24 @@ enddo END_PROVIDER + BEGIN_PROVIDER[double precision, mos_grad_in_r_array_transp_bis, (n_points_final_grid,mo_num,3)] + implicit none + BEGIN_DOC +! Transposed gradients +! + END_DOC + integer :: i,j,m + do m = 1, 3 + do j = 1, mo_num + do i = 1, n_points_final_grid + mos_grad_in_r_array_transp_bis(i,j,m) = mos_grad_in_r_array(j,i,m) + enddo + enddo + enddo + END_PROVIDER + + + BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)] &BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)] implicit none diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index e9cbd1a2..6cf0209e 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -15,7 +15,7 @@ subroutine routine_active_only double precision :: wee_aa_st_av, rdm_aa_st_av double precision :: wee_bb_st_av, rdm_bb_st_av double precision :: wee_ab_st_av, rdm_ab_st_av - double precision :: wee_tot_st_av, rdm_tot_st_av + double precision :: wee_tot_st_av, rdm_tot_st_av,spin_trace double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 wee_ab = 0.d0 @@ -38,6 +38,27 @@ subroutine routine_active_only provide act_2_rdm_ab_mo act_2_rdm_aa_mo act_2_rdm_bb_mo act_2_rdm_spin_trace_mo provide state_av_act_2_rdm_ab_mo state_av_act_2_rdm_aa_mo provide state_av_act_2_rdm_bb_mo state_av_act_2_rdm_spin_trace_mo + i = 1 + j = 2 +! print*,'testing stuffs' +! istate = 1 +! print*,'alpha/beta' +! print*,'',j,i,j,i +! print*,act_2_rdm_ab_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_ab_mo(i,j,i,j,istate) +! print*,'alpha/alpha' +! print*,'',j,i,j,i +! print*,act_2_rdm_aa_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_aa_mo(i,j,i,j,istate) +! print*,'spin_trace' +! print*,'',j,i,j,i +! print*,act_2_rdm_spin_trace_mo(j,i,j,i,istate) +! print*,'',i,j,i,j +! print*,act_2_rdm_spin_trace_mo(i,j,i,j,istate) +! stop +! print*,'**************************' print*,'**************************' do istate = 1, N_states @@ -51,6 +72,19 @@ subroutine routine_active_only korb = list_act(k) do l = 1, n_act_orb lorb = list_act(l) + if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)).gt.1.d-10)then + print*,'Error in act_2_rdm_spin_trace_mo' + print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l) - act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10" + print*,i,j,k,l + print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(j,i,l,k,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)) + endif + + if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)).gt.1.d-10)then + print*,'Error in act_2_rdm_spin_trace_mo' + print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate),istate).gt.1.d-10" + print*,i,j,k,l + print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(k,l,i,j,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)) + endif vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) @@ -59,7 +93,18 @@ subroutine routine_active_only rdmaa = act_2_rdm_aa_mo(l,k,j,i,istate) rdmbb = act_2_rdm_bb_mo(l,k,j,i,istate) rdmtot = act_2_rdm_spin_trace_mo(l,k,j,i,istate) + spin_trace = rdmaa + rdmbb + rdmab + if(dabs(rdmtot- spin_trace).gt.1.d-10)then + print*,'Error in non state average !!!!' + print*,l,k,j,i + print*,lorb,korb,jorb,iorb + print*,spin_trace,rdmtot,dabs(spin_trace - rdmtot) + print*,'rdmab,rdmaa,rdmbb' + print*, rdmab,rdmaa,rdmbb + + endif + wee_ab(istate) += vijkl * rdmab wee_aa(istate) += vijkl * rdmaa @@ -71,8 +116,8 @@ subroutine routine_active_only enddo enddo wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate) - wee_bb_st_av_2 += wee_aa(istate) * state_average_weight(istate) - wee_ab_st_av_2 += wee_aa(istate) * state_average_weight(istate) + wee_bb_st_av_2 += wee_bb(istate) * state_average_weight(istate) + wee_ab_st_av_2 += wee_ab(istate) * state_average_weight(istate) wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate) wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate) print*,'' @@ -87,7 +132,6 @@ subroutine routine_active_only print*,'Full energy ' print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) enddo - wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 wee_ab_st_av = 0.d0 @@ -103,10 +147,30 @@ subroutine routine_active_only vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10)then + print*,'Error in state_av_act_2_rdm_spin_trace_mo' + print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10" + print*,i,j,k,l + print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(j,i,l,k),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)) + endif + + if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10)then + print*,'Error in state_av_act_2_rdm_spin_trace_mo' + print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10" + print*,i,j,k,l + print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(k,l,i,j),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)) + endif + rdm_aa_st_av = state_av_act_2_rdm_aa_mo(l,k,j,i) rdm_bb_st_av = state_av_act_2_rdm_bb_mo(l,k,j,i) rdm_ab_st_av = state_av_act_2_rdm_ab_mo(l,k,j,i) + spin_trace = rdm_aa_st_av + rdm_bb_st_av + rdm_ab_st_av rdm_tot_st_av = state_av_act_2_rdm_spin_trace_mo(l,k,j,i) + if(dabs(spin_trace - rdm_tot_st_av).gt.1.d-10)then + print*,'Error !!!!' + print*,l,k,j,i + print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av) + endif wee_aa_st_av += vijkl * rdm_aa_st_av wee_bb_st_av += vijkl * rdm_bb_st_av diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index d85c3cdb..db793047 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -119,7 +119,11 @@ call wall_time(wall_1) double precision :: wall_1, wall_2 print*,'providing state_av_act_2_rdm_spin_trace_mo ' - call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) + state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo & + + state_av_act_2_rdm_aa_mo & + + state_av_act_2_rdm_bb_mo + +! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call wall_time(wall_2) print*,'Time to provide state_av_act_2_rdm_spin_trace_mo',wall_2 - wall_1 diff --git a/src/two_rdm_routines/update_rdm.irp.f b/src/two_rdm_routines/update_rdm.irp.f index 662b3ee9..8aeb0699 100644 --- a/src/two_rdm_routines/update_rdm.irp.f +++ b/src/two_rdm_routines/update_rdm.irp.f @@ -61,14 +61,24 @@ ! Therefore you don't necessayr have symmetry between electron 1 and 2 nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) + values(istate,nkeys) = 0.5d0 * c_1(istate) enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = h1 keys(4,nkeys) = h2 + + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 enddo enddo + else if (alpha_alpha)then do i = 1, n_occ_ab(1) i1 = occ(i,1) @@ -259,12 +269,20 @@ if(alpha_beta)then nkeys += 1 do istate = 1, N_st - values(istate,nkeys) = c_1(istate) * phase + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = p2 + nkeys += 1 + do istate = 1, N_st + values(istate,nkeys) = 0.5d0 * c_1(istate) * phase + enddo + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 else if(spin_trace)then nkeys += 1 do istate = 1, N_st @@ -278,10 +296,10 @@ do istate = 1, N_st values(istate,nkeys) = 0.5d0 * c_1(istate) * phase enddo - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 endif end diff --git a/src/two_rdm_routines/update_state_av_rdm.irp.f b/src/two_rdm_routines/update_state_av_rdm.irp.f index bf7b6148..0038e94b 100644 --- a/src/two_rdm_routines/update_state_av_rdm.irp.f +++ b/src/two_rdm_routines/update_state_av_rdm.irp.f @@ -60,11 +60,18 @@ ! If alpha/beta, electron 1 is alpha, electron 2 is beta ! Therefore you don't necessayr have symmetry between electron 1 and 2 nkeys += 1 - values(nkeys) = c_1 + values(nkeys) = 0.5d0 * c_1 keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = h1 keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = h1 enddo enddo else if (alpha_alpha)then @@ -236,11 +243,17 @@ p2 = list_orb_reverse(p2) if(alpha_beta)then nkeys += 1 - values(nkeys) = c_1 * phase + values(nkeys) = 0.5d0 * c_1 * phase keys(1,nkeys) = h1 keys(2,nkeys) = h2 keys(3,nkeys) = p1 keys(4,nkeys) = p2 + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 else if(spin_trace)then nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase @@ -250,10 +263,10 @@ keys(4,nkeys) = p2 nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 endif end diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index 7399b4a6..297a839e 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -3,6 +3,7 @@ integer, parameter :: SIMD_vector = 32 integer, parameter :: N_int_max = 32 double precision, parameter :: pi = dacos(-1.d0) +double precision, parameter :: inv_pi = 1.d0/dacos(-1.d0) double precision, parameter :: sqpi = dsqrt(dacos(-1.d0)) double precision, parameter :: pi_5_2 = 34.9868366552d0 double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 5b2f9067..405d2d20 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1589,9 +1589,9 @@ subroutine restore_symmetry(m,n,A,LDA,thresh) thresh2 = dsqrt(thresh) call nullify_small_elements(m,n,A,LDA,thresh) - if (.not.restore_symm) then - return - endif +! if (.not.restore_symm) then +! return +! endif ! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)): ! - copy all values in a 1D array From 0bf0513fb13623d3a2807c11a8316eea51a72705 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 11 Jun 2021 14:18:23 +0200 Subject: [PATCH 25/26] Added VERSION file --- VERSION | 1 + ocaml/qp_create_ezfio.ml | 27 ++++++++++++++++++++++++--- src/ao_basis/aos.irp.f | 10 +++++++--- src/basis/EZFIO.cfg | 10 ++++++++-- 4 files changed, 40 insertions(+), 8 deletions(-) create mode 100644 VERSION diff --git a/VERSION b/VERSION new file mode 100644 index 00000000..c043eea7 --- /dev/null +++ b/VERSION @@ -0,0 +1 @@ +2.2.1 diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index e4ead21b..a4865e2b 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -604,9 +604,30 @@ let run ?o b au c d m p cart xyz_file = ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ; Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ; - Ezfio.set_basis_shell_nucl (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| shell_num |] - ~data:(list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ; + Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] + ~data:( + list_map (fun (_,n) -> Nucl_number.to_int n) basis + |> List.fold_left (fun accu i -> + match accu with + | [] -> [] + | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((h+1,i)::(h+1,j)::rest) + ) [(0,0)] + |> List.rev + |> List.map fst + )) ; + Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] + ~data:( + list_map (fun (_,n) -> Nucl_number.to_int n) basis + |> List.fold_left (fun accu i -> + match accu with + | [] -> [(1,i)] + | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest) + ) [] + |> List.rev + |> List.map fst + )) ; Ezfio.set_basis_prim_coef (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| prim_num |] ~data:coef) ; Ezfio.set_basis_prim_expo (Ezfio.ezfio_array_of_list diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 2ff72898..17aa784d 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -40,9 +40,9 @@ END_PROVIDER do i=1,ao_num - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) + powA(1) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) + powA(2) = 0 + powA(3) = 0 ! Normalization of the primitives if (primitives_normalized) then @@ -57,6 +57,10 @@ END_PROVIDER enddo endif + powA(1) = ao_power(i,1) + powA(2) = ao_power(i,2) + powA(3) = ao_power(i,3) + ! Normalization of the contracted basis functions if (ao_normalized) then norm = 0.d0 diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index 42123373..7f2ede4c 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -13,6 +13,12 @@ type: integer doc: Number of shells interface: ezfio, provider +[nucleus_shell_num] +type: integer +doc: Number of shells per nucleus +size: (nuclei.nucl_num) +interface: ezfio, provider + [shell_normalization_factor] type: double precision doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ }$ @@ -37,10 +43,10 @@ doc: Max number of primitives in a shell size: (basis.shell_num) interface: ezfio, provider -[shell_nucl] +[basis_nucleus_index] type: integer doc: Index of the nucleus on which the shell is centered -size: (basis.shell_num) +size: (nuclei.nucl_num) interface: ezfio, provider [prim_normalization_factor] From b75ee98597a8769b050c6709a1c0f075f10a47f1 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Fri, 11 Jun 2021 16:19:32 +0200 Subject: [PATCH 26/26] added other transposed gradients of AOs --- src/dft_utils_in_r/ao_in_r.irp.f | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index a98afd5d..6fa6a4c7 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -139,3 +139,22 @@ END_PROVIDER + BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_3, (3,n_points_final_grid,ao_num)] + implicit none + BEGIN_DOC +! Transposed gradients +! + END_DOC + integer :: i,j,m + double precision :: aos_array(ao_num), r(3) + double precision :: aos_grad_array(3,ao_num) + do m = 1, 3 + do j = 1, ao_num + do i = 1, n_points_final_grid + aos_grad_in_r_array_transp_3(m,i,j) = aos_grad_in_r_array(j,i,m) + enddo + enddo + enddo + END_PROVIDER + +