From 0898c985758375a5ee58414866415775f1994f46 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 21 Jun 2021 14:30:29 +0200 Subject: [PATCH] Fast SVD --- ezfio_config/qmc.config | 15 ++- src/JASTROW/jastrow_simple.irp.f | 33 +++++- src/MAIN/qmcchem_info.irp.f | 3 + src/PROPERTIES/properties_ci.irp.f | 4 +- src/PROPERTIES/properties_energy.irp.f | 7 -- src/det.irp.f | 143 +++++++++++++++++++++---- src/ezfio_interface.irp.f | 7 +- src/svd.irp.f | 45 ++++++++ src/wf.irp.f | 6 +- 9 files changed, 227 insertions(+), 36 deletions(-) create mode 100644 src/svd.irp.f diff --git a/ezfio_config/qmc.config b/ezfio_config/qmc.config index 992033b..35ee84d 100644 --- a/ezfio_config/qmc.config +++ b/ezfio_config/qmc.config @@ -43,7 +43,20 @@ spindeterminants psi_coef_matrix_rows integer (spindeterminants_n_det) psi_coef_matrix_columns integer (spindeterminants_n_det) psi_coef_matrix_values double precision (spindeterminants_n_det,spindeterminants_n_states) - + n_svd_coefs_unique integer + n_svd_coefs integer + n_svd_selected integer + n_svd_toselect integer + psi_svd_alpha_unique double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs_unique,spindeterminants_n_states) + psi_svd_beta_unique double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs_unique,spindeterminants_n_states) + psi_svd_coefs_unique double precision (spindeterminants_n_svd_coefs_unique,spindeterminants_n_states) + psi_svd_alpha double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs,spindeterminants_n_states) + psi_svd_beta double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs,spindeterminants_n_states) + psi_svd_coefs double precision (spindeterminants_n_svd_coefs,spindeterminants_n_states) + psi_svd_alpha_numselected integer (spindeterminants_n_svd_selected,spindeterminants_n_states) + psi_svd_beta_numselected integer (spindeterminants_n_svd_selected,spindeterminants_n_states) + psi_svd_alpha_numtoselect integer (spindeterminants_n_svd_toselect,spindeterminants_n_states) + psi_svd_beta_numtoselect integer (spindeterminants_n_svd_toselect,spindeterminants_n_states) simulation do_run integer diff --git a/src/JASTROW/jastrow_simple.irp.f b/src/JASTROW/jastrow_simple.irp.f index bbd9c99..83085b1 100644 --- a/src/JASTROW/jastrow_simple.irp.f +++ b/src/JASTROW/jastrow_simple.irp.f @@ -20,10 +20,10 @@ implicit none enddo enddo - a = 0.5*jast_a_up_up + a = 0.5d0*jast_a_up_up b = jast_b_up_up - do j=1,elec_alpha_num + do j=1,elec_alpha_num-1 !DIR$ LOOP COUNT (50) do i=j+1,elec_alpha_num rij = elec_dist(i,j) @@ -33,7 +33,7 @@ implicit none enddo enddo - do j=elec_alpha_num+1,elec_num + do j=elec_alpha_num+1,elec_num-1 !DIR$ LOOP COUNT (50) do i=j+1,elec_num rij = elec_dist(i,j) @@ -43,7 +43,7 @@ implicit none enddo enddo - a = 0.5*jast_a_up_dn + a = 0.5d0*jast_a_up_dn b = jast_b_up_dn do j=1,elec_alpha_num @@ -190,5 +190,30 @@ BEGIN_PROVIDER [ double precision , jast_elec_Simple_lapl, (elec_num_8) ] enddo enddo +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, tmp1, tmp2 + + do j = 1, nucl_num + a = jast_pen(j) + tmp2 = 0.d0 + !DIR$ LOOP COUNT (100) + do i = 1, elec_num + rij = nucl_elec_dist(j,i) + tmp1 = (1.d0+a*rij)*(1.d0+a*rij)*(1.d0+a*rij) + tmp2 += rij*rij/tmp1 + end do + jast_elec_Simple_deriv_nucPar(j) = -2.d0 * a * tmp2 + end do END_PROVIDER diff --git a/src/MAIN/qmcchem_info.irp.f b/src/MAIN/qmcchem_info.irp.f index 6feb689..74d3109 100644 --- a/src/MAIN/qmcchem_info.irp.f +++ b/src/MAIN/qmcchem_info.irp.f @@ -13,6 +13,9 @@ program qmcchem_info endif print *, 'Number of determinants : ', det_num print *, 'Number of unique alpha/beta determinants : ', det_alpha_num, det_beta_num + if (use_svd) then + print *, 'SVD rank : ', n_svd_coefs + endif print *, 'Closed-shell MOs : ', mo_closed_num print *, 'Number of MOs in determinants : ', num_present_mos ! print *, 'Det alpha norm:' diff --git a/src/PROPERTIES/properties_ci.irp.f b/src/PROPERTIES/properties_ci.irp.f index a0d50fd..f2d5323 100644 --- a/src/PROPERTIES/properties_ci.irp.f +++ b/src/PROPERTIES/properties_ci.irp.f @@ -23,7 +23,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ] implicit none BEGIN_DOC - ! < Phi_0 | det(j) > + ! < Phi_0 | H | det(j) > ! ! Dimensions : det_num END_DOC @@ -54,7 +54,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, ci_overlap_matrix, (size_ci_overlap_matrix) ] implicit none BEGIN_DOC - ! < det(i) |H| det(j) > + ! < det(i) | det(j) > ! ! Dimensions : det_num*det_num END_DOC diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 84554a2..4deca87 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -276,10 +276,3 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ] ! E_loc_zv(:) = 0.d0 END_PROVIDER - - - - - - - diff --git a/src/det.irp.f b/src/det.irp.f index 6cfd595..ffd547b 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -1573,9 +1573,9 @@ END_PROVIDER endif det_alpha_value(det_i) = det_alpha_value_curr - det_alpha_grad_lapl(:,:,det_i) = det_alpha_grad_lapl_curr(:,:) + det_alpha_grad_lapl(1:4,1:elec_alpha_num,det_i) = det_alpha_grad_lapl_curr(1:4,1:elec_alpha_num) if (do_pseudo) then - det_alpha_pseudo(:,det_i) = det_alpha_pseudo_curr(:) + det_alpha_pseudo(1:elec_alpha_num,det_i) = det_alpha_pseudo_curr(1:elec_alpha_num) endif enddo @@ -1587,8 +1587,8 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [ double precision, det_beta_value, (det_beta_num_8) ] -&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_alpha_num+1:elec_num,det_beta_num) ] -&BEGIN_PROVIDER [ double precision, det_beta_pseudo, (elec_alpha_num+1:elec_num,det_beta_num_pseudo) ] +&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl, (4,elec_beta_num,det_beta_num) ] +&BEGIN_PROVIDER [ double precision, det_beta_pseudo, (det_beta_num_8,det_beta_num_pseudo) ] implicit none @@ -1622,9 +1622,11 @@ END_PROVIDER endif det_beta_value(det_j) = det_beta_value_curr - det_beta_grad_lapl(:,:,det_j) = det_beta_grad_lapl_curr(:,:) + det_beta_grad_lapl(1:4,1:elec_beta_num,det_j) = & + det_beta_grad_lapl_curr(1:4,elec_alpha_num+1:elec_num) if (do_pseudo) then - det_beta_pseudo(:,det_j) = det_beta_pseudo_curr(:) + det_beta_pseudo(1:elec_beta_num,det_j) = & + det_beta_pseudo_curr(elec_alpha_num+1:elec_num) endif enddo @@ -1729,10 +1731,93 @@ END_PROVIDER else - 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) - 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) + if (det_num == 1) then + + DaC(1) = det_alpha_value_curr + CDb(1) = det_beta_value_curr + + else if (use_svd) then + + DaC = 0.d0 + CDb = 0.d0 + double precision :: DaU(4), VtDb(4) + integer :: n_svd4 + + + n_svd4 = iand(n_svd_coefs,not(3)) + + do i=1,n_svd4,4 + + DaU = 0.d0 + do j=1,det_alpha_num + DaU(1) = DaU(1) + psi_svd_alpha(j,i+0) * det_alpha_value(j) + DaU(2) = DaU(2) + psi_svd_alpha(j,i+1) * det_alpha_value(j) + DaU(3) = DaU(3) + psi_svd_alpha(j,i+2) * det_alpha_value(j) + DaU(4) = DaU(4) + psi_svd_alpha(j,i+3) * det_alpha_value(j) + end do + DaU(1) = DaU(1)*psi_svd_coefs(i+0) + DaU(2) = DaU(2)*psi_svd_coefs(i+1) + DaU(3) = DaU(3)*psi_svd_coefs(i+2) + DaU(4) = DaU(4)*psi_svd_coefs(i+3) + + + VtDb = 0.d0 + do j=1,det_beta_num + VtDb(1) = VtDb(1) + psi_svd_beta(j,i+0) * det_beta_value(j) + VtDb(2) = VtDb(2) + psi_svd_beta(j,i+1) * det_beta_value(j) + VtDb(3) = VtDb(3) + psi_svd_beta(j,i+2) * det_beta_value(j) + VtDb(4) = VtDb(4) + psi_svd_beta(j,i+3) * det_beta_value(j) + + DaC(j) = DaC(j) + DaU(1) * psi_svd_beta(j,i+0) + & + DaU(2) * psi_svd_beta(j,i+1) + & + DaU(3) * psi_svd_beta(j,i+2) + & + DaU(4) * psi_svd_beta(j,i+3) + end do + VtDb(1) = VtDb(1)*psi_svd_coefs(i+0) + VtDb(2) = VtDb(2)*psi_svd_coefs(i+1) + VtDb(3) = VtDb(3)*psi_svd_coefs(i+2) + VtDb(4) = VtDb(4)*psi_svd_coefs(i+3) + + do j=1,det_alpha_num + CDb(j) = CDb(j) + VtDb(1) * psi_svd_alpha(j,i+0) + & + VtDb(2) * psi_svd_alpha(j,i+1) + & + VtDb(3) * psi_svd_alpha(j,i+2) + & + VtDb(4) * psi_svd_alpha(j,i+3) + end do + + end do + + + do i=n_svd4+1,n_svd_coefs + + DaU(1) = 0.d0 + do j=1,det_alpha_num + DaU(1) = DaU(1) + psi_svd_alpha(j,i) * det_alpha_value(j) + end do + DaU(1) = DaU(1)*psi_svd_coefs(i) + + VtDb(1) = 0.d0 + do j=1,det_beta_num + VtDb(1) = VtDb(1) + psi_svd_beta(j,i) * det_beta_value(j) + DaC(j) = DaC(j) + DaU(1) * psi_svd_beta(j,i) + end do + VtDb(1) = VtDb(1)*psi_svd_coefs(i) + + do j=1,det_alpha_num + CDb(j) = CDb(j) + VtDb(1) * psi_svd_alpha(j,i) + end do + + end do + + else + + 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) + + 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) + + endif endif @@ -1750,20 +1835,38 @@ END_PROVIDER endif psidet_inv = 1.d0/psidet_value + ! Gradients ! --------- - - 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), & - CDb, 1, 0.d0, psidet_grad_lapl, 1) - if (elec_beta_num /= 0) then - call dgemv('N',elec_beta_num*4,det_beta_num,1.d0, & - det_beta_grad_lapl(1,elec_alpha_num+1,1), & - size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2), & - DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1) + if(det_num .eq. 1) then + do i = 1, elec_alpha_num + psidet_grad_lapl(1,i) = det_alpha_grad_lapl_curr(1,i) * det_beta_value_curr + psidet_grad_lapl(2,i) = det_alpha_grad_lapl_curr(2,i) * det_beta_value_curr + psidet_grad_lapl(3,i) = det_alpha_grad_lapl_curr(3,i) * det_beta_value_curr + psidet_grad_lapl(4,i) = det_alpha_grad_lapl_curr(4,i) * det_beta_value_curr + enddo + do i = elec_alpha_num+1, elec_num + psidet_grad_lapl(1,i) = det_beta_grad_lapl_curr(1,i) * det_alpha_value_curr + psidet_grad_lapl(2,i) = det_beta_grad_lapl_curr(2,i) * det_alpha_value_curr + psidet_grad_lapl(3,i) = det_beta_grad_lapl_curr(3,i) * det_alpha_value_curr + psidet_grad_lapl(4,i) = det_beta_grad_lapl_curr(4,i) * det_alpha_value_curr + enddo + else + ! 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), & + CDb, 1, 0.d0, psidet_grad_lapl, 1) + if (elec_beta_num /= 0) then + call dgemv('N',elec_beta_num*4,det_beta_num,1.d0, & + det_beta_grad_lapl, & + size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2), & + DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1) + endif endif + if (do_pseudo) then call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv, & det_alpha_pseudo, size(det_alpha_pseudo,1), & diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 62c16af..09f032f 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -48,7 +48,12 @@ data = [ \ ("simulation_e_trial" , "double precision" , "" ), ("simulation_do_run" , "logical " , "" ), ("pseudo_do_pseudo" , "logical " , "" ), - +("spindeterminants_n_svd_coefs" , "integer", ""), +("spindeterminants_n_svd_selected" , "integer", ""), +("spindeterminants_n_svd_toselect" , "integer", ""), +("spindeterminants_psi_svd_alpha", "double precision", "(det_alpha_num,n_svd_coefs_full,n_states)"), +("spindeterminants_psi_svd_beta" , "double precision", "(det_beta_num,n_svd_coefs_full,n_states)"), +("spindeterminants_psi_svd_coefs", "double precision", "(n_svd_coefs_full,n_states)"), ] data_no_set = [\ diff --git a/src/svd.irp.f b/src/svd.irp.f new file mode 100644 index 0000000..382f2f9 --- /dev/null +++ b/src/svd.irp.f @@ -0,0 +1,45 @@ + BEGIN_PROVIDER [ logical, use_svd ] +&BEGIN_PROVIDER [ integer, n_svd_coefs_full ] + implicit none + BEGIN_DOC + ! If true, use SVD wave function + END_DOC + n_svd_coefs_full = -1 + call get_spindeterminants_n_svd_coefs(n_svd_coefs_full) + use_svd = n_svd_coefs_full > 0 + if (.not.use_SVD) then + n_svd_coefs_full = 1 + endif +END_PROVIDER + +BEGIN_PROVIDER [ integer, n_svd_coefs ] + implicit none + BEGIN_DOC + ! If true, use SVD wave function + END_DOC + integer :: i + double precision :: keep + keep = 1.d0 + n_svd_coefs = 0 + do i=1,n_svd_coefs_full + keep = keep - psi_svd_coefs(i)*psi_svd_coefs(i) + if (dsqrt(keep) < ci_threshold) then + exit + endif + n_svd_coefs = n_svd_coefs+1 + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, psi_svd_coefs, ( n_svd_coefs_full) ] +&BEGIN_PROVIDER [ double precision, psi_svd_alpha, (det_alpha_num, n_svd_coefs_full) ] +&BEGIN_PROVIDER [ double precision, psi_svd_beta , (det_beta_num , n_svd_coefs_full) ] + implicit none + BEGIN_DOC + ! !!! + ! truncated SVD + END_DOC + call get_spindeterminants_psi_svd_coefs(psi_svd_coefs) + call get_spindeterminants_psi_svd_alpha(psi_svd_alpha) + call get_spindeterminants_psi_svd_beta(psi_svd_beta) + END_PROVIDER + diff --git a/src/wf.irp.f b/src/wf.irp.f index 4625193..c922b62 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 (use_svd) then + t = -1.d0 + else + t = ci_threshold + endif ! Compute the norm of the alpha and beta determinants d_alpha = 0.d0