From 583e8859d0602841d430b360b02a0b1c22878878 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 25 May 2014 01:18:41 +0200 Subject: [PATCH] CISD is cleaned. Perturbation/selection is broken --- scripts/generate_h_apply.py | 55 ++++++++++++- src/Bitmask/README.rst | 6 ++ src/Bitmask/bitmasks.ezfio_config | 4 +- src/Bitmask/bitmasks.irp.f | 93 ++++++++++++++++++++-- src/Bitmask/bitmasks_module.f90 | 6 ++ src/CISD/H_apply.irp.f | 70 ++--------------- src/CISD/README.rst | 21 +++-- src/CISD/cisd.irp.f | 42 ++-------- src/CISD/cisd_lapack.irp.f | 28 +++---- src/CISD/cisd_sc2.irp.f | 57 ++++---------- src/CISD/h_apply.py | 11 --- src/Dets/H_apply.irp.f | 41 +++++++++- src/Dets/H_apply_template.f | 99 ++++++++++++++++-------- src/Dets/README.rst | 3 + src/Dets/determinants.ezfio_config | 6 +- src/Dets/determinants.irp.f | 37 ++++++++- src/Dets/diagonalize_CI.irp.f | 64 +++++++++++++++ src/Dets/utils.irp.f | 3 +- src/Perturbation/README.rst | 30 +++---- src/Perturbation/epstein_nesbet.irp.f | 7 +- src/Perturbation/perturbation_template.f | 2 +- src/Perturbation/perturbation_test.irp.f | 6 -- src/Perturbation/temporary_stuff.irp.f | 72 ----------------- src/Utils/LinearAlgebra.irp.f | 52 ++++++++----- 24 files changed, 476 insertions(+), 339 deletions(-) delete mode 100755 src/CISD/h_apply.py create mode 100644 src/Dets/diagonalize_CI.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index d6187b3e..ab43f743 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -8,10 +8,14 @@ file.close() keywords = """ subroutine parameters +params_main initialization declarations +decls_main keys_work +copy_buffer finalization +generate_psi_guess """.split() class H_apply(object): @@ -21,9 +25,8 @@ class H_apply(object): for k in keywords: s[k] = "" s["subroutine"] = "H_apply_%s"%(sub) + s["params_post"] = "" self.openmp = openmp - if openmp: - s["subroutine"] += "_OpenMP" self.selection_pt2 = None self.perturbation = None @@ -47,6 +50,33 @@ class H_apply(object): s["omp_do"] = "!$OMP DO SCHEDULE (static)" s["omp_enddo"] = "!$OMP ENDDO NOWAIT" + s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)" + + s["generate_psi_guess"] = """ + ! Sort H_jj to find the N_states lowest states + integer :: i,k + integer, allocatable :: iorder(:) + double precision, allocatable :: H_jj(:) + double precision, external :: diag_h_mat_elem + allocate(H_jj(N_det),iorder(N_det)) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(psi_det,N_int,H_jj,iorder,N_det) & + !$OMP PRIVATE(i) + !$OMP DO + do i = 1, N_det + H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) + iorder(i) = i + enddo + !$OMP END DO + !$OMP END PARALLEL + + call dsort(H_jj,iorder,N_det) + do k=1,N_states + psi_coef(iorder(k),k) = 1.d0 + enddo + deallocate(H_jj,iorder) + """ + if not openmp: for k in s: s[k] = "" @@ -54,6 +84,7 @@ class H_apply(object): s["set_i_H_j_threshold"] = """ thresh = H_apply_threshold """ + s["copy_buffer"] = "call copy_h_apply_buffer_to_wf" self.data = s def __setitem__(self,key,value): @@ -101,6 +132,21 @@ class H_apply(object): sum_norm_pert_in = sum_norm_pert sum_H_pert_diag_in = sum_H_pert_diag """ + self.data["copy_buffer"] = "" + self.data["generate_psi_guess"] = "" + + self.data["params_main"] = "pt2, norm_pert, H_pert_diag, N_st" + self.data["params_post"] = ","+self.data["params_main"] +", N_int" + self.data["decls_main"] = """ integer, intent(in) :: N_st + double precision, intent(inout):: pt2(N_st) + double precision, intent(inout):: norm_pert(N_st) + double precision, intent(inout):: H_pert_diag + PROVIDE reference_energy N_det_generators key_pattern_not_in_ref + pt2 = 0.d0 + norm_pert = 0.d0 + H_pert_diag = 0.d0 + """ + if self.openmp: self.data["omp_parallel"] += """& !$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) & @@ -113,6 +159,11 @@ class H_apply(object): self.selection_pt2 = pert if pert is not None: self.data["size_max"] = str(1024*128) + self.data["copy_buffer"] = """ + call copy_h_apply_buffer_to_wf + selection_criterion_min = selection_criterion_min*0.1d0 + selection_criterion = selection_criterion_min + """ self.data["keys_work"] = """ e_2_pert_buffer = 0.d0 coef_pert_buffer = 0.d0 diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index ada78e34..772c5b1c 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -69,9 +69,15 @@ Documentation `n_int `_ Number of 64-bit integers needed to represent determinants as binary strings +`n_reference_bitmask `_ + Number of bitmasks for reference + `ref_bitmask `_ Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask +`reference_bitmask `_ + Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference) + `bitstring_to_hexa `_ Transform a bit string to a string in hexadecimal format for printing diff --git a/src/Bitmask/bitmasks.ezfio_config b/src/Bitmask/bitmasks.ezfio_config index f402ef02..79f1ae6c 100644 --- a/src/Bitmask/bitmasks.ezfio_config +++ b/src/Bitmask/bitmasks.ezfio_config @@ -2,7 +2,7 @@ bitmasks N_int integer bit_kind integer N_mask_gen integer - generators integer (bitmasks_N_int*bitmasks_bit_kind/4,2,2,bitmasks_N_mask_gen) + generators integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen) N_mask_ref integer - reference integer (bitmasks_N_int*bitmasks_bit_kind/4,2,2,bitmasks_N_mask_ref) + reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_ref) diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 3aa71421..71833017 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -88,10 +88,17 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ] END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,2,N_generators_bitmask) ] +BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_bitmask) ] implicit none BEGIN_DOC - ! Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator) + ! Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator). + ! 3rd index is : + ! * 1 : hole for single exc + ! * 1 : particle for single exc + ! * 3 : hole for 1st exc of double + ! * 4 : particle for 1st exc of double + ! * 5 : hole for 2dn exc of double + ! * 6 : particle for 2dn exc of double END_DOC logical :: exists PROVIDE ezfio_filename @@ -100,12 +107,86 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,2,N_generators_ if (exists) then call ezfio_get_bitmasks_generators(generators_bitmask) else - generators_bitmask(:,:,1,1) = HF_bitmask - generators_bitmask(:,1,2,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1)) - generators_bitmask(:,2,2,1) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2)) + generators_bitmask(:,:,s_hole ,1) = HF_bitmask + generators_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + generators_bitmask(:,:,d_hole1,1) = HF_bitmask + generators_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + generators_bitmask(:,:,d_hole2,1) = HF_bitmask + generators_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) call ezfio_set_bitmasks_generators(generators_bitmask) endif - ASSERT (N_generators_bitmask > 0) END_PROVIDER +BEGIN_PROVIDER [ integer, N_reference_bitmask ] + implicit none + BEGIN_DOC + ! Number of bitmasks for reference + END_DOC + logical :: exists + PROVIDE ezfio_filename + + call ezfio_has_bitmasks_N_mask_ref(exists) + if (exists) then + call ezfio_get_bitmasks_N_mask_ref(N_reference_bitmask) + integer :: N_int_check + integer :: bit_kind_check + call ezfio_get_bitmasks_bit_kind(bit_kind_check) + if (bit_kind_check /= bit_kind) then + print *, bit_kind_check, bit_kind + print *, 'Error: bit_kind is not correct in EZFIO file' + endif + call ezfio_get_bitmasks_N_int(N_int_check) + if (N_int_check /= N_int) then + print *, N_int_check, N_int + print *, 'Error: N_int is not correct in EZFIO file' + endif + else + N_reference_bitmask = 1 + call ezfio_set_bitmasks_N_int(N_int) + call ezfio_set_bitmasks_bit_kind(bit_kind) + call ezfio_set_bitmasks_N_mask_ref(N_reference_bitmask) + endif + ASSERT (N_reference_bitmask > 0) + +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), reference_bitmask, (N_int,2,2,N_reference_bitmask) ] + implicit none + BEGIN_DOC + ! Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference) + END_DOC + logical :: exists + PROVIDE ezfio_filename + + call ezfio_has_bitmasks_reference(exists) + if (exists) then + call ezfio_get_bitmasks_reference(reference_bitmask) + else + reference_bitmask(:,:,s_hole ,1) = HF_bitmask + reference_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + reference_bitmask(:,:,d_hole1,1) = HF_bitmask + reference_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + reference_bitmask(:,:,d_hole2,1) = HF_bitmask + reference_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + call ezfio_set_bitmasks_reference(reference_bitmask) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, i_bitmask_gen ] + implicit none + BEGIN_DOC + ! Current bitmask for the generators + END_DOC + i_bitmask_gen = 1 +END_PROVIDER + +BEGIN_PROVIDER [ integer, i_bitmask_ref ] + implicit none + BEGIN_DOC + ! Current bitmask for the reference + END_DOC + i_bitmask_ref = 1 +END_PROVIDER + diff --git a/src/Bitmask/bitmasks_module.f90 b/src/Bitmask/bitmasks_module.f90 index 28363ff0..5e5e2f02 100644 --- a/src/Bitmask/bitmasks_module.f90 +++ b/src/Bitmask/bitmasks_module.f90 @@ -2,4 +2,10 @@ module bitmasks integer, parameter :: bit_kind_shift = 6 ! 5: 32 bits, 6: 64 bits integer, parameter :: bit_kind_size = 64 integer, parameter :: bit_kind = 64/8 + integer, parameter :: s_hole = 1 + integer, parameter :: s_part = 2 + integer, parameter :: d_hole1 = 3 + integer, parameter :: d_part1 = 4 + integer, parameter :: d_hole2 = 5 + integer, parameter :: d_part2 = 6 end module diff --git a/src/CISD/H_apply.irp.f b/src/CISD/H_apply.irp.f index 15582375..a7dd568f 100644 --- a/src/CISD/H_apply.irp.f +++ b/src/CISD/H_apply.irp.f @@ -1,65 +1,9 @@ -BEGIN_SHELL [ /bin/bash ] -./h_apply.py +! Generates subroutine H_apply_cisd +! ---------------------------------- + +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import H_apply +H = H_apply("cisd",openmp=True) +print H END_SHELL - -subroutine fill_H_apply_buffer(n_selected,det_buffer,Nint,iproc) - use bitmasks - implicit none - BEGIN_DOC - ! Fill the H_apply buffer with determiants for CISD - END_DOC - - integer, intent(in) :: n_selected, Nint, iproc - integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) - integer :: i,j,k - integer :: new_size - PROVIDE H_apply_buffer_allocated - 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 - 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) - enddo - do i=1,n_selected - do j=1,N_int - H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i) - H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i) - enddo - ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num) - ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num) - H_apply_buffer(iproc)%coef(i,:) = 0.d0 - enddo - H_apply_buffer(iproc)%N_det = new_size - 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) - enddo -end - - -subroutine H_apply_cisd - implicit none - BEGIN_DOC - ! Calls H_apply on the HF determinant and selects all connected single and double - ! excitations (of the same symmetry). - END_DOC - - integer(bit_kind) :: hole_mask(N_int,2) - integer(bit_kind) :: particle_mask(N_int,2) - - ASSERT (N_det_generators == 1) - PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map - - call H_apply_cisd_OpenMP_monoexc(HF_bitmask, & - generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1)) - call H_apply_cisd_OpenMP_diexc(HF_bitmask, & - generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1), & - generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1) ) - - call copy_h_apply_buffer_to_wf -end - - diff --git a/src/CISD/README.rst b/src/CISD/README.rst index a49b2be1..5d4c9b4d 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -1,3 +1,14 @@ +CISD +==== + +This is a test directory which builds a CISD by setting the follwoing rules: + +* The only generator determinant is the Hartee-Fock (single-reference method) +* All generated determinants are included in the wave function (no perturbative + selection) + +These rules are set in the ``H_apply.irp.f`` file. + Needed Modules ============== @@ -24,13 +35,6 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fill_h_apply_buffer `_ - Fill the H_apply buffer with determiants for CISD - -`h_apply_cisd `_ - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). - `cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br @@ -50,5 +54,8 @@ Documentation `repeat_excitation `_ Undocumented +`cisd `_ + Undocumented + diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index e6e7adcd..e46b8bf8 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -1,46 +1,18 @@ program cisd implicit none - integer :: i,k - double precision, allocatable :: eigvalues(:),eigvectors(:,:) - PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map + integer :: i - double precision :: pt2(10), norm_pert(10), H_pert_diag - double precision,allocatable :: H_jj(:) - double precision :: diag_h_mat_elem - integer,allocatable :: iorder(:) - -! N_states = 3 -! touch N_states - call H_apply_cisd - allocate(eigvalues(n_states),eigvectors(n_det,n_states)) - print *, 'N_det = ', N_det + print *, 'HF = ', HF_energy print *, 'N_states = ', N_states - psi_coef = - 1.d-4 - allocate(H_jj(n_det),iorder(n_det)) - do i = 1, N_det - H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) - iorder(i) = i - enddo - call dsort(H_jj,iorder,n_det) - - do k=1,N_states - psi_coef(iorder(k),k) = 1.d0 - enddo - call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD) - do i = 1, N_states - print*,'eigvalues(i) = ',eigvalues(i) + call H_apply_cisd + print *, 'N_det = ', N_det + do i = 1,N_states + print *, 'energy = ',CI_energy(i) + print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy enddo - print *, '---' - print *, 'HF:', HF_energy - print *, '---' - do i = 1,1 - print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion - print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy - enddo ! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) ! do i = 1, N_states ! print*,'eigvalues(i) = ',eigvalues(i) ! enddo - deallocate(eigvalues,eigvectors) end diff --git a/src/CISD/cisd_lapack.irp.f b/src/CISD/cisd_lapack.irp.f index b8f2295b..c57d11f1 100644 --- a/src/CISD/cisd_lapack.irp.f +++ b/src/CISD/cisd_lapack.irp.f @@ -1,21 +1,21 @@ program cisd implicit none integer :: i - PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map + + N_states = 3 + diag_algorithm = "Lapack" + touch N_states diag_algorithm + print *, 'HF = ', HF_energy + print *, 'N_states = ', N_states call H_apply_cisd - double precision, allocatable :: eigvalues(:),eigvectors(:,:) - allocate(eigvalues(n_det),eigvectors(n_det,n_det)) print *, 'N_det = ', N_det - call lapack_diag(eigvalues,eigvectors,H_matrix_all_dets,n_det,n_det) - -! print *, H_matrix_all_dets - print *, '---' - print *, 'HF:', HF_energy - print *, '---' - do i = 1,20 - print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion + do i = 1,N_states + print *, 'energy = ',CI_energy(i) + print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy enddo -! print *, eigvectors(:,1) - deallocate(eigvalues,eigvectors) -end +! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) +! do i = 1, N_states +! print*,'eigvalues(i) = ',eigvalues(i) +! enddo +end diff --git a/src/CISD/cisd_sc2.irp.f b/src/CISD/cisd_sc2.irp.f index 2012855d..6e4b1bd3 100644 --- a/src/CISD/cisd_sc2.irp.f +++ b/src/CISD/cisd_sc2.irp.f @@ -1,46 +1,23 @@ program cisd implicit none - integer :: i,k - double precision, allocatable :: eigvalues(:),eigvectors(:,:) - PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map + integer :: i, j - double precision :: pt2(10), norm_pert(10), H_pert_diag - double precision,allocatable :: H_jj(:) - double precision :: diag_h_mat_elem - integer,allocatable :: iorder(:) - -! N_states = 3 -! touch N_states - call H_apply_cisd - allocate(eigvalues(n_states),eigvectors(n_det,n_states)) - print *, 'N_det = ', N_det + print *, 'HF = ', HF_energy print *, 'N_states = ', N_states - psi_coef = - 1.d-4 - allocate(H_jj(n_det),iorder(n_det)) - do i = 1, N_det - H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) - iorder(i) = i - enddo - call dsort(H_jj,iorder,n_det) + call H_apply_cisd + print *, 'N_det = ', N_det + do i = 1,N_states + print *, 'energy = ',CI_energy(i) + print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy + do j=1,N_det + psi_coef(j,i) = CI_eigenvectors(j,i) + enddo + enddo + SOFT_TOUCH CI_electronic_energy CI_eigenvectors + call CISD_SC2(psi_det,psi_coef,CI_electronic_energy,size(psi_coef,1),N_det,N_states,N_int) + TOUCH CI_electronic_energy CI_eigenvectors - do k=1,N_states - psi_coef(iorder(k),k) = 1.d0 - enddo - call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD) - do i = 1, N_states - print*,'eigvalues(i) = ',eigvalues(i) - enddo - - print *, '---' - print *, 'HF:', HF_energy - print *, '---' - do i = 1,1 - print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion - print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy - enddo - call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD) - do i = 1, N_states - print*,'eigvalues(i) = ',eigvalues(i) - enddo - deallocate(eigvalues,eigvectors) + do i = 1, N_states + print*,'eigvalues(i) = ',CI_energy(i) + enddo end diff --git a/src/CISD/h_apply.py b/src/CISD/h_apply.py deleted file mode 100755 index ac71a5e1..00000000 --- a/src/CISD/h_apply.py +++ /dev/null @@ -1,11 +0,0 @@ -#!/usr/bin/env python - -from generate_h_apply import * - -s = H_apply("cisd",openmp=True) -s["keys_work"] += """ -call fill_H_apply_buffer(key_idx,keys_out,N_int,iproc) -""" -print s - - diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index f3896623..b8a551b5 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -194,8 +194,45 @@ subroutine copy_H_apply_buffer_to_wf end - - +subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + BEGIN_DOC + ! Fill the H_apply buffer with determiants for CISD + END_DOC + + integer, intent(in) :: n_selected, Nint, iproc + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k + integer :: new_size + PROVIDE H_apply_buffer_allocated + 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 + 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) + enddo + do i=1,n_selected + do j=1,N_int + H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i) + H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i) + enddo + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num) + enddo + do j=1,N_states + do i=1,N_selected + H_apply_buffer(iproc)%coef(i,j) = 0.d0 + enddo + enddo + H_apply_buffer(iproc)%N_det = new_size + 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) + enddo +end diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 63044da3..fe9d7394 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -37,7 +37,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame iproc = 0 $omp_parallel -!$ iproc = omp_get_thread_num() + !$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max),hij_tab(size_max)) !print*,'key_in !!' @@ -69,7 +69,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame i_a = occ_hole(ii,ispin) ASSERT (i_a > 0) ASSERT (i_a <= mo_tot_num) - + do jj=1,N_elec_in_key_part_1(ispin) !particle j_a = occ_particle(jj,ispin) ASSERT (j_a > 0) @@ -142,20 +142,20 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) -! call i_H_j(key,key_in,N_int,hij_elec) -! if(dabs(hij_elec)>=thresh) then - key_idx += 1 - do k=1,N_int - keys_out(k,1,key_idx) = key(k,1) - keys_out(k,2,key_idx) = key(k,2) - enddo - hij_tab(key_idx) = hij_elec - ASSERT (key_idx <= size_max) - if (key_idx == size_max) then - $keys_work - key_idx = 0 - endif -! endif + ! call i_H_j(key,key_in,N_int,hij_elec) + ! if(dabs(hij_elec)>=thresh) then + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + hij_tab(key_idx) = hij_elec + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + $keys_work + key_idx = 0 + endif + ! endif enddo enddo endif @@ -181,20 +181,20 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame key(k,ispin) = ibset(key(k,ispin),l) !! a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin) -! call i_H_j(key,key_in,N_int,hij_elec) -! if(dabs(hij_elec)>=thresh) then - key_idx += 1 - do k=1,N_int - keys_out(k,1,key_idx) = key(k,1) - keys_out(k,2,key_idx) = key(k,2) - enddo - hij_tab(key_idx) = hij_elec - ASSERT (key_idx <= size_max) - if (key_idx == size_max) then - $keys_work - key_idx = 0 - endif -! endif + ! call i_H_j(key,key_in,N_int,hij_elec) + ! if(dabs(hij_elec)>=thresh) then + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + hij_tab(key_idx) = hij_elec + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + $keys_work + key_idx = 0 + endif + ! endif enddo enddo! kk enddo ! ii @@ -246,7 +246,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) iproc = 0 $omp_parallel -!$ iproc = omp_get_thread_num() + !$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max),hij_tab(size_max)) !!!! First couple hole particle do j = 1, N_int @@ -255,13 +255,13 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1)) particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2)) enddo - + call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int) call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int) call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int) call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int) allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2)) - + do ispin=1,2 i=0 do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole @@ -315,3 +315,36 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) end + +subroutine $subroutine($params_main) + implicit none + use bitmasks + BEGIN_DOC + ! Calls H_apply on the HF determinant and selects all connected single and double + ! excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + END_DOC + + $decls_main + + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators + integer :: imask + + do imask=1,N_det_generators + call $subroutine_monoexc(psi_generators(1,1,imask), & + generators_bitmask(1,1,s_hole ,i_bitmask_gen), & + generators_bitmask(1,1,s_part ,i_bitmask_gen) & + $params_post) + call $subroutine_diexc(psi_generators(1,1,imask), & + generators_bitmask(1,1,d_hole1,i_bitmask_gen), & + generators_bitmask(1,1,d_part1,i_bitmask_gen), & + generators_bitmask(1,1,d_hole2,i_bitmask_gen), & + generators_bitmask(1,1,d_part2,i_bitmask_gen) & + $params_post) + enddo + + $copy_buffer + $generate_psi_guess + SOFT_TOUCH psi_det psi_coef + +end + diff --git a/src/Dets/README.rst b/src/Dets/README.rst index aaf90cfe..7c7ffade 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -53,6 +53,9 @@ Documentation `copy_h_apply_buffer_to_wf `_ Undocumented +`fill_h_apply_buffer_no_selection `_ + Fill the H_apply buffer with determiants for CISD + `h_apply_buffer_allocated `_ Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index a7a86f43..0a6ead21 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -1,9 +1,9 @@ determinants N_int integer bit_kind integer - n_dets integer + n_det integer n_states integer - psi_coef double precision (determinants_n_dets,determinants_n_states) - psi_det integer (determinants_N_int*determinants_bit_kind/4,2,determinants_n_dets) + psi_coef double precision (determinants_n_det,determinants_n_states) + psi_det integer*8 (determinants_N_int*determinants_bit_kind/8,2,determinants_n_det) H_apply_threshold double precision diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index a73a29ae..216b4f06 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -5,7 +5,16 @@ BEGIN_PROVIDER [ integer, N_states ] BEGIN_DOC ! Number of states to consider END_DOC - N_states = 1 + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_determinants_n_states(exists) + if (exists) then + call ezfio_get_determinants_n_states(N_states) + else + N_states = 1 + call ezfio_set_determinants_n_states(N_states) + endif + ASSERT (N_states > 0) END_PROVIDER BEGIN_PROVIDER [ integer, N_det ] @@ -13,22 +22,33 @@ BEGIN_PROVIDER [ integer, N_det ] BEGIN_DOC ! Number of determinants in the wave function END_DOC - N_det = 1 + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_determinants_n_det(exists) + if (exists) then + call ezfio_get_determinants_n_det(N_det) + else + N_det = 1 + call ezfio_set_determinants_n_det(N_det) + endif + ASSERT (N_det > 0) END_PROVIDER + BEGIN_PROVIDER [ integer, psi_det_size ] implicit none BEGIN_DOC ! Size of the psi_det/psi_coef arrays END_DOC - psi_det_size = 1000 + psi_det_size = 1000*N_states END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] implicit none BEGIN_DOC - ! The wave function. Initialized with Hartree-Fock + ! The wave function. Initialized with Hartree-Fock if the EZFIO file + ! is empty END_DOC integer, save :: ifirst = 0 @@ -52,3 +72,12 @@ END_PROVIDER END_PROVIDER + +BEGIN_PROVIDER [ integer, N_det_reference ] + implicit none + BEGIN_DOC + ! Number of determinants in the reference wave function + END_DOC + N_det_reference = N_det + ASSERT (N_det_reference > 0) +END_PROVIDER diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f new file mode 100644 index 00000000..7e0eef49 --- /dev/null +++ b/src/Dets/diagonalize_CI.irp.f @@ -0,0 +1,64 @@ +BEGIN_PROVIDER [ character*(64), diag_algorithm ] + implicit none + BEGIN_DOC + ! Diagonalization algorithm (Davidson or Lapack) + END_DOC + if (N_det > 500) then + diag_algorithm = "Davidson" + else + diag_algorithm = "Lapack" + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + do j=1,min(N_states,N_det) + CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ] + implicit none + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + integer :: i,j + + do j=1,N_states + do i=1,N_det + CI_eigenvectors(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, & + size(CI_eigenvectors,1),N_det,N_states,N_int,output_Dets) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + CI_electronic_energy(:) = 0.d0 + do j=1,min(N_states,N_det) + do i=1,N_det + CI_eigenvectors(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy(j) = eigenvalues(j) + enddo + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + diff --git a/src/Dets/utils.irp.f b/src/Dets/utils.irp.f index b16d44c9..0b89f339 100644 --- a/src/Dets/utils.irp.f +++ b/src/Dets/utils.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(n_det,n_det) ] +BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] implicit none BEGIN_DOC ! H matrix on the basis of the slater deter;inants defined by psi_det @@ -13,3 +13,4 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(n_det,n_det) ] enddo enddo END_PROVIDER + diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 8a5e66dc..df2d272a 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,8 +82,15 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd `_ - Undocumented +`pt2_moller_plesset `_ + compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states. + .br + c_pert(i) = /(difference of orbital energies) + .br + e_2_pert(i) = ^2/(difference of orbital energies) + .br `pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution @@ -95,7 +102,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various n_st states. @@ -105,7 +112,7 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`pt2_epstein_nesbet_2x2_sc2 `_ +`pt2_epstein_nesbet_2x2_sc2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various n_st states. @@ -127,7 +134,7 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`pt2_epstein_nesbet_sc2 `_ +`pt2_epstein_nesbet_sc2 `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various n_st states, @@ -149,7 +156,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ( ) ) .br -`pt2_epstein_nesbet_sc2_projected `_ +`pt2_epstein_nesbet_sc2_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various n_st states, @@ -175,6 +182,9 @@ Documentation NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !! .br +`cisd `_ + Undocumented + `fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection @@ -190,14 +200,6 @@ Documentation `diagonalize `_ Undocumented -`h_apply_cisd_pt2 `_ - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). - -`h_apply_cisd_selection `_ - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). - `n_det_ref `_ Undocumented diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index ac3ed0c4..35cd12d1 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -24,6 +24,7 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s ASSERT (Nint > 0) call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array) H_pert_diag = diag_H_mat_elem(det_pert,Nint) + print *, H_pert_diag do i =1,n_st c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag + eps) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) @@ -72,7 +73,7 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet integer(bit_kind), intent(in) :: det_pert(Nint,2) double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(ndet) + integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution @@ -122,7 +123,7 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint, integer(bit_kind), intent(in) :: det_pert(Nint,2) double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(ndet) + integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution @@ -174,7 +175,7 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag integer(bit_kind), intent(in) :: det_pert(Nint,2) double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag double precision :: i_H_psi_array(N_st) - integer :: idx_repeat(ndet) + integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 735c58c5..ebcb3048 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -24,7 +24,7 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buf ASSERT (N_st > 0) do i = 1,buffer_size - c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_generators,N_det,h_apply_threshold) + c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_generators,N_det_reference,h_apply_threshold) if (c_ref /= 0) then cycle diff --git a/src/Perturbation/perturbation_test.irp.f b/src/Perturbation/perturbation_test.irp.f index 9862df14..8f8aadfc 100644 --- a/src/Perturbation/perturbation_test.irp.f +++ b/src/Perturbation/perturbation_test.irp.f @@ -2,11 +2,6 @@ program cisd implicit none integer :: i,k double precision, allocatable :: eigvalues(:),eigvectors(:,:) - PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map - - -! N_states = 8 -! TOUCH N_states double precision, allocatable :: pt2(:), norm_pert(:) double precision :: H_pert_diag integer :: N_st @@ -14,7 +9,6 @@ program cisd allocate (pt2(N_st), norm_pert(N_st)) call H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st) - allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'pt2 = ', pt2 diff --git a/src/Perturbation/temporary_stuff.irp.f b/src/Perturbation/temporary_stuff.irp.f index e628d08f..f6734082 100644 --- a/src/Perturbation/temporary_stuff.irp.f +++ b/src/Perturbation/temporary_stuff.irp.f @@ -60,75 +60,3 @@ END_PROVIDER END_PROVIDER -subroutine H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st) - implicit none - BEGIN_DOC - ! Calls H_apply on the HF determinant and selects all connected single and double - ! excitations (of the same symmetry). - END_DOC - - integer, intent(in) :: N_st - double precision, intent(inout):: pt2(N_st) - double precision, intent(inout):: norm_pert(N_st) - double precision, intent(inout):: H_pert_diag - integer(bit_kind) :: hole_mask(N_int,2) - integer(bit_kind) :: particle_mask(N_int,2) - hole_mask(:,1) = HF_bitmask(:,1) - hole_mask(:,2) = HF_bitmask(:,2) - particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1)) - particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2)) - - - PROVIDE reference_energy N_det_generators key_pattern_not_in_ref - PROVIDE mo_bielec_integrals_in_map - PROVIDE H_apply_buffer_allocated - pt2 = 0.d0 - call H_apply_cisd_pt2_OpenMP_monoexc(HF_bitmask, & - hole_mask, particle_mask, & - pt2,norm_pert,H_pert_diag,N_st,N_int ) - call H_apply_cisd_pt2_OpenMP_diexc(HF_bitmask, & - hole_mask, particle_mask, & - hole_mask, particle_mask, & - pt2,norm_pert,H_pert_diag,N_st,N_int ) - -end - - -subroutine H_apply_cisd_selection(pt2, norm_pert, H_pert_diag, N_st) - implicit none - BEGIN_DOC - ! Calls H_apply on the HF determinant and selects all connected single and double - ! excitations (of the same symmetry). - END_DOC - - integer, intent(in) :: N_st - double precision, intent(inout):: pt2(N_st) - double precision, intent(inout):: norm_pert(N_st) - double precision, intent(inout):: H_pert_diag - integer(bit_kind) :: hole_mask(N_int,2) - integer(bit_kind) :: particle_mask(N_int,2) - hole_mask(:,1) = HF_bitmask(:,1) - hole_mask(:,2) = HF_bitmask(:,2) - particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1)) - particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2)) - - - PROVIDE reference_energy N_det_generators key_pattern_not_in_ref - PROVIDE mo_bielec_integrals_in_map selection_criterion - PROVIDE H_apply_buffer_allocated - pt2 = 0.d0 - norm_pert = 0.d0 - H_pert_diag = 0.d0 - call H_apply_cisd_selection_OpenMP_monoexc(HF_bitmask, & - hole_mask, particle_mask, & - pt2,norm_pert,H_pert_diag,N_st,N_int ) - call H_apply_cisd_selection_OpenMP_diexc(HF_bitmask, & - hole_mask, particle_mask, & - hole_mask, particle_mask, & - pt2,norm_pert,H_pert_diag,N_st,N_int ) - call copy_h_apply_buffer_to_wf - selection_criterion_min = selection_criterion_min*0.1d0 - selection_criterion = selection_criterion_min -end - - diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 22ca38a8..a5e12b31 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -196,31 +196,42 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) double precision, intent(in) :: H(nmax,n) double precision,allocatable :: eigenvalues(:) double precision,allocatable :: work(:) + integer ,allocatable :: iwork(:) double precision,allocatable :: A(:,:) - allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax)) - integer :: LWORK, info, i,j,l,k + integer :: lwork, info, i,j,l,k, liwork + + allocate(A(nmax,n),eigenvalues(n)) + A=H + lwork = 2*n*n + 6*n+ 1 + liwork = 5*n + 3 + allocate (work(lwork),iwork(liwork)) -! if (n<30) then -! do i=1,n -! do j=1,n -! print *, j,i, H(j,i) -! enddo -! print *, '---' -! enddo -! print *, '---' -! endif - - LWORK = 4*nmax - call dsyev( 'V', 'U', n, A, nmax, eigenvalues, work, LWORK, info ) + lwork = -1 + liwork = -1 + call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + iwork, liwork, info ) if (info < 0) then - print *, irp_here, ': the ',-info,'-th argument had an illegal value' - stop 1 - else if (info > 0) then - print *, irp_here, ': the algorithm failed to converge; ',info,' off-diagonal' - print *, 'elements of an intermediate tridiagonal form did not converge to zero.' - stop 1 + print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value' + stop 2 endif + lwork = int( work( 1 ) ) + liwork = iwork(1) + deallocate (work,iwork) + + allocate (work(lwork),iwork(liwork)) + call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + iwork, liwork, info ) + deallocate(work,iwork) + + if (info < 0) then + print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'DSYEVD Failed' + stop 1 + end if + eigvectors = 0.d0 eigvalues = 0.d0 do j = 1, n @@ -229,4 +240,5 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) eigvectors(i,j) = A(i,j) enddo enddo + deallocate(A,eigenvalues) end