From 67e2818ee24ec40fc6504e01022b1aa53d16b846 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 23 Jul 2019 15:12:25 +0200 Subject: [PATCH] Added QMCChem module --- devel/qmcchem/EZFIO.cfg | 6 + devel/qmcchem/NEED | 4 + devel/qmcchem/README.rst | 12 ++ devel/qmcchem/expand_spindets_qmcchem.irp.f | 8 ++ devel/qmcchem/pot_ao_pseudo_ints.irp.f | 131 ++++++++++++++++++++ devel/qmcchem/pseudo.irp.f | 8 ++ devel/qmcchem/qmc_e_curve.irp.f | 113 +++++++++++++++++ devel/qmcchem/save_for_qmcchem.irp.f | 45 +++++++ devel/qmcchem/truncate_wf_qmcchem.irp.f | 116 +++++++++++++++++ 9 files changed, 443 insertions(+) create mode 100644 devel/qmcchem/EZFIO.cfg create mode 100644 devel/qmcchem/NEED create mode 100644 devel/qmcchem/README.rst create mode 100644 devel/qmcchem/expand_spindets_qmcchem.irp.f create mode 100644 devel/qmcchem/pot_ao_pseudo_ints.irp.f create mode 100644 devel/qmcchem/pseudo.irp.f create mode 100644 devel/qmcchem/qmc_e_curve.irp.f create mode 100644 devel/qmcchem/save_for_qmcchem.irp.f create mode 100644 devel/qmcchem/truncate_wf_qmcchem.irp.f diff --git a/devel/qmcchem/EZFIO.cfg b/devel/qmcchem/EZFIO.cfg new file mode 100644 index 0000000..c9d238f --- /dev/null +++ b/devel/qmcchem/EZFIO.cfg @@ -0,0 +1,6 @@ +[ci_threshold] +type: Threshold +doc: Threshold on the CI coefficients as computed in QMCChem +interface: ezfio,provider,ocaml +default: 1.e-8 + diff --git a/devel/qmcchem/NEED b/devel/qmcchem/NEED new file mode 100644 index 0000000..afb53ca --- /dev/null +++ b/devel/qmcchem/NEED @@ -0,0 +1,4 @@ +determinants +davidson_undressed +fci +zmq diff --git a/devel/qmcchem/README.rst b/devel/qmcchem/README.rst new file mode 100644 index 0000000..0c72b8e --- /dev/null +++ b/devel/qmcchem/README.rst @@ -0,0 +1,12 @@ +============== +QmcChem Module +============== + +For multi-state calculations, to extract state 2 use: + +`` +QP_STATE=2 qp_run save_for_qmcchem x.ezfio +`` +(state 1 is the ground state). + + diff --git a/devel/qmcchem/expand_spindets_qmcchem.irp.f b/devel/qmcchem/expand_spindets_qmcchem.irp.f new file mode 100644 index 0000000..2e1894b --- /dev/null +++ b/devel/qmcchem/expand_spindets_qmcchem.irp.f @@ -0,0 +1,8 @@ +program densify + implicit none + read_wf = .True. + touch read_wf + call generate_all_alpha_beta_det_products() + call diagonalize_ci + call save_wavefunction +end diff --git a/devel/qmcchem/pot_ao_pseudo_ints.irp.f b/devel/qmcchem/pot_ao_pseudo_ints.irp.f new file mode 100644 index 0000000..97ac748 --- /dev/null +++ b/devel/qmcchem/pot_ao_pseudo_ints.irp.f @@ -0,0 +1,131 @@ + +BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ] + implicit none + BEGIN_DOC +! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C +! +! + END_DOC + ! l,m : Y(l,m) parameters + ! c(3) : pseudopotential center + ! a(3) : Atomic Orbital center + ! n_a(3) : Powers of x,y,z in the Atomic Orbital + ! g_a : Atomic Orbital exponent + ! r : Distance between the Atomic Orbital center and the considered point + double precision, external :: ylm_orb + integer :: n_a(3) + double precision :: a(3), c(3), g_a + integer :: i,j,k,l,m,n,p + double precision :: dr, Ulc + double precision :: y + double precision, allocatable :: r(:) + + allocate (r(pseudo_grid_size)) + dr = pseudo_grid_rmax/dble(pseudo_grid_size) + r(1) = 0.d0 + do j=2,pseudo_grid_size + r(j) = r(j-1) + dr + enddo + + ao_pseudo_grid = 0.d0 + do j=1,pseudo_grid_size + do k=1,nucl_num + c(1:3) = nucl_coord(k,1:3) + do l=0,pseudo_lmax + do i=1,ao_num + a(1:3) = nucl_coord(ao_nucl(i),1:3) + n_a(1:3) = ao_power(i,1:3) + do p=1,ao_prim_num(i) + g_a = ao_expo_ordered_transp(p,i) + do m=-l,l + y = ylm_orb(l,m,c,a,n_a,g_a,r(j)) + ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + & + ao_coef_normalized_ordered_transp(p,i)*y + enddo + enddo + enddo + enddo + enddo + enddo + deallocate(r) + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ] + implicit none + BEGIN_DOC +! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C +! +! + END_DOC + ! l,m : Y(l,m) parameters + ! c(3) : pseudopotential center + ! a(3) : Atomic Orbital center + ! n_a(3) : Powers of x,y,z in the Atomic Orbital + ! g_a : Atomic Orbital exponent + ! r : Distance between the Atomic Orbital center and the considered point + double precision, external :: ylm_orb + integer :: n_a(3) + double precision :: a(3), c(3), g_a + integer :: i,j,k,l,m,n,p + double precision :: dr, Ulc + double precision :: y + double precision, allocatable :: r(:) + + allocate (r(pseudo_grid_size)) + + dr = pseudo_grid_rmax/dble(pseudo_grid_size) + r(1) = 0.d0 + do j=2,pseudo_grid_size + r(j) = r(j-1) + dr + enddo + + mo_pseudo_grid = 0.d0 + do n=1,pseudo_grid_size + do k=1,nucl_num + do l=0,pseudo_lmax + do m=-l,l + do i=1,ao_num + do j=1,mo_num + if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then + cycle + endif + if (dabs(mo_coef(i,j)) < 1.e-8) then + cycle + endif + mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + & + ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j) + enddo + enddo + enddo + enddo + enddo + enddo + deallocate(r) + +END_PROVIDER + +double precision function test_pseudo_grid_ao(i,j) + implicit none + integer, intent(in) :: i,j + integer :: k,l,m,n + double precision :: r, dr,u + dr = pseudo_grid_rmax/dble(pseudo_grid_size) + + test_pseudo_grid_ao = 0.d0 + r = 0.d0 + do k=1,pseudo_grid_size + do n=1,nucl_num + do l = 0,pseudo_lmax + u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr + do m=-l,l + test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u + enddo + enddo + enddo + r = r+dr + enddo +end diff --git a/devel/qmcchem/pseudo.irp.f b/devel/qmcchem/pseudo.irp.f new file mode 100644 index 0000000..1398100 --- /dev/null +++ b/devel/qmcchem/pseudo.irp.f @@ -0,0 +1,8 @@ +subroutine write_pseudopotential + implicit none + BEGIN_DOC +! Write the pseudo_potential into the EZFIO file + END_DOC + call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid) +end + diff --git a/devel/qmcchem/qmc_e_curve.irp.f b/devel/qmcchem/qmc_e_curve.irp.f new file mode 100644 index 0000000..c0a3b7d --- /dev/null +++ b/devel/qmcchem/qmc_e_curve.irp.f @@ -0,0 +1,113 @@ +program e_curve + use bitmasks + implicit none + integer :: i,j,k, kk, nab, m, l + double precision :: norm, E, hij, num, ci, cj + integer, allocatable :: iorder(:) + double precision , allocatable :: norm_sort(:) + double precision :: e_0(N_states) + PROVIDE mo_two_e_integrals_in_map + + nab = n_det_alpha_unique+n_det_beta_unique + allocate ( norm_sort(0:nab), iorder(0:nab) ) + + double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:) + double precision, allocatable :: u_0(:,:), v_0(:,:) + allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det)) + allocate(u_0(N_states,N_det),v_0(N_states,N_det)) + + + + norm_sort(0) = 0.d0 + iorder(0) = 0 + do i=1,n_det_alpha_unique + norm_sort(i) = det_alpha_norm(i) + iorder(i) = i + enddo + + do i=1,n_det_beta_unique + norm_sort(i+n_det_alpha_unique) = det_beta_norm(i) + iorder(i+n_det_alpha_unique) = -i + enddo + + call dsort(norm_sort(1),iorder(1),nab) + + if (.not.read_wf) then + stop 'Please set read_wf to true' + endif + + PROVIDE psi_bilinear_matrix_values nuclear_repulsion + print *, '' + print *, '==============================' + print *, 'Energies at different cut-offs' + print *, '==============================' + print *, '' + print *, '==========================================================' + print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)', 'Thresh.', 'Ndet', 'Cost', 'Norm', 'E' + print *, '==========================================================' + double precision :: thresh + integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:) + thresh = 1.d-10 + do j=0,nab + i = iorder(j) + if (i<0) then + do k=1,n_det + if (psi_bilinear_matrix_columns(k) == -i) then + psi_bilinear_matrix_values(k,1) = 0.d0 + endif + enddo + else + do k=1,n_det + if (psi_bilinear_matrix_rows(k) == i) then + psi_bilinear_matrix_values(k,1) = 0.d0 + endif + enddo + endif + if (thresh > norm_sort(j)) then + cycle + endif + + u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states) + v_t = 0.d0 + s_t = 0.d0 + call dtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_states) + call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1) + call dtranspose( & + v_t, & + size(v_t, 1), & + v_0, & + size(v_0, 1), & + N_states, N_det) + + double precision, external :: u_dot_u, u_dot_v + do i=1,N_states + e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det) + enddo + + m = 0 + do k=1,n_det + if (psi_bilinear_matrix_values(k,1) /= 0.d0) then + m = m+1 + endif + enddo + + if (m == 0) then + exit + endif + E = E_0(1) + nuclear_repulsion + norm = u_dot_u(u_0(1,1),N_det) + print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', thresh, m, & + dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / & + dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E + thresh = thresh * dsqrt(10.d0) + enddo + print *, '==========================================================' + + deallocate (iorder, norm_sort) +end + diff --git a/devel/qmcchem/save_for_qmcchem.irp.f b/devel/qmcchem/save_for_qmcchem.irp.f new file mode 100644 index 0000000..c7bcf46 --- /dev/null +++ b/devel/qmcchem/save_for_qmcchem.irp.f @@ -0,0 +1,45 @@ +program save_for_qmc + + integer :: iunit + logical :: exists + double precision :: e_ref + + ! Determinants + read_wf = .True. + TOUCH read_wf + print *, "N_det = ", N_det + call write_spindeterminants + + ! Reference Energy + if (do_pseudo) then + call write_pseudopotential + endif + call system( & + 'mkdir -p '//trim(ezfio_filename)//'/simulation ;' // & + 'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/simulation/.version ; ' // & + 'mkdir -p '//trim(ezfio_filename)//'/properties ;' // & + 'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/properties/.version ; ' // & + 'echo T > '//trim(ezfio_filename)//'/properties/e_loc' & + ) + iunit = 13 + open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write') + call ezfio_has_fci_energy_pt2(exists) + if (exists) then + call ezfio_get_fci_energy_pt2(e_ref) + else + call ezfio_has_fci_energy(exists) + if (exists) then + call ezfio_get_fci_energy(e_ref) + else + call ezfio_has_hartree_fock_energy(exists) + if (exists) then + call ezfio_get_hartree_fock_energy(e_ref) + else + e_ref = 0.d0 + endif + endif + endif + write(iunit,*) e_ref + close(iunit) + +end diff --git a/devel/qmcchem/truncate_wf_qmcchem.irp.f b/devel/qmcchem/truncate_wf_qmcchem.irp.f new file mode 100644 index 0000000..bc33499 --- /dev/null +++ b/devel/qmcchem/truncate_wf_qmcchem.irp.f @@ -0,0 +1,116 @@ +program truncate + read_wf = .True. + SOFT_TOUCH read_wf + call run +end + +subroutine run + use bitmasks + implicit none + integer :: i,j,k, kk, nab, m, l + double precision :: norm, E, hij, num, ci, cj + integer, allocatable :: iorder(:) + double precision , allocatable :: norm_sort(:) + double precision :: e_0(N_states) + + PROVIDE mo_two_e_integrals_in_map H_apply_buffer_allocated + + nab = n_det_alpha_unique+n_det_beta_unique + allocate ( norm_sort(0:nab), iorder(0:nab) ) + + integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:) + double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:) + double precision, allocatable :: u_0(:,:), v_0(:,:) + allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det)) + allocate(u_0(N_det,N_states),v_0(N_det,N_states)) + + norm_sort(0) = 0.d0 + iorder(0) = 0 + do i=1,n_det_alpha_unique + norm_sort(i) = det_alpha_norm(i) + iorder(i) = i + enddo + + do i=1,n_det_beta_unique + norm_sort(i+n_det_alpha_unique) = det_beta_norm(i) + iorder(i+n_det_alpha_unique) = -i + enddo + + call dsort(norm_sort(1),iorder(1),nab) + + + PROVIDE psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns + PROVIDE nuclear_repulsion + print *, '' + do j=0,nab + i = iorder(j) + if (i<0) then + !$OMP PARALLEL DO PRIVATE(k) + do k=1,n_det + if (psi_bilinear_matrix_columns(k) == -i) then + do l=1,N_states + psi_bilinear_matrix_values(k,l) = 0.d0 + enddo + endif + enddo + !$OMP END PARALLEL DO + else + !$OMP PARALLEL DO PRIVATE(k) + do k=1,n_det + if (psi_bilinear_matrix_rows(k) == i) then + do l=1,N_states + psi_bilinear_matrix_values(k,l) = 0.d0 + enddo + endif + enddo + !$OMP END PARALLEL DO + endif + if (ci_threshold <= norm_sort(j)) then + exit + endif + enddo + + m = 0 + do k=1,n_det + if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then + m = m+1 + endif + enddo + + do k=1,N_states + E = E_0(k) + nuclear_repulsion + enddo + print *, 'Number of determinants:', m + call wf_of_psi_bilinear_matrix(.True.) + call save_wavefunction + + u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states) + v_0(1:N_det,1:N_states) = 0.d0 + u_t(1:N_states,1:N_det) = 0.d0 + v_t(1:N_states,1:N_det) = 0.d0 + s_t(1:N_states,1:N_det) = 0.d0 + call dtranspose( & + u_0, & + size(u_0, 1), & + u_t, & + size(u_t, 1), & + N_det, N_states) + print *, 'Computing H|Psi> ...' + call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1) + print *, 'Done' + call dtranspose( & + v_t, & + size(v_t, 1), & + v_0, & + size(v_0, 1), & + N_states, N_det) + + double precision, external :: u_dot_u, u_dot_v + do i=1,N_states + e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det) + print *, 'E(',i,') = ', e_0(i) + nuclear_repulsion + enddo + + deallocate (iorder, norm_sort) +end +