LCPQ
/
qmcchem
mirror of https://gitlab.com/scemama/qmcchem.git
10
1
Fork 1
Code Releases Activity

Fast SVD

This commit is contained in:
Anthony Scemama 2021-06-21 14:30:29 +02:00
parent 3815c4ef35
commit 0898c98575
9 changed files with 227 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
ezfio_config/qmc.config
View File

@ -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

33
src/JASTROW/jastrow_simple.irp.f
View File

@ -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

3
src/MAIN/qmcchem_info.irp.f
View File

@ -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:'

4
src/PROPERTIES/properties_ci.irp.f
View File

@ -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

7
src/PROPERTIES/properties_energy.irp.f
View File

@ -276,10 +276,3 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
! E_loc_zv(:) = 0.d0
END_PROVIDER

143
src/det.irp.f
View File

@ -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), &

7
src/ezfio_interface.irp.f
View File

@ -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 = [\

45
src/svd.irp.f Normal file
View File

@ -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

6
src/wf.irp.f
View File

@ -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