From e9962280dc101247c72a505c0a85de06682bc2ba Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Tue, 27 Apr 2021 03:32:12 +0200 Subject: [PATCH] save_cluster --- src/JASTROW/jastrow_simple.irp.f | 26 ++ src/PROPERTIES/properties_ci_SVD.irp.f | 221 ++++++++++++ src/PROPERTIES/properties_ci_postSVD.irp.f | 235 +++++++++++++ src/PROPERTIES/properties_deriv_jast.irp.f | 124 +++++++ src/PROPERTIES/properties_energy.irp.f | 1 - src/PROPERTIES/properties_energy_SVD.irp.f | 49 +++ src/det.irp.f | 370 +++++++++++++++++++++ src/ezfio_interface.irp.f | 5 +- src/psi.irp.f | 101 ++++-- src/psi_SVD.irp.f | 85 +++++ src/wf.irp.f | 6 +- 11 files changed, 1199 insertions(+), 24 deletions(-) create mode 100644 src/PROPERTIES/properties_ci_SVD.irp.f create mode 100644 src/PROPERTIES/properties_ci_postSVD.irp.f create mode 100644 src/PROPERTIES/properties_deriv_jast.irp.f create mode 100644 src/PROPERTIES/properties_energy_SVD.irp.f create mode 100644 src/psi_SVD.irp.f diff --git a/src/JASTROW/jastrow_simple.irp.f b/src/JASTROW/jastrow_simple.irp.f index bbd9c99..6bd471f 100644 --- a/src/JASTROW/jastrow_simple.irp.f +++ b/src/JASTROW/jastrow_simple.irp.f @@ -192,3 +192,29 @@ BEGIN_PROVIDER [ double precision , jast_elec_Simple_lapl, (elec_num_8) ] END_PROVIDER + +BEGIN_PROVIDER [ double precision , jast_elec_Simple_deriv_nucPar, (nucl_num) ] + implicit none + BEGIN_DOC +! Variation of the Jastrow factor with respect to nuclear parameters + END_DOC + + integer :: i,j + double precision :: a, rij, tmp + + do j=1,nucl_num + jast_elec_Simple_deriv_nucPar(j) = 0.d0 + enddo + do j=1,nucl_num + a = jast_pen(j) + tmp = 0.d0 + !DIR$ LOOP COUNT (100) + do i=1,elec_num + rij = nucl_elec_dist(i,j) + tmp += rij*rij/((1.d0+a*rij)*(1.d0+a*rij)*(1.d0+a*rij)) + enddo + jast_elec_Simple_deriv_nucPar(j) = -2.d0*a*tmp + enddo + + +END_PROVIDER diff --git a/src/PROPERTIES/properties_ci_SVD.irp.f b/src/PROPERTIES/properties_ci_SVD.irp.f new file mode 100644 index 0000000..c2c822f --- /dev/null +++ b/src/PROPERTIES/properties_ci_SVD.irp.f @@ -0,0 +1,221 @@ + + + +BEGIN_PROVIDER [ double precision, ci_overlap_psidet_SVD, (size_ci_overlap_psidet_SVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < psi_0 | det(j) > + ! Dimensions : n_svd_coefs + END_DOC + + integer :: k + do k = 1, n_svd_coefs + ci_overlap_psidet_SVD(k) = det_alpha_value_SVD(k) * det_beta_value_SVD(k) * psidet_inv_SVD + enddo + + ci_overlap_psidet_SVD_min = min(ci_overlap_psidet_SVD_min,minval(ci_overlap_psidet_SVD)) + ci_overlap_psidet_SVD_max = max(ci_overlap_psidet_SVD_max,maxval(ci_overlap_psidet_SVD)) + SOFT_TOUCH ci_overlap_psidet_SVD_min ci_overlap_psidet_SVD_max +END_PROVIDER + + + + + +BEGIN_PROVIDER [ double precision, ci_h_psidet_SVD, (size_ci_h_psidet_SVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < psi_0 |H| det(j) > + ! Dimensions : n_svd_coefs + END_DOC + + integer :: k, e + double precision :: T + + do k = 1, n_svd_coefs + T = 0.d0 + do e = 1, elec_alpha_num + T += det_alpha_grad_lapl_SVD(4,e,k) * det_beta_value_SVD(k) + enddo + do e = elec_beta_num+1, elec_num + T += det_beta_grad_lapl_SVD(4,e,k) * det_alpha_value_SVD(k) + enddo + ci_h_psidet_SVD(k) = -0.5d0*T + E_pot * det_alpha_value_SVD(k) * det_beta_value_SVD(k) + ci_h_psidet_SVD(k) *= psidet_inv_SVD + enddo + + ci_h_psidet_SVD_min = min(ci_h_psidet_SVD_min,minval(ci_h_psidet_SVD)) + ci_h_psidet_SVD_max = max(ci_h_psidet_SVD_max,maxval(ci_h_psidet_SVD)) + SOFT_TOUCH ci_h_psidet_SVD_min ci_h_psidet_SVD_max +END_PROVIDER + + + + + + + +BEGIN_PROVIDER [ double precision, ci_overlap_matrix_SVD, (size_ci_overlap_matrix_SVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < det(i) | det(j) > + ! Dimensions : n_svd_coefs * n_svd_coefs + END_DOC + + integer :: k, l + double precision :: f + + do k = 1, n_svd_coefs + f = det_alpha_value_SVD(k) * det_beta_value_SVD(k) * psidet_inv_SVD * psidet_inv_SVD + do l = 1, n_svd_coefs + ci_overlap_matrix_SVD( n_svd_coefs*(k-1) + l) = det_alpha_value_SVD(l) * det_beta_value_SVD(l) * f + enddo + enddo + + ci_overlap_matrix_SVD_min = min(ci_overlap_matrix_SVD_min,minval(ci_overlap_matrix_SVD)) + ci_overlap_matrix_SVD_max = max(ci_overlap_matrix_SVD_max,maxval(ci_overlap_matrix_SVD)) + SOFT_TOUCH ci_overlap_matrix_SVD_min ci_overlap_matrix_SVD_max +END_PROVIDER + + + + + + +BEGIN_PROVIDER [ double precision, ci_h_matrix_SVD, (size_ci_h_matrix_SVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < det(i) | H | det(j) > + ! Dimensions : n_svd_coefs * n_svd_coefs + END_DOC + + integer :: k, l, e + double precision :: f, g, h, T, V + + do l = 1, n_svd_coefs + ! Lapl D + g = 0.d0 + do e = 1, elec_alpha_num + g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(l) + enddo + do e = elec_alpha_num+1, elec_num + g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,l) + enddo + T = g + ! D (Lapl J)/J + g = 0.d0 + do e = 1, elec_num + g += jast_lapl_jast_inv(e) + enddo + T += det_alpha_value_SVD(l) * det_beta_value_SVD(l) * g + ! 2 (grad D).(Grad J)/J + g = 0.d0 + do e = 1, elec_alpha_num + g += & + det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + h = 0.d0 + do e = elec_alpha_num+1, elec_num + h += & + det_beta_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_beta_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_beta_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + T += 2.d0 * ( g * det_beta_value_SVD(l) + h * det_alpha_value_SVD(l) ) + g = det_alpha_value_SVD(l) * det_beta_value_SVD(l) + V = E_pot * g + do e = 1, elec_alpha_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(l) + enddo + do e = elec_alpha_num+1, elec_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,l) + enddo + f = -0.5d0*T + V + f *= psidet_inv_SVD * psidet_inv_SVD + do k = 1, n_svd_coefs + ci_h_matrix_SVD( n_svd_coefs*(l-1) + k) = f * det_alpha_value_SVD(k) * det_beta_value_SVD(k) + enddo + enddo + + ci_h_matrix_SVD_min = min(ci_h_matrix_SVD_min,minval(ci_h_matrix_SVD)) + ci_h_matrix_SVD_max = max(ci_h_matrix_SVD_max,maxval(ci_h_matrix_SVD)) + SOFT_TOUCH ci_h_matrix_SVD_min ci_h_matrix_SVD_max +END_PROVIDER + + + + + + +BEGIN_PROVIDER [ double precision, ci_h_matrix_diag_SVD, (size_ci_h_matrix_diag_SVD) ] + implicit none + BEGIN_DOC + ! < det(i) |H| det(j) > + ! + ! Dimensions : n_svd_coefs + END_DOC + + integer :: l, e + double precision :: f, g, h, T, V + + do l = 1, n_svd_coefs + ! Lapl D + g = 0.d0 + do e = 1, elec_alpha_num + g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(l) + enddo + do e = elec_alpha_num+1, elec_num + g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,l) + enddo + T = g + ! D (Lapl J)/J + g = 0.d0 + do e = 1, elec_num + g += jast_lapl_jast_inv(e) + enddo + T += det_alpha_value(l) * det_beta_value_SVD(l) * g + ! 2 (grad D).(Grad J)/J + g = 0.d0 + do e = 1 , elec_alpha_num + g += & + det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + h = 0.d0 + do e = elec_alpha_num+1, elec_num + h += & + det_beta_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_beta_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_beta_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + T += 2.d0 * ( g * det_beta_value_SVD(l) + h * det_alpha_value_SVD(l) ) + g = det_alpha_value_SVD(l) * det_beta_value_SVD(l) + V = E_pot * g + do e = 1, elec_alpha_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(l) + enddo + do e = elec_alpha_num+1, elec_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,l) + enddo + f = -0.5d0*T + V + f *= psidet_inv_SVD * psidet_inv_SVD + ci_h_matrix_diag_SVD(l) = f * det_alpha_value_SVD(l) * det_beta_value_SVD(l) + enddo + + ci_h_matrix_diag_SVD_min = min(ci_h_matrix_diag_SVD_min,minval(ci_h_matrix_diag_SVD)) + ci_h_matrix_diag_SVD_max = max(ci_h_matrix_diag_SVD_max,maxval(ci_h_matrix_diag_SVD)) + SOFT_TOUCH ci_h_matrix_diag_SVD_min ci_h_matrix_diag_SVD_max +END_PROVIDER + + diff --git a/src/PROPERTIES/properties_ci_postSVD.irp.f b/src/PROPERTIES/properties_ci_postSVD.irp.f new file mode 100644 index 0000000..8222a1b --- /dev/null +++ b/src/PROPERTIES/properties_ci_postSVD.irp.f @@ -0,0 +1,235 @@ + + + +BEGIN_PROVIDER [ double precision, ci_overlap_psidet_postSVD, (size_ci_overlap_psidet_postSVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < psi_0 | det(j) > + ! Dimensions : n_svd_coefs2 + END_DOC + + integer :: k, kp + do k = 1, n_svd_coefs + do kp = 1, n_svd_coefs + ci_overlap_psidet_postSVD(kp+(k-1)*n_svd_coefs) = det_alpha_value_SVD(k) * det_beta_value_SVD(kp) * psidet_inv_SVD + enddo + enddo + + ci_overlap_psidet_postSVD_min = min(ci_overlap_psidet_postSVD_min,minval(ci_overlap_psidet_postSVD)) + ci_overlap_psidet_postSVD_max = max(ci_overlap_psidet_postSVD_max,maxval(ci_overlap_psidet_postSVD)) + SOFT_TOUCH ci_overlap_psidet_postSVD_min ci_overlap_psidet_postSVD_max +END_PROVIDER + + + + + +BEGIN_PROVIDER [ double precision, ci_h_psidet_postSVD, (size_ci_h_psidet_postSVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < psi_0 |H| det(j) > + ! Dimensions : n_svd_coefs2 + END_DOC + + integer :: k, kp, e + double precision :: T + + do k = 1, n_svd_coefs + do kp = 1, n_svd_coefs + T = 0.d0 + do e = 1, elec_alpha_num + T += det_alpha_grad_lapl_SVD(4,e,k) * det_beta_value_SVD(kp) + enddo + do e = elec_beta_num+1, elec_num + T += det_alpha_value_SVD(k) * det_beta_grad_lapl_SVD(4,e,kp) + enddo + ci_h_psidet_postSVD(kp+(k-1)*n_svd_coefs) = -0.5d0*T + E_pot * det_alpha_value_SVD(k) * det_beta_value_SVD(kp) + ci_h_psidet_postSVD(kp+(k-1)*n_svd_coefs) *= psidet_inv_SVD + enddo + enddo + + ci_h_psidet_postSVD_min = min(ci_h_psidet_postSVD_min,minval(ci_h_psidet_postSVD)) + ci_h_psidet_postSVD_max = max(ci_h_psidet_postSVD_max,maxval(ci_h_psidet_postSVD)) + SOFT_TOUCH ci_h_psidet_postSVD_min ci_h_psidet_postSVD_max +END_PROVIDER + + + + + + + +BEGIN_PROVIDER [ double precision, ci_overlap_matrix_postSVD, (size_ci_overlap_matrix_postSVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < det(i) | det(j) > + ! Dimensions : n_svd_coefs2 * n_svd_coefs2 + END_DOC + + integer :: k, kp, l, lp + double precision :: f + + do k = 1, n_svd_coefs + do kp = 1, n_svd_coefs + f = det_alpha_value_SVD(k) * det_beta_value_SVD(kp) * psidet_inv_SVD * psidet_inv_SVD + do l = 1, n_svd_coefs + do lp = 1, n_svd_coefs + ci_overlap_matrix_postSVD(lp+(l-1)*n_svd_coefs+(kp-1)*n_svd_coefs2+(k-1)*n_svd_coefs3) = det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * f + enddo + enddo + enddo + enddo + + ci_overlap_matrix_postSVD_min = min(ci_overlap_matrix_postSVD_min,minval(ci_overlap_matrix_postSVD)) + ci_overlap_matrix_postSVD_max = max(ci_overlap_matrix_postSVD_max,maxval(ci_overlap_matrix_postSVD)) + SOFT_TOUCH ci_overlap_matrix_postSVD_min ci_overlap_matrix_postSVD_max +END_PROVIDER + + + + + + +BEGIN_PROVIDER [ double precision, ci_h_matrix_postSVD, (size_ci_h_matrix_postSVD) ] + implicit none + BEGIN_DOC + ! !!! + ! < det(i) | H | det(j) > + ! Dimensions : n_svd_coefs2 * n_svd_coefs2 + END_DOC + + integer :: k, kp, l, lp, e + double precision :: f, g, h, T, V + + do l = 1, n_svd_coefs + do lp = 1, n_svd_coefs + ! Lapl D + g = 0.d0 + do e = 1, elec_alpha_num + g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(lp) + enddo + do e = elec_alpha_num+1, elec_num + g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,lp) + enddo + T = g + ! D (Lapl J)/J + g = 0.d0 + do e = 1, elec_num + g += jast_lapl_jast_inv(e) + enddo + T += det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * g + ! 2 (grad D).(Grad J)/J + g = 0.d0 + do e = 1, elec_alpha_num + g += & + det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + h = 0.d0 + do e = elec_alpha_num+1, elec_num + h += & + det_beta_grad_lapl_SVD(1,e,lp) * jast_grad_jast_inv_x(e) + & + det_beta_grad_lapl_SVD(2,e,lp) * jast_grad_jast_inv_y(e) + & + det_beta_grad_lapl_SVD(3,e,lp) * jast_grad_jast_inv_z(e) + enddo + T += 2.d0 * ( g * det_beta_value_SVD(lp) + h * det_alpha_value_SVD(l) ) + g = det_alpha_value_SVD(l) * det_beta_value_SVD(lp) + V = E_pot * g + do e = 1, elec_alpha_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(lp) + enddo + do e = elec_alpha_num+1, elec_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,lp) + enddo + f = -0.5d0*T + V + f *= psidet_inv_SVD * psidet_inv_SVD + do k = 1, n_svd_coefs + do kp = 1, n_svd_coefs + ci_h_matrix_postSVD(kp+(k-1)*n_svd_coefs+(lp-1)*n_svd_coefs2+(l-1)*n_svd_coefs3) = f * det_alpha_value_SVD(k) * det_beta_value_SVD(kp) + enddo + enddo + enddo + enddo + + ci_h_matrix_postSVD_min = min(ci_h_matrix_postSVD_min,minval(ci_h_matrix_postSVD)) + ci_h_matrix_postSVD_max = max(ci_h_matrix_postSVD_max,maxval(ci_h_matrix_postSVD)) + SOFT_TOUCH ci_h_matrix_postSVD_min ci_h_matrix_postSVD_max +END_PROVIDER + + + + + + +BEGIN_PROVIDER [ double precision, ci_h_matrix_diag_postSVD, (size_ci_h_matrix_diag_postSVD) ] + implicit none + BEGIN_DOC + ! < det(i) |H| det(j) > + ! + ! Dimensions : n_svd_coefs2 + END_DOC + + integer :: l, lp, e + double precision :: f, g, h, T, V + + do l = 1, n_svd_coefs + do lp = 1, n_svd_coefs + ! Lapl D + g = 0.d0 + do e = 1, elec_alpha_num + g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(lp) + enddo + do e = elec_alpha_num+1, elec_num + g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,lp) + enddo + T = g + ! D (Lapl J)/J + g = 0.d0 + do e = 1, elec_num + g += jast_lapl_jast_inv(e) + enddo + T += det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * g + ! 2 (grad D).(Grad J)/J + g = 0.d0 + do e = 1, elec_alpha_num + g += & + det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + & + det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + & + det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e) + enddo + h = 0.d0 + do e = elec_alpha_num+1, elec_num + h += & + det_beta_grad_lapl_SVD(1,e,lp) * jast_grad_jast_inv_x(e) + & + det_beta_grad_lapl_SVD(2,e,lp) * jast_grad_jast_inv_y(e) + & + det_beta_grad_lapl_SVD(3,e,lp) * jast_grad_jast_inv_z(e) + enddo + T += 2.d0 * ( g * det_beta_value_SVD(lp) + h * det_alpha_value_SVD(l) ) + g = det_alpha_value_SVD(l) * det_beta_value_SVD(lp) + V = E_pot * g + do e = 1, elec_alpha_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(lp) + enddo + do e = elec_alpha_num+1, elec_num + V -= pseudo_non_local_SVD(e) * g + V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,lp) + enddo + f = -0.5d0*T + V + f *= psidet_inv_SVD * psidet_inv_SVD + ci_h_matrix_diag_postSVD(lp+(l-1)*n_svd_coefs) = f * det_alpha_value_SVD(l) * det_beta_value_SVD(lp) + enddo + enddo + + ci_h_matrix_diag_postSVD_min = min(ci_h_matrix_diag_postSVD_min,minval(ci_h_matrix_diag_postSVD)) + ci_h_matrix_diag_postSVD_max = max(ci_h_matrix_diag_postSVD_max,maxval(ci_h_matrix_diag_postSVD)) + SOFT_TOUCH ci_h_matrix_diag_postSVD_min ci_h_matrix_diag_postSVD_max +END_PROVIDER + + diff --git a/src/PROPERTIES/properties_deriv_jast.irp.f b/src/PROPERTIES/properties_deriv_jast.irp.f new file mode 100644 index 0000000..b785067 --- /dev/null +++ b/src/PROPERTIES/properties_deriv_jast.irp.f @@ -0,0 +1,124 @@ + +BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc1, (size_E_deriv_nucPar_loc1) ] + implicit none + BEGIN_DOC + ! Local energy variation with respect to nuclear parameters + ! + ! Dimensions : nucl_num + END_DOC + + integer :: i + + do i=1,nucl_num + E_deriv_nucPar_loc1(i) = E_loc*jast_elec_Simple_deriv_nucPar(i) + enddo + + E_deriv_nucPar_loc1_min = min(E_deriv_nucPar_loc1_min,minval(E_deriv_nucPar_loc1)) + E_deriv_nucPar_loc1_max = max(E_deriv_nucPar_loc1_max,maxval(E_deriv_nucPar_loc1)) + SOFT_TOUCH E_deriv_nucPar_loc1_min E_deriv_nucPar_loc1_max +END_PROVIDER + +BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc2, (size_E_deriv_nucPar_loc2) ] + implicit none + BEGIN_DOC + ! Local energy variation with respect to nuclear parameters + ! + ! Dimensions : nucl_num + END_DOC + + integer :: i + + do i=1,nucl_num + E_deriv_nucPar_loc2(i) = jast_elec_Simple_deriv_nucPar(i) + enddo + + E_deriv_nucPar_loc2_min = min(E_deriv_nucPar_loc2_min,minval(E_deriv_nucPar_loc2)) + E_deriv_nucPar_loc2_max = max(E_deriv_nucPar_loc2_max,maxval(E_deriv_nucPar_loc2)) + SOFT_TOUCH E_deriv_nucPar_loc2_min E_deriv_nucPar_loc2_max +END_PROVIDER + +BEGIN_PROVIDER [ double precision, J_deriv_nucPar_verif, (size_J_deriv_nucPar_verif) ] + implicit none + BEGIN_DOC + ! Jastrow variation with respect to nuclear parameters + ! + ! Dimensions : nucl_num + END_DOC + + integer :: i,j + double precision :: sqt_eps = 1d-3, der1, der2, pos1, pos2, pos0 + + do j=1,nucl_num + !! + pos0 = sqt_eps * jast_pen(j) + !! + jast_pen(j) = jast_pen(j) + pos0 + TOUCH jast_pen + pos1 = jast_pen(j) + !DIR$ LOOP COUNT (100) + do i=1,elec_num + der1 += jast_elec_Simple_value(i) + enddo + !! + jast_pen(j) = jast_pen(j) - 2.d0 * pos0 + TOUCH jast_pen + pos2 = jast_pen(j) + !DIR$ LOOP COUNT (100) + do i=1,elec_num + der2 += jast_elec_Simple_value(i) + enddo + !! + J_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0 + !! + enddo + + J_deriv_nucPar_verif_min = min(J_deriv_nucPar_verif_min,minval(J_deriv_nucPar_verif)) + J_deriv_nucPar_verif_max = max(J_deriv_nucPar_verif_max,maxval(J_deriv_nucPar_verif)) + SOFT_TOUCH J_deriv_nucPar_verif_min J_deriv_nucPar_verif_max +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, E_deriv_nucPar_verif, (size_E_deriv_nucPar_verif) ] + implicit none + BEGIN_DOC + ! Local energy variation with respect to nuclear parameters: verification + ! + ! Dimensions : nucl_num + END_DOC + + integer :: j + double precision :: sqt_eps = 1d-3, der1, der2, pos1, pos2, pos0 + + do j=1,nucl_num + !! + pos0 = sqt_eps * jast_pen(j) + !! + jast_pen(j) = jast_pen(j) + pos0 + TOUCH jast_pen + pos1 = jast_pen(j) + der1 = E_loc + !!DIR$ LOOP COUNT (100) + !do i=1,elec_num + ! der1 += jast_elec_Simple_value(i) + !enddo + !! + jast_pen(j) = jast_pen(j) - 2.d0 * pos0 + TOUCH jast_pen + pos2 = jast_pen(j) + der2 = E_loc + !!DIR$ LOOP COUNT (100) + !do i=1,elec_num + ! der2 += jast_elec_Simple_value(i) + !enddo + !! + E_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0 + !! + enddo + + E_deriv_nucPar_verif_min = min(E_deriv_nucPar_verif_min,minval(E_deriv_nucPar_verif)) + E_deriv_nucPar_verif_max = max(E_deriv_nucPar_verif_max,maxval(E_deriv_nucPar_verif)) + SOFT_TOUCH E_deriv_nucPar_verif_min E_deriv_nucPar_verif_max + +END_PROVIDER + diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 0f01449..95ca9e3 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -281,4 +281,3 @@ END_PROVIDER - diff --git a/src/PROPERTIES/properties_energy_SVD.irp.f b/src/PROPERTIES/properties_energy_SVD.irp.f new file mode 100644 index 0000000..0c3c9eb --- /dev/null +++ b/src/PROPERTIES/properties_energy_SVD.irp.f @@ -0,0 +1,49 @@ +!==========================================================================! +! DIMENSIONS ! +!==========================================================================! + +BEGIN_PROVIDER [ double precision, E_kin_elec_SVD, (elec_num) ] + implicit none + BEGIN_DOC + ! Electronic Kinetic energy : -1/2 (Lapl.Psi_SVD)/Psi_SVD + END_DOC + integer :: i + do i = 1, elec_num + E_kin_elec_SVD(i) = -0.5d0 * psi_lapl_psi_inv_SVD(i) + enddo +END_PROVIDER + + + +!==========================================================================! +! PROPERTIES ! +!==========================================================================! + + +BEGIN_PROVIDER [ double precision, E_loc_SVD ] + implicit none + include '../types.F' + BEGIN_DOC + ! Local energy : E_kin + E_pot + E_nucl + END_DOC + + integer :: i + + E_loc_SVD = E_nucl + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT(200) + do i = 1, elec_num + E_loc_SVD += E_kin_elec_SVD(i) + E_pot_elec(i) + enddo + + ! Avoid divergence of E_loc_SVD and population explosion + if (do_pseudo) then + double precision :: delta_e + E_loc_SVD = max(2.d0*E_ref, E_loc_SVD) + endif + + E_loc_SVD_min = min(E_loc_SVD,E_loc_SVD_min) + E_loc_SVD_max = max(E_loc_SVD,E_loc_SVD_max) + SOFT_TOUCH E_loc_SVD_min E_loc_SVD_max + +END_PROVIDER diff --git a/src/det.irp.f b/src/det.irp.f index a7c0fca..69e529c 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -1561,8 +1561,10 @@ END_PROVIDER else + ! DaC(det_beta_num) = Trans( det_coef_matrix_dense(det_alpha_num,det_beta_num) ) @ det_alpha_value(det_alpha_num) call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1) + ! CDb(det_alpha_num) = det_coef_matrix_dense(det_alpha_num,det_beta_num) @ det_beta_value(det_beta_num) call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, & size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1) @@ -1585,6 +1587,7 @@ END_PROVIDER ! Gradients ! --------- + ! psidet_grad_lapl(4,elec_alpha_num) = det_alpha_grad_lapl(4,elec_alpha_num,det_alpha_num) @ CDb(det_alpha_num) call dgemv('N',elec_alpha_num*4,det_alpha_num,1.d0, & det_alpha_grad_lapl, & size(det_alpha_grad_lapl,1)*size(det_alpha_grad_lapl,2), & @@ -1609,6 +1612,367 @@ END_PROVIDER END_PROVIDER + + + + + + + + + + + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ + + BEGIN_PROVIDER [ logical, utilise_SVD ] +&BEGIN_PROVIDER [ integer, n_svd_coefs ] + implicit none + BEGIN_DOC + ! SVD rank + END_DOC + n_svd_coefs = -1 + call get_spindeterminants_n_svd_coefs(n_svd_coefs) + utilise_SVD = n_svd_coefs > 0 + if (.not.utilise_SVD) then + n_svd_coefs = 1 + endif +END_PROVIDER + + BEGIN_PROVIDER [ integer, n_svd_coefs2 ] +&BEGIN_PROVIDER [ integer, n_svd_coefs3 ] + implicit none + BEGIN_DOC + ! square and cube of n_svd_coefs + END_DOC + n_svd_coefs2 = n_svd_coefs * n_svd_coefs + n_svd_coefs3 = n_svd_coefs * n_svd_coefs * n_svd_coefs +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, psi_svd_alpha, (det_alpha_num, n_svd_coefs, n_states) ] +&BEGIN_PROVIDER [ double precision, psi_svd_coefs, (n_svd_coefs, n_states) ] +&BEGIN_PROVIDER [ double precision, psi_svd_beta, (det_beta_num, n_svd_coefs, n_states) ] + implicit none + BEGIN_DOC + ! SVD U + ! SVD sigma + ! SVD Vt + END_DOC + call get_spindeterminants_psi_svd_alpha(psi_svd_alpha) + call get_spindeterminants_psi_svd_coefs(psi_svd_coefs) + call get_spindeterminants_psi_svd_beta(psi_svd_beta) +END_PROVIDER + + + +!!!!!!!!!!!!!!!! +! ! ! ! ! ! ! ! +!!!!!!!!!!!!!!!! + + + + + BEGIN_PROVIDER [ double precision, det_alpha_value_SVD, (n_svd_coefs) ] +&BEGIN_PROVIDER [ double precision, det_beta_value_SVD , (n_svd_coefs) ] + + implicit none + BEGIN_DOC + ! !!! + ! D_alpha_SVD_{k} = sum_i U_{i,k} D_alpha_{i} + ! D_beta_SVD_{k } = sum_j V_{j,k} D_beta_{j } + ! !!! + END_DOC + + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + det_alpha_value_SVD = 0.d0 + det_beta_value_SVD = 0.d0 + endif + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + ! det_alpha_value_SVD = psi_svd_alpha.T @ det_alpha_value + call dgemv('T', det_alpha_num, n_svd_coefs & + , 1.d0, psi_svd_alpha(:,:,1), size(psi_svd_alpha,1), det_alpha_value, 1 & + , 0.d0, det_alpha_value_SVD, 1) + ! !!! + ! det_beta_value_SVD = psi_svd_beta.T @ det_beta_value + call dgemv('T', det_beta_num, n_svd_coefs & + , 1.d0, psi_svd_beta(:,:,1), size(psi_svd_beta,1), det_beta_value, 1 & + , 0.d0, det_beta_value_SVD, 1) + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + + + !do l = 1, n_svd_coefs + ! tmp = 0.d0 + ! do ii = 1, det_alpha_num + ! tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_value(ii) + ! enddo + ! det_alpha_value_SVD(l) = tmp + ! tmp = 0.d0 + ! do jj = 1, det_beta_num + ! tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_value(jj) + ! enddo + ! det_beta_value_SVD(l) = tmp + !enddo + + +END_PROVIDER + + + + + BEGIN_PROVIDER [ double precision, psidet_value_SVD ] +&BEGIN_PROVIDER [ double precision, psidet_inv_SVD ] +&BEGIN_PROVIDER [ double precision, det_alpha_pseudo_SVD, (elec_alpha_num, n_svd_coefs) ] +&BEGIN_PROVIDER [ double precision, det_beta_pseudo_SVD, (elec_alpha_num+1:elec_num, n_svd_coefs) ] + + implicit none + BEGIN_DOC + ! !!! + END_DOC + + integer :: l + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + det_alpha_pseudo_SVD = 0.d0 + det_beta_pseudo_SVD = 0.d0 + endif + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + psidet_value_SVD = 0.d0 + do l = 1, n_svd_coefs + psidet_value_SVD = psidet_value_SVD + det_beta_value_SVD(l) * psi_svd_coefs(l,1) * det_alpha_value_SVD(l) + enddo + if (psidet_value_SVD == 0.d0) then + call abrt(irp_here,'Determinantal component of the SVD wave function is zero.') + endif + psidet_inv_SVD = 1.d0 / psidet_value_SVD + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + if (do_pseudo) then + ! det_alpha_pseudo_SVD = det_alpha_pseudo @ psi_svd_alpha / psidet_inv_SVD + call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, 1.d0 & + , det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) & + , 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) ) + ! !!! + if (elec_beta_num /= 0) then + ! det_beta_pseudo_SVD = det_beta_pseudo @ psi_svd_beta / psidet_inv_SVD + call dgemm('N', 'N', elec_beta_num, n_svd_coefs , det_beta_num, 1.d0 & + , det_beta_pseudo(elec_alpha_num+1,1), size(det_beta_pseudo,1), psi_svd_beta(:,:,1), size(psi_svd_beta,1) & + , 0.d0, det_beta_pseudo_SVD(elec_alpha_num+1,1), size(det_beta_pseudo_SVD,1) ) + endif + endif + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + + + + !do l = 1, n_svd_coefs + ! do ee = 1, elec_alpha_num + ! tmp = 0.d0 + ! do ii = 1, det_alpha_num + ! tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_pseudo(ee,ii) + ! enddo + ! det_alpha_pseudo_SVD(ee,l) = tmp + ! enddo + ! do ee = elec_alpha_num+1, elec_num + ! tmp = 0.d0 + ! do jj = 1, det_beta_num + ! tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_pseudo(ee,jj) + ! enddo + ! det_beta_pseudo_SVD(ee,l) = tmp + ! enddo + !enddo + + +END_PROVIDER + + + +!!!!!!!!!!!!!!!! +! ! ! ! ! ! ! ! +!!!!!!!!!!!!!!!! + + + + BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_SVD, (4, elec_alpha_num , n_svd_coefs) ] +&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_SVD , (4, elec_alpha_num+1:elec_num, n_svd_coefs) ] + + implicit none + BEGIN_DOC + ! !!! + ! det_alpha_grad_lapl_SVD_{k} = sum_i U_{i,k} det_alpha_grad_lapl_{i} + ! det_beta_grad_lapl_SVD_{k } = sum_j V_{j,k} det_beta_grad_lapl_{j } + ! !!! + END_DOC + + integer :: mm + integer, save :: ifirst = 0 + if (ifirst == 0) then + ifirst = 1 + det_alpha_grad_lapl_SVD = 0.d0 + det_beta_grad_lapl_SVD = 0.d0 + endif + + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + do mm = 1, 4 + ! !!! + call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, 1.d0 & + , det_alpha_grad_lapl(mm,:,:), size(det_alpha_grad_lapl,2), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) & + , 0.d0, det_alpha_grad_lapl_SVD(mm,:,:), size(det_alpha_grad_lapl_SVD,2) ) + ! !!! + if (elec_beta_num /= 0) then + call dgemm('N', 'N', elec_beta_num, n_svd_coefs , det_beta_num, 1.d0 & + , det_beta_grad_lapl(mm,:,:), size(det_beta_grad_lapl,2), psi_svd_beta(:,:,1), size(psi_svd_beta,1) & + , 0.d0, det_beta_grad_lapl_SVD(mm,:,:), size(det_beta_grad_lapl_SVD,2) ) + endif + ! !!! + enddo + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + + !do l = 1, n_svd_coefs + ! do mm = 1, 4 + ! do ee = 1, elec_alpha_num + ! tmp = 0.d0 + ! do ii = 1, det_alpha_num + ! tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_grad_lapl(mm,ee,ii) + ! enddo + ! det_alpha_grad_lapl_SVD(mm,ee,l) = tmp + ! enddo + ! do ee = elec_alpha_num+1, elec_num + ! tmp = 0.d0 + ! do jj = 1, det_beta_num + ! tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_grad_lapl(mm,ee,jj) + ! enddo + ! det_beta_grad_lapl_SVD(mm,ee,l) = tmp + ! enddo + ! enddo + !enddo + +END_PROVIDER + + + + + BEGIN_PROVIDER [ double precision, psidet_grad_lapl_SVD, (4,elec_num) ] +&BEGIN_PROVIDER [ double precision, pseudo_non_local_SVD, (elec_num) ] + + implicit none + BEGIN_DOC + ! !!! + ! Sigma is diagonal matrix: Sigma = psi_svd_coefs + ! !!! + ! Gradient of the determinantal part of the wave function after SVD + ! Laplacian of determinantal part of the wave function after SVD + ! Non-local component of the pseudopotentials after SVD + ! !!! + ! for each electron: + ! for mm = 1, 2, 3 : grad_x Psi, grad_y Psi, grad_z Psi + ! for mm = 4: first term (only) of Laplacian Psi + ! !!! + END_DOC + + integer :: l, mm, ee + integer, save :: ifirst=0 + double precision :: tmp + if (ifirst == 0) then + ifirst = 1 + psidet_grad_lapl_SVD = 0.d0 + pseudo_non_local_SVD = 0.d0 + endif + + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + do mm = 1, 4 + ! !!! + do ee = 1, elec_alpha_num + tmp = 0.d0 + do l = 1, n_svd_coefs + tmp = tmp + det_alpha_grad_lapl_SVD(mm,ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l) + enddo + psidet_grad_lapl_SVD(mm,ee) = tmp + enddo + ! !!! + if (elec_beta_num /= 0) then + do ee = elec_alpha_num+1, elec_num + tmp = 0.d0 + do l = 1, n_svd_coefs + tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) * det_beta_grad_lapl_SVD(mm,ee,l) + enddo + psidet_grad_lapl_SVD(mm,ee) = tmp + enddo + endif + ! !!! + enddo + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + + + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + ! !!! + if (do_pseudo) then + ! !!! + do ee = 1, elec_alpha_num + tmp = 0.d0 + do l = 1, n_svd_coefs + tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l) + enddo + pseudo_non_local_SVD(ee) = tmp * psi_value_inv_svd + enddo + ! !!! + if (elec_beta_num /= 0) then + do ee = elec_alpha_num+1, elec_num + tmp = 0.d0 + do l = 1, n_svd_coefs + tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) * det_beta_pseudo_SVD(ee,l) + enddo + pseudo_non_local_SVD(ee) = tmp * psi_value_inv_svd + enddo + endif + ! !!! + endif + ! !!! + ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~ + +END_PROVIDER + + + + +! ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + + + + + + + + + + + + + BEGIN_PROVIDER [ double precision, det_alpha_pseudo_curr, (elec_alpha_num) ] implicit none BEGIN_DOC @@ -1869,3 +2233,9 @@ END_PROVIDER enddo END_PROVIDER + + + + + + diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 62c16af..a79c481 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -48,7 +48,10 @@ data = [ \ ("simulation_e_trial" , "double precision" , "" ), ("simulation_do_run" , "logical " , "" ), ("pseudo_do_pseudo" , "logical " , "" ), - +("spindeterminants_n_svd_coefs", "integer", ""), +("spindeterminants_psi_svd_alpha", "double precision", "(det_alpha_num,n_svd_coefs,n_states)"), +("spindeterminants_psi_svd_beta", "double precision", "(det_beta_num,n_svd_coefs,n_states)"), +("spindeterminants_psi_svd_coefs", "double precision", "(n_svd_coefs,n_states)") ] data_no_set = [\ diff --git a/src/psi.irp.f b/src/psi.irp.f index 3f5f5c7..21d85f8 100644 --- a/src/psi.irp.f +++ b/src/psi.irp.f @@ -4,10 +4,28 @@ ! Value of the wave function END_DOC - psi_value = psidet_value*jast_value + double precision :: norm_Err, psi_value2 + + if (utilise_svd) then + psi_value = psidet_value_svd*jast_value + else + psi_value = psidet_value*jast_value + endif + if (psi_value == 0.d0) then call abrt(irp_here,"Value of the wave function is 0.") endif + + !psi_value2 = psidet_value_svd * jast_value + !norm_Err = (psi_value2 - psi_value)**2 / psi_value**2 + !if (norm_Err > 1.d-6) then + ! print *, 'probleme dans PROVIDER [ psi_value ]: norm_Err = ', norm_Err + !print *, psi_value2 + !print *, psi_value + !print *, irp_here + !stop + !endif + END_PROVIDER BEGIN_PROVIDER [ double precision, psi_value_inv ] @@ -32,17 +50,37 @@ BEGIN_PROVIDER [ double precision, psi_lapl, (elec_num_8) ] BEGIN_DOC ! Laplacian of the wave function END_DOC - - integer :: i, j - !DIR$ VECTOR ALIGNED - !DIR$ LOOP COUNT (100) - do j=1,elec_num - psi_lapl(j) = jast_value*(psidet_grad_lapl(4,j) + psidet_value*jast_lapl_jast_inv(j) + 2.d0*(& - psidet_grad_lapl(1,j)*jast_grad_jast_inv_x(j) + & - psidet_grad_lapl(2,j)*jast_grad_jast_inv_y(j) + & - psidet_grad_lapl(3,j)*jast_grad_jast_inv_z(j) )) - enddo - + double precision :: psi_lapl2(elec_num) + double precision :: norm_Err + integer :: i, j + if (utilise_svd) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_lapl(j) = jast_value*(psidet_grad_lapl_svd(4,j) + psidet_value_svd*jast_lapl_jast_inv(j) + 2.d0*(& + psidet_grad_lapl_svd(1,j)*jast_grad_jast_inv_x(j) + & + psidet_grad_lapl_svd(2,j)*jast_grad_jast_inv_y(j) + & + psidet_grad_lapl_svd(3,j)*jast_grad_jast_inv_z(j) )) + enddo + else + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_lapl(j) = jast_value*(psidet_grad_lapl(4,j) + psidet_value*jast_lapl_jast_inv(j) + 2.d0*(& + psidet_grad_lapl(1,j)*jast_grad_jast_inv_x(j) + & + psidet_grad_lapl(2,j)*jast_grad_jast_inv_y(j) + & + psidet_grad_lapl(3,j)*jast_grad_jast_inv_z(j) )) + enddo + endif + + !norm_Err = sum( (psi_lapl2(1:elec_num) - psi_lapl(1:elec_num))**2 ) / sum( psi_lapl(1:elec_num)**2 ) + !if (norm_Err > 1.d-6) then + ! print *, 'probleme dans PROVIDER [ psi_lapl ]: norm_Err = ', norm_Err + !print *, psi_lapl2(1:elec_num) + !print *, psi_lapl(1:elec_num) + !print *, irp_here + !stop + !endif END_PROVIDER BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x, (elec_num_8) ] @@ -53,15 +91,36 @@ END_PROVIDER ! grad(psi)/psi END_DOC - integer :: j - !DIR$ VECTOR ALIGNED - !DIR$ LOOP COUNT (100) - do j=1,elec_num - psi_grad_psi_inv_x(j) = psidet_grad_lapl(1,j)*psidet_inv + jast_grad_jast_inv_x(j) - psi_grad_psi_inv_y(j) = psidet_grad_lapl(2,j)*psidet_inv + jast_grad_jast_inv_y(j) - psi_grad_psi_inv_z(j) = psidet_grad_lapl(3,j)*psidet_inv + jast_grad_jast_inv_z(j) - enddo - + double precision :: psi_grad_psi_inv_x2(elec_num), psi_grad_psi_inv_y2(elec_num), psi_grad_psi_inv_z2(elec_num) + double precision :: norm_Err + integer :: j + if (utilise_svd) then + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_grad_psi_inv_x(j) = psidet_grad_lapl_svd(1,j)*psidet_inv_svd + jast_grad_jast_inv_x(j) + psi_grad_psi_inv_y(j) = psidet_grad_lapl_svd(2,j)*psidet_inv_svd + jast_grad_jast_inv_y(j) + psi_grad_psi_inv_z(j) = psidet_grad_lapl_svd(3,j)*psidet_inv_svd + jast_grad_jast_inv_z(j) + enddo + else + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j=1,elec_num + psi_grad_psi_inv_x(j) = psidet_grad_lapl(1,j)*psidet_inv + jast_grad_jast_inv_x(j) + psi_grad_psi_inv_y(j) = psidet_grad_lapl(2,j)*psidet_inv + jast_grad_jast_inv_y(j) + psi_grad_psi_inv_z(j) = psidet_grad_lapl(3,j)*psidet_inv + jast_grad_jast_inv_z(j) + enddo + endif + !norm_Err = sum( (psi_grad_psi_inv_x2(1:elec_num) - psi_grad_psi_inv_x(1:elec_num))**2 ) / sum( psi_grad_psi_inv_x(1:elec_num)**2 ) + !norm_Err = norm_Err + sum( (psi_grad_psi_inv_y2(1:elec_num) - psi_grad_psi_inv_y(1:elec_num))**2 ) / sum( psi_grad_psi_inv_y(1:elec_num)**2 ) + !norm_Err = norm_Err + sum( (psi_grad_psi_inv_z2(1:elec_num) - psi_grad_psi_inv_z(1:elec_num))**2 ) / sum( psi_grad_psi_inv_z(1:elec_num)**2 ) + !if (norm_Err > 1.d-6) then + ! print *, 'probleme dans PROVIDER [ psi_grad_psi_inv_xyz ]: norm_Err = ', norm_Err + !print *, irp_here + !stop + !endif + + END_PROVIDER diff --git a/src/psi_SVD.irp.f b/src/psi_SVD.irp.f new file mode 100644 index 0000000..e218c65 --- /dev/null +++ b/src/psi_SVD.irp.f @@ -0,0 +1,85 @@ + BEGIN_PROVIDER [ double precision, psi_value_SVD ] + implicit none + BEGIN_DOC + ! Value of the wave function + END_DOC + psi_value_SVD = psidet_value_SVD * jast_value + if (psi_value_SVD == 0.d0) then + call abrt(irp_here,"Value of the wave function is 0.") + endif +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, psi_value_inv_SVD ] + implicit none + BEGIN_DOC + ! 1. / psi_value_SVD + END_DOC + psi_value_inv_SVD = 1.d0 / psi_value_SVD +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, psi_value_inv2_SVD ] + implicit none + BEGIN_DOC + ! 1./(psi_value_SVD)**2 + END_DOC + psi_value_inv2_SVD = psi_value_inv_SVD * psi_value_inv_SVD +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, psi_lapl_SVD, (elec_num_8) ] + implicit none + BEGIN_DOC + ! Laplacian of the wave function + END_DOC + integer :: i, j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j = 1, elec_num + psi_lapl_SVD(j) = jast_value * ( & + psidet_grad_lapl_SVD(4,j) + & + psidet_value_SVD * jast_lapl_jast_inv(j) + 2.d0 * ( & + psidet_grad_lapl_SVD(1,j) * jast_grad_jast_inv_x(j) + & + psidet_grad_lapl_SVD(2,j) * jast_grad_jast_inv_y(j) + & + psidet_grad_lapl_SVD(3,j) * jast_grad_jast_inv_z(j) ) ) + enddo +END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x_SVD, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_y_SVD, (elec_num_8) ] +&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_z_SVD, (elec_num_8) ] + implicit none + BEGIN_DOC + ! grad(psi_SVD) / psi_SVD + END_DOC + integer :: j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j = 1, elec_num + psi_grad_psi_inv_x_SVD(j) = psidet_grad_lapl_SVD(1,j) * psidet_inv_SVD + jast_grad_jast_inv_x(j) + psi_grad_psi_inv_y_SVD(j) = psidet_grad_lapl_SVD(2,j) * psidet_inv_SVD + jast_grad_jast_inv_y(j) + psi_grad_psi_inv_z_SVD(j) = psidet_grad_lapl_SVD(3,j) * psidet_inv_SVD + jast_grad_jast_inv_z(j) + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, psi_lapl_psi_inv_SVD, (elec_num_8) ] + implicit none + BEGIN_DOC + ! (Laplacian psi_SVD) / psi_SVD + END_DOC + integer :: i, j + !DIR$ VECTOR ALIGNED + !DIR$ LOOP COUNT (100) + do j = 1, elec_num + psi_lapl_psi_inv_SVD(j) = psi_lapl_SVD(j) * psi_value_inv_SVD + enddo +END_PROVIDER + diff --git a/src/wf.irp.f b/src/wf.irp.f index 4625193..fe00dd5 100644 --- a/src/wf.irp.f +++ b/src/wf.irp.f @@ -112,7 +112,11 @@ END_PROVIDER allocate (psi_det_tmp (N_int,max(det_alpha_num,det_beta_num))) - t = ci_threshold + if (utilise_svd) then + t = -1.d0 + else + t = ci_threshold + endif ! Compute the norm of the alpha and beta determinants d_alpha = 0.d0