From 5c66e4b99f5d6ba7b0e46fdc506949d24e305c5f Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 19 Feb 2020 17:59:27 -0600 Subject: [PATCH] complex determinants --- src/determinants/occ_pattern.irp.f | 4 + src/determinants/prune_wf.irp.f | 4 + src/determinants/psi_cas.irp.f | 16 ++- src/determinants/psi_cas_complex.irp.f | 145 +++++++++++++++++++++++++ src/utils_complex/qp2-pbc-diff.txt | 4 +- 5 files changed, 171 insertions(+), 2 deletions(-) create mode 100644 src/determinants/psi_cas_complex.irp.f diff --git a/src/determinants/occ_pattern.irp.f b/src/determinants/occ_pattern.irp.f index 781c7d75..9e0ccd33 100644 --- a/src/determinants/occ_pattern.irp.f +++ b/src/determinants/occ_pattern.irp.f @@ -529,7 +529,11 @@ subroutine make_s2_eigenfunction if (update) then call copy_H_apply_buffer_to_wf + if (is_complex) then + TOUCH N_det psi_coef_complex psi_det psi_occ_pattern N_occ_pattern + else TOUCH N_det psi_coef psi_det psi_occ_pattern N_occ_pattern + endif endif call write_time(6) diff --git a/src/determinants/prune_wf.irp.f b/src/determinants/prune_wf.irp.f index c3cd8d12..136d4ec1 100644 --- a/src/determinants/prune_wf.irp.f +++ b/src/determinants/prune_wf.irp.f @@ -25,7 +25,11 @@ BEGIN_PROVIDER [ logical, pruned, (N_det) ] else ndet_new = max(1,int( dble(N_det) * (1.d0 - pruning) + 0.5d0 )) + if (is_complex) then + thr = psi_average_norm_contrib_sorted_complex(ndet_new) + else thr = psi_average_norm_contrib_sorted(ndet_new) + endif do i=1, N_det pruned(i) = psi_average_norm_contrib(i) < thr enddo diff --git a/src/determinants/psi_cas.irp.f b/src/determinants/psi_cas.irp.f index 19a1c260..d262622c 100644 --- a/src/determinants/psi_cas.irp.f +++ b/src/determinants/psi_cas.irp.f @@ -150,7 +150,20 @@ END_PROVIDER double precision :: hij,norm,u_dot_v psi_cas_energy = 0.d0 - + if (is_complex) then + complex*16 :: hij_c + do k = 1, N_states + norm = 0.d0 + do i = 1, N_det_cas_complex + norm += cdabs(psi_cas_coef_complex(i,k) * psi_cas_coef_complex(i,k)) + do j = 1, N_det_cas_complex + !TODO: accum imag parts to ensure that sum is zero? + psi_cas_energy(k) += dble(dconjg(psi_cas_coef_complex(i,k)) * psi_cas_coef_complex(j,k) * H_matrix_cas_complex(i,j)) + enddo + enddo + psi_cas_energy(k) = psi_cas_energy(k) /norm + enddo + else do k = 1, N_states norm = 0.d0 do i = 1, N_det_cas @@ -161,6 +174,7 @@ END_PROVIDER enddo psi_cas_energy(k) = psi_cas_energy(k) /norm enddo + endif END_PROVIDER diff --git a/src/determinants/psi_cas_complex.irp.f b/src/determinants/psi_cas_complex.irp.f new file mode 100644 index 00000000..9e8ded87 --- /dev/null +++ b/src/determinants/psi_cas_complex.irp.f @@ -0,0 +1,145 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas_complex, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ complex*16, psi_cas_coef_complex, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas_complex, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_cas_complex ] + implicit none + BEGIN_DOC + ! |CAS| wave function, defined from the application of the |CAS| bitmask on the + ! determinants. idx_cas gives the indice of the |CAS| determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + n_det_cas_complex = 0 + do i=1,N_det + do l = 1, N_states + psi_cas_coef_complex(i,l) = (0.d0,0.d0) + enddo + good = .True. + do k=1,N_int + good = good .and. ( & + iand(not(act_bitmask(k,1)), psi_det(k,1,i)) == & + iand(not(act_bitmask(k,1)), hf_bitmask(k,1)) ) .and. ( & + iand(not(act_bitmask(k,2)), psi_det(k,2,i)) == & + iand(not(act_bitmask(k,2)), hf_bitmask(k,2)) ) + enddo + if (good) then + exit + endif + if (good) then + n_det_cas_complex = n_det_cas_complex+1 + do k=1,N_int + psi_cas_complex(k,1,n_det_cas_complex) = psi_det(k,1,i) + psi_cas_complex(k,2,n_det_cas_complex) = psi_det(k,2,i) + enddo + idx_cas(n_det_cas_complex) = i + do k=1,N_states + psi_cas_coef_complex(n_det_cas_complex,k) = psi_coef_complex(i,k) + enddo + endif + enddo + call write_int(6,n_det_cas_complex, 'Number of determinants in the CAS') + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas_sorted_bit_complex, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ complex*16, psi_cas_coef_sorted_bit_complex, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! |CAS| determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key_complex(n_det_cas_complex, psi_cas_complex, psi_cas_coef_complex, size(psi_cas_coef_complex,1), & + psi_cas_sorted_bit_complex, psi_cas_coef_sorted_bit_complex, N_states) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas_complex, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ complex*16, psi_non_cas_coef,_complex (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_non_cas_complex, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_non_cas_complex ] + implicit none + BEGIN_DOC + ! Set of determinants which are not part of the |CAS|, defined from the application + ! of the |CAS| bitmask on the determinants. + ! idx_non_cas gives the indice of the determinant in psi_det. + END_DOC + integer :: i_non_cas,j,k + integer :: degree + logical :: in_cas + i_non_cas =0 + do k=1,N_det + in_cas = .False. + do j=1,N_det_cas_complex + call get_excitation_degree(psi_cas_complex(1,1,j), psi_det(1,1,k), degree, N_int) + if (degree == 0) then + in_cas = .True. + exit + endif + enddo + if (.not.in_cas) then + double precision :: hij + i_non_cas += 1 + do j=1,N_int + psi_non_cas_complex(j,1,i_non_cas) = psi_det(j,1,k) + psi_non_cas_complex(j,2,i_non_cas) = psi_det(j,2,k) + enddo + do j=1,N_states + psi_non_cas_coef_complex(i_non_cas,j) = psi_coef_complex(k,j) + enddo + idx_non_cas_complex(i_non_cas) = k + endif + enddo + N_det_non_cas_complex = i_non_cas +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas_sorted_bit_complex, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ complex*16, psi_non_cas_coef_sorted_bit_complex, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! |CAS| determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + !TODO: should this be n_det_non_cas_complex? + call sort_dets_by_det_search_key_complex(N_det_cas_complex, psi_non_cas_complex, psi_non_cas_coef_complex, size(psi_non_cas_coef_complex,1), & + psi_non_cas_sorted_bit_complex, psi_non_cas_coef_sorted_bit_complex, N_states) + +END_PROVIDER + + +BEGIN_PROVIDER [complex*16, H_matrix_cas_complex, (N_det_cas_complex,N_det_cas_complex)] + implicit none + integer :: i,j + complex*16 :: hij + do i = 1, N_det_cas_complex + do j = 1, N_det_cas_complex + call i_h_j_complex(psi_cas_complex(1,1,i),psi_cas_complex(1,1,j),N_int,hij) + H_matrix_cas_complex(i,j) = hij + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [complex*16, psi_coef_cas_diagonalized_complex, (N_det_cas_complex,N_states)] +&BEGIN_PROVIDER [double precision, psi_cas_energy_diagonalized_complex, (N_states)] + implicit none + integer :: i,j + double precision, allocatable :: eigenvalues(:) + complex*16, allocatable :: eigenvectors(:,:) + allocate (eigenvectors(size(H_matrix_cas,1),N_det_cas)) + allocate (eigenvalues(N_det_cas)) + call lapack_diag_complex(eigenvalues,eigenvectors, & + H_matrix_cas_complex,size(H_matrix_cas_complex,1),N_det_cas_complex) + do i = 1, N_states + psi_cas_energy_diagonalized_complex(i) = eigenvalues(i) + do j = 1, N_det_cas_complex + psi_coef_cas_diagonalized_complex(j,i) = eigenvectors(j,i) + enddo + enddo + + + END_PROVIDER + diff --git a/src/utils_complex/qp2-pbc-diff.txt b/src/utils_complex/qp2-pbc-diff.txt index ad83db6d..06d74811 100644 --- a/src/utils_complex/qp2-pbc-diff.txt +++ b/src/utils_complex/qp2-pbc-diff.txt @@ -13,12 +13,14 @@ determinants: ezfio_set_determinants_psi_coef_complex_qp_edit? (need ocaml?) psi_coef_{max,min}? save_wavefunction_specified{,_complex} qp_edit save? - + psi_energy_mono_elec diag_h_mat_elem for complex + ... DONE create_excitations build_singly_excited_wavefunction{_complex} + ... -------------------------------------------------------------------------------------