From 7fa6c5f60749ba354b7e526558e73d5494c2bda6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 5 May 2015 11:44:04 +0200 Subject: [PATCH 1/5] Equations in LaTeX format in pseudo --- src/Pseudo_integrals/int.f90 | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/Pseudo_integrals/int.f90 b/src/Pseudo_integrals/int.f90 index be806b3f..cf27b6a9 100644 --- a/src/Pseudo_integrals/int.f90 +++ b/src/Pseudo_integrals/int.f90 @@ -1,7 +1,15 @@ -!! Vps= +!! INFO : You can display equations using : http://www.codecogs.com/latex/eqneditor.php + +!! +!! {\tt Vps}(C) = \langle \Phi_A|{\tt Vloc}(C)+{\tt Vpp}(C)| \Phi_B \rangle +!! +!! with : +!! +!! {\tt Vloc}(C)=\sum_{k=1}^{\tt klocmax} v_k r_C^{n_k} \exp(-dz_k r_C^2) \\ +!! +!! {\tt Vpp}(C)=\sum_{l=0}^{\tt lmax}\left( \sum_{k=1}^{\tt kmax} v_{kl} +!! r_C^{n_{kl}} \exp(-dz_{kl} r_C)^2 \right) |l\rangle \langle l| !! -!! with: Vloc(C)=\sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) -!! Vpp(C)=\sum_{l=0}^lmax\sum_{k=1}^kmax v_kl rC**n_kl exp(-dz_kl rC**2)*|l> Date: Wed, 6 May 2015 15:05:08 +0200 Subject: [PATCH 2/5] Added EZFIO.cfg and s2_eig test in diag --- src/Bitmask/README.rst | 7 ++ src/Bitmask/bitmasks_routines.irp.f | 35 ++++++- src/Determinants/EZFIO.cfg | 2 +- src/Determinants/README.rst | 38 +++----- src/Determinants/determinants.irp.f | 3 - src/Determinants/diagonalize_CI.irp.f | 48 ++++----- src/Determinants/diagonalize_CI_mono.irp.f | 52 ++++++---- src/Determinants/s2.irp.f | 13 ++- src/Determinants/save_for_qmcchem.irp.f | 6 +- src/Determinants/spindeterminants.irp.f | 108 +++++++++++---------- src/MRCC/mrcc_utils.irp.f | 40 +++++--- src/MonoInts/EZFIO.cfg | 6 ++ 12 files changed, 213 insertions(+), 145 deletions(-) create mode 100644 src/MonoInts/EZFIO.cfg diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index 395efc52..7eeac8c9 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -128,6 +128,10 @@ Documentation Subroutine to print the content of a determinant in '+-' notation and hexadecimal representation. +`debug_spindet `_ + Subroutine to print the content of a determinant in '+-' notation and + hexadecimal representation. + `list_to_bitstring `_ Returns the physical string "string(N_int,2)" from the array of occupations "list(N_int*bit_kind_size,2) @@ -135,5 +139,8 @@ Documentation `print_det `_ Subroutine to print the content of a determinant using the '+-' notation +`print_spindet `_ + Subroutine to print the content of a determinant using the '+-' notation + diff --git a/src/Bitmask/bitmasks_routines.irp.f b/src/Bitmask/bitmasks_routines.irp.f index 8694f93f..1b3186b1 100644 --- a/src/Bitmask/bitmasks_routines.irp.f +++ b/src/Bitmask/bitmasks_routines.irp.f @@ -126,7 +126,7 @@ subroutine debug_det(string,Nint) END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint,2) - character*(512) :: output(2) + character*(2048) :: output(2) call bitstring_to_hexa( output(1), string(1,1), Nint ) call bitstring_to_hexa( output(2), string(1,2), Nint ) print *, trim(output(1)) , '|', trim(output(2)) @@ -143,7 +143,7 @@ subroutine print_det(string,Nint) END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint,2) - character*(512) :: output(2) + character*(2048) :: output(2) call bitstring_to_str( output(1), string(1,1), Nint ) call bitstring_to_str( output(2), string(1,2), Nint ) @@ -151,3 +151,34 @@ subroutine print_det(string,Nint) print *, trim(output(2)) end + +subroutine debug_spindet(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant in '+-' notation and + ! hexadecimal representation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + character*(2048) :: output(1) + call bitstring_to_hexa( output(1), string(1,1), Nint ) + print *, trim(output(1)) + call print_spindet(string,Nint) + +end + +subroutine print_spindet(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant using the '+-' notation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint,2) + character*(2048) :: output(1) + + call bitstring_to_str( output(1), string(1,1), Nint ) + print *, trim(output(1)) + +end diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 5f63404b..eab4ea0e 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -100,4 +100,4 @@ size: (determinants_det_num) [expected_s2] interface: OCaml doc: expcted_s2 -type: double precision \ No newline at end of file +type: double precision diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 69f7abc3..bf1e8c13 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -532,7 +532,7 @@ Documentation `save_dets_qmcchem `_ Undocumented -`save_for_qmc `_ +`save_for_qmc `_ Undocumented `save_natorb `_ @@ -623,61 +623,49 @@ Documentation `n_con_int `_ Number of integers to represent the connections between determinants -`create_wf_of_psi_svd_matrix `_ +`create_wf_of_psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`generate_all_alpha_beta_det_products `_ +`generate_all_alpha_beta_det_products `_ Create a wave function from all possible alpha x beta determinants -`get_index_in_psi_det_alpha_unique `_ +`get_index_in_psi_det_alpha_unique `_ Returns the index of the determinant in the ``psi_det_alpha_unique`` array -`get_index_in_psi_det_beta_unique `_ +`get_index_in_psi_det_beta_unique `_ Returns the index of the determinant in the ``psi_det_beta_unique`` array -`n_det_alpha_unique `_ - Unique alpha determinants - -`n_det_beta_unique `_ - Unique beta determinants - `psi_det_alpha `_ List of alpha determinants of psi_det -`psi_det_alpha_unique `_ - Unique alpha determinants - `psi_det_beta `_ List of beta determinants of psi_det -`psi_det_beta_unique `_ - Unique beta determinants - -`psi_svd_alpha `_ +`psi_svd_alpha `_ SVD wave function -`psi_svd_beta `_ +`psi_svd_beta `_ SVD wave function -`psi_svd_coefs `_ +`psi_svd_coefs `_ SVD wave function -`psi_svd_matrix `_ +`psi_svd_matrix `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_columns `_ +`psi_svd_matrix_columns `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_rows `_ +`psi_svd_matrix_rows `_ Matrix of wf coefficients. Outer product of alpha and beta determinants -`psi_svd_matrix_values `_ +`psi_svd_matrix_values `_ Matrix of wf coefficients. Outer product of alpha and beta determinants `spin_det_search_key `_ Return an integer*8 corresponding to a determinant index for searching -`write_spindeterminants `_ +`write_spindeterminants `_ Undocumented `cisd `_ diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index a70d0fe8..8cc545f5 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -147,9 +147,7 @@ END_PROVIDER !$DIR FORCEINLINE bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) enddo - print*,'passed 1' call i8sort(bit_tmp,iorder,N_det) - print*,'passed 2' !DIR$ IVDEP do i=1,N_det do k=1,N_int @@ -189,7 +187,6 @@ END_PROVIDER endif enddo enddo - print*,'passed 3' N_occ_pattern=0 do i=1,N_det diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f index 0e697ab3..c651ab96 100644 --- a/src/Determinants/diagonalize_CI.irp.f +++ b/src/Determinants/diagonalize_CI.irp.f @@ -66,28 +66,32 @@ END_PROVIDER enddo integer :: i_state double precision :: s2 - i_state = 0 - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy(i_state) = eigenvalues(j) - CI_eigenvectors_s2(i_state) = s2 - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo -! if(i_state < min(N_states_diag,N_det))then -! print *, 'pb with the number of states' -! print *, 'i_state = ',i_state -! print *, 'N_states_diag ',N_states_diag -! print *,'stopping ...' -! stop -! endif + if (s2_eig) then + i_state = 0 + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state) = eigenvalues(j) + CI_eigenvectors_s2(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + else + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy(j) = eigenvalues(j) + CI_eigenvectors_s2(j) = s2 + enddo + endif deallocate(eigenvectors,eigenvalues) endif diff --git a/src/Determinants/diagonalize_CI_mono.irp.f b/src/Determinants/diagonalize_CI_mono.irp.f index 1c9a4de3..ad54e430 100644 --- a/src/Determinants/diagonalize_CI_mono.irp.f +++ b/src/Determinants/diagonalize_CI_mono.irp.f @@ -32,25 +32,39 @@ integer :: i_state double precision :: s2 i_state = 0 - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - print*,'j = ',j - print*,'e = ',eigenvalues(j) - print*,'c = ',dabs(eigenvectors(1,j)) - if(dabs(eigenvectors(1,j)).gt.0.9d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_mono(i_state) = eigenvalues(j) - CI_eigenvectors_s2_mono(i_state) = s2 - endif - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo + if (s2_eig) then + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + print*,'j = ',j + print*,'e = ',eigenvalues(j) + print*,'c = ',dabs(eigenvectors(1,j)) + if(dabs(eigenvectors(1,j)).gt.0.9d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_mono(i_state) = eigenvalues(j) + CI_eigenvectors_s2_mono(i_state) = s2 + endif + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + else + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(eigenvectors(1,j)).gt.0.9d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_mono(i_state) = eigenvalues(j) + CI_eigenvectors_s2_mono(i_state) = s2 + endif + enddo + endif deallocate(eigenvectors,eigenvalues) endif diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index cd1d9fda..376528d7 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -92,14 +92,13 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) & !$OMP REDUCTION(+:s2) SCHEDULE(dynamic) - do i = 1, n + do i=1,n call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) -! print*,'s2_tmp = ',s2_tmp - do j = 1, n - call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) - if (s2_tmp == 0.d0) cycle - s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp - enddo + s2 += psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp + do j=i+1,n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) + s2 += (psi_coefs_tmp(i)+psi_coefs_tmp(i))*psi_coefs_tmp(j)*s2_tmp + enddo enddo !$OMP END PARALLEL DO end diff --git a/src/Determinants/save_for_qmcchem.irp.f b/src/Determinants/save_for_qmcchem.irp.f index b707ff7c..93a9778d 100644 --- a/src/Determinants/save_for_qmcchem.irp.f +++ b/src/Determinants/save_for_qmcchem.irp.f @@ -7,8 +7,6 @@ subroutine save_dets_qmcchem integer, allocatable :: occ(:,:,:), occ_tmp(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp - read_wf = .True. - TOUCH read_wf call ezfio_set_determinants_det_num(N_det) call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1)) @@ -46,6 +44,8 @@ subroutine save_dets_qmcchem end program save_for_qmc - call save_dets_qmcchem + read_wf = .True. + TOUCH read_wf +! call save_dets_qmcchem call write_spindeterminants end diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index afe88880..5daf4e2c 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -50,80 +50,88 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta, (N_int,psi_det_size) ] enddo END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha_unique, (N_int,psi_det_size) ] -&BEGIN_PROVIDER [ integer, N_det_alpha_unique ] + +BEGIN_TEMPLATE + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_$alpha_unique, (N_int,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_$alpha_unique ] implicit none BEGIN_DOC - ! Unique alpha determinants + ! Unique $alpha determinants END_DOC - integer :: i,k + integer :: i,j,k integer, allocatable :: iorder(:) integer*8, allocatable :: bit_tmp(:) integer*8 :: last_key integer*8, external :: spin_det_search_key + logical,allocatable :: duplicate(:) - allocate ( iorder(N_det), bit_tmp(N_det)) + allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) ) do i=1,N_det iorder(i) = i - bit_tmp(i) = spin_det_search_key(psi_det_alpha(1,i),N_int) + bit_tmp(i) = spin_det_search_key(psi_det_$alpha(1,i),N_int) enddo call i8sort(bit_tmp,iorder,N_det) - N_det_alpha_unique = 0 + N_det_$alpha_unique = 0 last_key = 0_8 do i=1,N_det - if (bit_tmp(i) /= last_key) then - last_key = bit_tmp(i) - N_det_alpha_unique += 1 - do k=1,N_int - psi_det_alpha_unique(k,N_det_alpha_unique) = psi_det_alpha(k,iorder(i)) - enddo + last_key = bit_tmp(i) + N_det_$alpha_unique += 1 + do k=1,N_int + psi_det_$alpha_unique(k,N_det_$alpha_unique) = psi_det_$alpha(k,iorder(i)) + enddo + duplicate(i) = .False. + enddo + + j=1 + do i=1,N_det_$alpha_unique-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j += 1 + cycle + endif + duplicate(j) = .True. + do k=1,N_int + if (psi_det_$alpha_unique(k,i) /= psi_det_$alpha_unique(k,j)) then + duplicate(j) = .False. + exit + endif + enddo + j+=1 + if (j > N_det_$alpha_unique) then + exit + endif + enddo + enddo + + j=1 + do i=2,N_det_$alpha_unique + if (duplicate(i)) then + cycle + else + j += 1 + psi_det_$alpha_unique(:,j) = psi_det_$alpha_unique(:,i) endif enddo + N_det_$alpha_unique = j - deallocate (iorder, bit_tmp) + deallocate (iorder, bit_tmp, duplicate) END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta_unique, (N_int,psi_det_size) ] -&BEGIN_PROVIDER [ integer, N_det_beta_unique ] - implicit none - BEGIN_DOC - ! Unique beta determinants - END_DOC +SUBST [ alpha ] - integer :: i,k - integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8 :: last_key - integer*8, external :: spin_det_search_key - - allocate ( iorder(N_det), bit_tmp(N_det)) - - do i=1,N_det - iorder(i) = i - bit_tmp(i) = spin_det_search_key(psi_det_beta(1,i),N_int) - enddo - - call i8sort(bit_tmp,iorder,N_det) - - N_det_beta_unique = 0 - last_key = 0_8 - do i=1,N_det - if (bit_tmp(i) /= last_key) then - last_key = bit_tmp(i) - N_det_beta_unique += 1 - do k=1,N_int - psi_det_beta_unique(k,N_det_beta_unique) = psi_det_beta(k,iorder(i)) - enddo - endif - enddo - - deallocate (iorder, bit_tmp) -END_PROVIDER +alpha ;; +beta ;; +END_TEMPLATE @@ -150,6 +158,7 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) !DIR$ FORCEINLINE det_ref = spin_det_search_key(key,Nint) + !DIR$ FORCEINLINE det_search = spin_det_search_key(psi_det_alpha_unique(1,1),Nint) @@ -439,6 +448,7 @@ BEGIN_PROVIDER [ double precision, psi_svd_matrix_values, (N_det,N_states) ] do k=1,N_det i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int) j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int) + do l=1,N_states psi_svd_matrix_values(k,l) = psi_coef(k,l) enddo diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 9b4add38..3c250155 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -110,20 +110,32 @@ END_PROVIDER integer :: i_state double precision :: s2 i_state = 0 - do j=1,N_det - call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) - if(dabs(s2-expected_s2).le.0.3d0)then - i_state += 1 - do i=1,N_det - CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state) = s2 - endif - if (i_state.ge.N_states_diag) then - exit - endif - enddo + if (s2_eig) then + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + else + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + i_state += 1 + do i=1,N_det + CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state) = s2 + enddo + endif deallocate(eigenvectors,eigenvalues) endif diff --git a/src/MonoInts/EZFIO.cfg b/src/MonoInts/EZFIO.cfg new file mode 100644 index 00000000..33a8d619 --- /dev/null +++ b/src/MonoInts/EZFIO.cfg @@ -0,0 +1,6 @@ +[do_pseudo] +type: logical +doc: Using pseudo potential integral of not +interface: input +default: False + From d6ec67e8408ab2408dafd2db1d167b4503647184 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 May 2015 17:01:45 +0200 Subject: [PATCH 3/5] Fixed Full CAS-SD --- src/CAS_SD/cas_sd.irp.f | 14 ----- src/Determinants/H_apply.irp.f | 95 ++++++++++++++++++++++---------- src/Determinants/README.rst | 40 +++++++------- src/Determinants/s2.irp.f | 3 +- src/MRCC/README.rst | 4 +- src/MonoInts/pot_ao_ints.irp.f | 2 +- src/Perturbation/selection.irp.f | 2 +- 7 files changed, 92 insertions(+), 68 deletions(-) diff --git a/src/CAS_SD/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f index 144eec88..0b34e7bd 100644 --- a/src/CAS_SD/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -69,18 +69,4 @@ program full_ci print *, '' enddo print *, 'Max excitation degree in the CAS :', exc_max - do i=1,N_det - in_cas = .False. - degree_min = 1000 - do k=1,N_det_generators - call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int) - degree_min = min(degree_min,degree) - enddo - if (degree_min > 2) then - print *, 'Error : This is not a CAS-SD : ' - print *, 'Excited determinant:', degree_min - call debug_det(psi_det(1,1,k),N_int) - stop - endif - enddo end diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index d230c765..dab4390a 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -53,13 +53,17 @@ subroutine resize_H_apply_buffer(new_size,iproc) double precision, pointer :: buffer_e2(:,:) integer :: i,j,k integer :: Ndet + + BEGIN_DOC +! Resizes the H_apply buffer of proc iproc. The buffer lock should +! be set before calling this function. + END_DOC PROVIDE H_apply_buffer_allocated ASSERT (new_size > 0) ASSERT (iproc >= 0) ASSERT (iproc < nproc) - call omp_set_lock(H_apply_buffer_lock(1,iproc)) allocate ( buffer_det(N_int,2,new_size), & buffer_coef(new_size,N_states), & buffer_e2(new_size,N_states) ) @@ -93,7 +97,6 @@ subroutine resize_H_apply_buffer(new_size,iproc) H_apply_buffer(iproc)%sze = new_size H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det) - call omp_unset_lock(H_apply_buffer_lock(1,iproc)) end @@ -101,8 +104,7 @@ subroutine copy_H_apply_buffer_to_wf use omp_lib implicit none BEGIN_DOC -! Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det -! after calling this function. +! Copies the H_apply buffer to psi_coef. ! After calling this subroutine, N_det, psi_det and psi_coef need to be touched END_DOC integer(bit_kind), allocatable :: buffer_det(:,:,:) @@ -181,42 +183,77 @@ subroutine copy_H_apply_buffer_to_wf call normalize(psi_coef,N_det) SOFT_TOUCH N_det psi_det psi_coef - call debug_unicity_of_determinants + logical :: found_duplicates + call remove_duplicates_in_psi_det(found_duplicates) end -subroutine debug_unicity_of_determinants +subroutine remove_duplicates_in_psi_det(found_duplicates) implicit none + logical, intent(out) :: found_duplicates BEGIN_DOC -! This subroutine checks that there are no repetitions in the wave function +! Removes duplicate determinants in the wave function. END_DOC - logical :: same, failed - integer :: i,k - print *, "======= DEBUG UNICITY =========" - failed = .False. - do i=2,N_det - same = .True. - do k=1,N_int - if ( psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,i-1) ) then - same = .False. - exit + integer :: i,j,k + integer(bit_kind), allocatable :: bit_tmp(:) + logical,allocatable :: duplicate(:) + + allocate (duplicate(N_det), bit_tmp(N_det)) + + do i=1,N_det + integer, external :: det_search_key + !$DIR FORCEINLINE + bit_tmp(i) = det_search_key(psi_det_sorted_bit(1,1,i),N_int) + duplicate(i) = .False. + enddo + + do i=1,N_det-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j += 1 + cycle endif - if ( psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,i-1) ) then - same = .False. + duplicate(j) = .True. + do k=1,N_int + if ( (psi_det_sorted_bit(k,1,i) /= psi_det_sorted_bit(k,1,j) ) & + .or. (psi_det_sorted_bit(k,2,i) /= psi_det_sorted_bit(k,2,j) ) ) then + duplicate(j) = .False. + exit + endif + enddo + j += 1 + if (j > N_det) then exit endif enddo - if (same) then - failed = .True. - call debug_det(psi_det_sorted_bit(1,1,i)) + enddo + + found_duplicates = .False. + do i=1,N_det + if (duplicate(i)) then + found_duplicates = .True. + exit endif enddo - if (failed) then - print *, '======= Determinants not unique : Failed ! =========' - stop - else - print *, '======= Determinants are unique : OK ! =========' + if (found_duplicates) then + call write_bool(output_determinants,found_duplicates,'Found duplicate determinants') + k=0 + do i=1,N_det + if (.not.duplicate(i)) then + k += 1 + psi_det(:,:,k) = psi_det_sorted_bit (:,:,i) + psi_coef(k,:) = psi_coef_sorted_bit(i,:) + endif + enddo + print *, N_det,k + N_det = k + TOUCH N_det psi_det psi_coef endif + deallocate (duplicate,bit_tmp) end subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) @@ -231,11 +268,11 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) integer :: i,j,k integer :: new_size PROVIDE H_apply_buffer_allocated + call omp_set_lock(H_apply_buffer_lock(1,iproc)) new_size = H_apply_buffer(iproc)%N_det + n_selected if (new_size > H_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc) endif - call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) @@ -250,7 +287,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) enddo do j=1,N_states do i=1,N_selected - H_apply_buffer(iproc)%coef(i,j) = 0.d0 + H_apply_buffer(iproc)%coef(i+H_apply_buffer(iproc)%N_det,j) = 0.d0 enddo enddo H_apply_buffer(iproc)%N_det = new_size diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index bf1e8c13..4b8c0472 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -187,10 +187,10 @@ Documentation `det_svd `_ Computes the SVD of the Alpha x Beta determinant coefficient matrix -`filter_3_highest_electrons `_ +`filter_3_highest_electrons `_ Returns a determinant with only the 3 highest electrons -`int_of_3_highest_electrons `_ +`int_of_3_highest_electrons `_ Returns an integer*8 as : .br |_<--- 21 bits ---><--- 21 bits ---><--- 21 bits --->| @@ -207,26 +207,26 @@ Documentation `n_det `_ Number of determinants in the wave function -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_ab `_ +`psi_coef_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave @@ -239,46 +239,46 @@ Documentation `psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_ab `_ +`psi_det_sorted_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave function. -`psi_det_sorted_next_ab `_ +`psi_det_sorted_next_ab `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`save_wavefunction_general `_ +`save_wavefunction_general `_ Save the wave function into the EZFIO file -`save_wavefunction_unsorted `_ +`save_wavefunction_unsorted `_ Save the wave function into the EZFIO file -`sort_dets_by_3_highest_electrons `_ +`sort_dets_by_3_highest_electrons `_ Determinants on which we apply . They are sorted by the 3 highest electrons in the alpha part, then by the 3 highest electrons in the beta part to accelerate the research of connected determinants. -`sort_dets_by_det_search_key `_ +`sort_dets_by_det_search_key `_ Determinants are sorted are sorted according to their det_search_key. Useful to accelerate the search of a random determinant in the wave function. @@ -316,7 +316,7 @@ Documentation `diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix @@ -345,7 +345,7 @@ Documentation `ci_electronic_energy_mono `_ Eigenvectors/values of the CI matrix -`diagonalize_ci_mono `_ +`diagonalize_ci_mono `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index 376528d7..760893e0 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -88,7 +88,7 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) double precision, intent(out) :: s2 integer :: i,j,l double precision :: s2_tmp - s2 = S_z2_Sz + s2 = 0.d0 !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) & !$OMP REDUCTION(+:s2) SCHEDULE(dynamic) @@ -101,5 +101,6 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) enddo enddo !$OMP END PARALLEL DO + s2 += S_z2_Sz end diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 38137667..5bb2009c 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -51,7 +51,7 @@ Documentation `ci_electronic_energy_dressed `_ Eigenvectors/values of the CI matrix -`ci_energy_dressed `_ +`ci_energy_dressed `_ N_states lowest eigenvalues of the dressed CI matrix `delta_ij `_ @@ -60,7 +60,7 @@ Documentation `delta_ij_non_cas `_ Dressing matrix in SD basis -`diagonalize_ci_dressed `_ +`diagonalize_ci_dressed `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index c93301e5..a9e72e57 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -10,7 +10,7 @@ integer :: i,j,k,l,n_pt_in,m double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - if (do_pseudo.eqv..TRUE.) then + if (do_pseudo) then ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo else ao_nucl_elec_integral = 0.d0 diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 77313888..84cc59ae 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -19,6 +19,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c ASSERT (Nint > 0) ASSERT (N_int == N_int) ASSERT (N_selected >= 0) + call omp_set_lock(H_apply_buffer_lock(1,iproc)) smax = selection_criterion smin = selection_criterion_min new_size = H_apply_buffer(iproc)%N_det + n_selected @@ -26,7 +27,6 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c if (new_size > h_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(h_apply_buffer(iproc)%sze*2,new_size),iproc) endif - call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) From 8ee7b87954a7e7caee8431eeb5c3d6c9423260c8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 May 2015 17:02:59 +0200 Subject: [PATCH 4/5] Removed print --- src/Determinants/H_apply.irp.f | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index dab4390a..6e2b3a5a 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -249,7 +249,6 @@ subroutine remove_duplicates_in_psi_det(found_duplicates) psi_coef(k,:) = psi_coef_sorted_bit(i,:) endif enddo - print *, N_det,k N_det = k TOUCH N_det psi_det psi_coef endif From 0a1ba724aa51983121b6c7f15262e08295b05c7b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 May 2015 18:26:36 +0200 Subject: [PATCH 5/5] Fixed bugs --- src/CAS_SD/cas_sd.irp.f | 8 +++-- src/CAS_SD/cas_sd_selected.irp.f | 36 ++++++++++++---------- src/Determinants/README.rst | 16 +++++----- src/Determinants/diagonalize_CI.irp.f | 2 +- src/Determinants/diagonalize_CI_mono.irp.f | 2 +- src/MRCC/mrcc_utils.irp.f | 2 +- 6 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/CAS_SD/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f index 0b34e7bd..e28334b2 100644 --- a/src/CAS_SD/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -1,6 +1,7 @@ program full_ci implicit none integer :: i,k + integer :: N_det_old double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) @@ -9,6 +10,7 @@ program full_ci allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) character*(64) :: perturbation + N_det_old = 0 pt2 = 1.d0 diag_algorithm = "Lapack" if (N_det > n_det_max_cas_sd) then @@ -29,6 +31,7 @@ program full_ci endif do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) + N_det_old = N_det call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef @@ -53,10 +56,11 @@ program full_ci if (abort_all) then exit endif + if (N_det == N_det_old) then + exit + endif enddo - ! Check that it is a CAS-SD - logical :: in_cas integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_generators diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f index 86cac969..c3abadf4 100644 --- a/src/CAS_SD/cas_sd_selected.irp.f +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -50,14 +50,30 @@ program full_ci print *, 'E = ', CI_energy print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' - call ezfio_set_full_ci_energy(CI_energy) + call ezfio_set_cas_sd_energy(CI_energy(1)) if (abort_all) then exit endif enddo + call diagonalize_CI + + if(do_pt2_end)then + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + call H_apply_CAS_SD_PT2(pt2, norm_pert, H_pert_diag, N_st) + + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_cas_sd_energy_pt2(CI_energy(1)+pt2(1)) + endif + - ! Check that it is a CAS-SD - logical :: in_cas integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_cas @@ -70,18 +86,4 @@ program full_ci print *, '' enddo print *, 'Max excitation degree in the CAS :', exc_max - do i=1,N_det - in_cas = .False. - degree_min = 1000 - do k=1,N_det_cas - call get_excitation_degree(psi_cas(1,1,k),psi_det(1,1,i),degree,N_int) - degree_min = min(degree_min,degree) - enddo - if (degree_min > 2) then - print *, 'Error : This is not a CAS-SD : ' - print *, 'Excited determinant:', degree_min - call debug_det(psi_det(1,1,k),N_int) - stop - endif - enddo end diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index 4b8c0472..e9307357 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -40,15 +40,11 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ - Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det - after calling this function. +`copy_h_apply_buffer_to_wf `_ + Copies the H_apply buffer to psi_coef. After calling this subroutine, N_det, psi_det and psi_coef need to be touched -`debug_unicity_of_determinants `_ - This subroutine checks that there are no repetitions in the wave function - -`fill_h_apply_buffer_no_selection `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD `h_apply_buffer_allocated `_ @@ -59,8 +55,12 @@ Documentation Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. +`remove_duplicates_in_psi_det `_ + Removes duplicate determinants in the wave function. + `resize_h_apply_buffer `_ - Undocumented + Resizes the H_apply buffer of proc iproc. The buffer lock should + be set before calling this function. `cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f index c651ab96..7c017956 100644 --- a/src/Determinants/diagonalize_CI.irp.f +++ b/src/Determinants/diagonalize_CI.irp.f @@ -83,7 +83,7 @@ END_PROVIDER endif enddo else - do j=1,N_det + do j=1,N_states_diag call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) do i=1,N_det CI_eigenvectors(i,j) = eigenvectors(i,j) diff --git a/src/Determinants/diagonalize_CI_mono.irp.f b/src/Determinants/diagonalize_CI_mono.irp.f index ad54e430..3f9b94ec 100644 --- a/src/Determinants/diagonalize_CI_mono.irp.f +++ b/src/Determinants/diagonalize_CI_mono.irp.f @@ -53,7 +53,7 @@ endif enddo else - do j=1,N_det + do j=1,N_states_diag call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) if(dabs(eigenvectors(1,j)).gt.0.9d0)then i_state += 1 diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 3c250155..611a70be 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -126,7 +126,7 @@ END_PROVIDER endif enddo else - do j=1,N_det + do j=1,N_states_diag call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) i_state += 1 do i=1,N_det