diff --git a/fnmf/EZFIO.cfg b/fnmf/EZFIO.cfg index 86b2eb0..033d0f1 100644 --- a/fnmf/EZFIO.cfg +++ b/fnmf/EZFIO.cfg @@ -3,15 +3,15 @@ type: double precision doc: DMC energy interface: ezfio, provider -[dmc_h] +[h_dmc] type: double precision doc: Dressing matrix obtained from DMC size: (determinants.n_det) interface: ezfio, provider -[dmc_s] +[s_dmc] type: double precision doc: Dressing matrix obtained from H_TC size: (determinants.n_det) -interface: ezfio +interface: ezfio, provider diff --git a/fnmf/dmc_data.irp.f b/fnmf/dmc_data.irp.f index a275aa2..84257ce 100644 --- a/fnmf/dmc_data.irp.f +++ b/fnmf/dmc_data.irp.f @@ -5,17 +5,17 @@ ! Data sampled with QMC=Chem END_DOC -! h_dmc_row(:) = h_dmc(:) -! s_dmc_row(:) = s_dmc(:) -! call dset_order(h_dmc_row,psi_bilinear_matrix_order_reverse,N_det) -! call dset_order(s_dmc_row,psi_bilinear_matrix_order_reverse,N_det) + h_dmc_row(:) = h_dmc(:) + s_dmc_row(:) = s_dmc(:) + call dset_order(h_dmc_row,psi_bilinear_matrix_order_reverse,N_det) + call dset_order(s_dmc_row,psi_bilinear_matrix_order_reverse,N_det) - integer :: i - do i=1,N_det - s_dmc_row(i) = psi_coef(i,1) - call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, & - N_det, 1, h_dmc_row(i) ) - enddo +! integer :: i +! do i=1,N_det +! s_dmc_row(i) = psi_coef(i,1) +! call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, & +! N_det, 1, h_dmc_row(i) ) +! enddo END_PROVIDER @@ -44,7 +44,7 @@ BEGIN_PROVIDER [ double precision, H_dmc_mat, (mat_size, mat_size) ] N_det, 1, H_dmc_mat(i,N_det+1) ) H_dmc_mat(N_det+1,i) = h_dmc_row(i) enddo - H_dmc_mat(mat_size,mat_size) = E_dmc + H_dmc_mat(mat_size,mat_size) = E_dmc - nuclear_repulsion END_PROVIDER diff --git a/fnmf/fnmf.irp.f b/fnmf/fnmf.irp.f index bb73d9e..fff7a6a 100644 --- a/fnmf/fnmf.irp.f +++ b/fnmf/fnmf.irp.f @@ -13,15 +13,64 @@ subroutine run implicit none integer :: i, n_real double precision, allocatable :: beta(:), vr(:,:), vl(:,:) + double precision, allocatable :: htmp(:,:), stmp(:,:) - allocate(beta(mat_size), vr(mat_size,mat_size), vl(mat_size, mat_size)) - call lapack_g_non_sym_real(mat_size, H_dmc_mat, size(H_dmc_mat,1), & - S_dmc_mat, size(S_dmc_mat,1), beta, n_real, vl, size(vl,1), vr, size(vr,1)) + allocate(Htmp(N_det,N_det), Stmp(N_det,N_det)) - print *, 'EV VR VL' - print *, '---------------------------' - do i=1,mat_size - print '(3(F16.12,X))', beta(i), vr(i,1), vl(i,1) + htmp(:, :) = H_dmc_mat(1:N_det, 1:N_det) + stmp(:, :) = S_dmc_mat(1:N_det, 1:N_det) + + htmp(1, 1) = H_dmc_mat(N_det+1, N_det+1) + stmp(1, 1) = S_dmc_mat(N_det+1, N_det+1) + + htmp(2:N_det, 1) = H_dmc_mat(2:N_det, N_det+1) + stmp(2:N_det, 1) = S_dmc_mat(2:N_det, N_det+1) + + htmp(1, 2:N_det) = H_dmc_mat(N_det+1, 2:N_det) + stmp(1, 2:N_det) = S_dmc_mat(N_det+1, 2:N_det) + + + allocate(beta(N_det), vr(N_det,N_det), vl(N_det, N_det)) + + call lapack_g_non_sym_real(N_det, Htmp, size(Htmp,1), & + Stmp, size(Stmp,1), beta, n_real, vl, size(vl,1), vr, size(vr,1)) + + + vl(:,1) = psi_coef(:,1) + + psi_coef(:,1) *= vr(1,1) + do i=2,N_det + psi_coef(i,1) += vr(i,1) + end do + double precision, external :: dnrm2, ddot + double precision :: norm + + norm = dnrm2(N_det, psi_coef, 1) + psi_coef(:,1) *= 1.d0/norm + SOFT_TOUCH psi_coef + + print *, '-------------------------------------------------' + print *, ' psi_old psi_new' + print *, '-------------------------------------------------' + do i=1,N_det + print '(2(E16.8,X))', vl(i,1), psi_coef(i,1) enddo + print *, '-------------------------------------------------' + + print *, '' + print *, 'Overlap : ', ddot(N_det, psi_coef, 1, vl, 1) + print *, 'DMC energy : ', htmp(1,1) + nuclear_repulsion + print *, 'Updated energy : ', beta(1) + nuclear_repulsion + + print *, 'H' + do i=1,N_det + print *, htmp(i, 1:N_det) + enddo + print *, 'S' + do i=1,N_det + print *, stmp(i, 1:N_det) + enddo + call save_wavefunction + end diff --git a/fnmf/non_hermit.irp.f b/fnmf/non_hermit.irp.f index 6479c05..058754c 100644 --- a/fnmf/non_hermit.irp.f +++ b/fnmf/non_hermit.irp.f @@ -7,21 +7,24 @@ subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, & integer, intent(out) :: n_real integer :: lwork, info, i,j - double precision, allocatable :: work(:) + double precision, allocatable :: work(:), Htmp(:,:), Stmp(:,:) double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:) - integer, allocatable :: iorder(:) + integer, allocatable :: iorder(:), list_good(:) lwork = -1 - allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n)) + allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), & + Htmp(n,n), Stmp(n,n), list_good(n)) - call dggev('V', 'V', n, H, size(H,1), S, size(S,1), alphar, alphai, beta, & + Htmp(1:n,1:n) = H(1:n,1:n) + Stmp(1:n,1:n) = S(1:n,1:n) + call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, & vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info) lwork = int(work(1)) deallocate(work) allocate(work(lwork)) - call dggev('V', 'V', n, H, size(H,1), S, size(S,1), alphar, alphai, beta, & + call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, & vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info) deallocate(work) @@ -32,15 +35,31 @@ subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, & allocate(iorder(n)) n_real = 0 do i=1,n - iorder(i) = i if (dabs(alphai(i)) < 1.d-10) then n_real += 1 + list_good(n_real) = i + else alphar(i) = dble(huge(1.0)) endif enddo - beta(:) = alphar(:)/beta(:) - call dsort(beta, iorder, n) + do i=1,n + if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then + beta(i) = 0.d0 + else + beta(i) = alphar(i)/beta(i) + endif + enddo + + do i=1,n_real + iorder(i) = i + do j=1,n + vrtmp(j,i) = vrtmp(j,list_good(i)) + vltmp(j,i) = vltmp(j,list_good(i)) + end do + end do + + call dsort(beta, iorder, n_real) do i=1,n_real do j=1,n @@ -49,7 +68,7 @@ subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, & end do end do - deallocate(vrtmp, vltmp, iorder) + deallocate(vrtmp, vltmp, iorder, Htmp, Stmp) end