From 80f863ef3ae455184d12a6e3062ec5738046a11b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 17 May 2014 14:20:55 +0200 Subject: [PATCH] Accelerated Davidson --- scripts/generate_h_apply.py | 63 ++++++- scripts/perturbation.py | 26 +++ src/CISD/H_apply.irp.f | 10 +- src/CISD/cisd.irp.f | 3 + src/CISD/h_apply.py | 12 +- src/Dets/H_apply_template.f | 18 +- src/Dets/README.rst | 18 +- src/Dets/connected_to_ref.irp.f | 254 ++++++++++++++++++++++++++++ src/Dets/davidson.irp.f | 12 +- src/Dets/filter_connected.irp.f | 174 +++++++++++++++++++ src/Dets/slater_rules.irp.f | 211 +++-------------------- src/NEEDED_MODULES | 2 +- src/Perturbation/README.rst | 105 ++++++++++++ src/Perturbation/perturbation.irp.f | 13 ++ 14 files changed, 683 insertions(+), 238 deletions(-) create mode 100755 scripts/perturbation.py create mode 100644 src/Dets/connected_to_ref.irp.f create mode 100644 src/Dets/filter_connected.irp.f create mode 100644 src/Perturbation/README.rst create mode 100644 src/Perturbation/perturbation.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 5fa3e436..74c35116 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -14,10 +14,17 @@ keys_work finalization """.split() -def new_dict(openmp=True): - s ={} +class H_apply(object): + + def __init__(self,sub,openmp=True): + s = {} for k in keywords: s[k] = "" + s["subroutine"] = "H_apply_%s"%(sub) + self.openmp = openmp + if openmp: + s["subroutine"] += "_OpenMP" + self.perturbation = None #s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) & s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, & @@ -29,7 +36,7 @@ def new_dict(openmp=True): !$OMP N_elec_in_key_hole_2,ia_ja_pairs) & !$OMP SHARED(key_in,N_int,elec_num_tab, & !$OMP hole_1, particl_1, hole_2, particl_2, & - !$OMP lck,thresh,elec_alpha_num,E_ref)""" + !$OMP lck,thresh,elec_alpha_num)""" s["omp_init_lock"] = "call omp_init_lock(lck)" s["omp_set_lock"] = "call omp_set_lock(lck)" s["omp_unset_lock"] = "call omp_unset_lock(lck)" @@ -50,16 +57,54 @@ def new_dict(openmp=True): s["set_i_H_j_threshold"] = """ thresh = H_apply_threshold """ - return s + self.data = s + def __setitem__(self,key,value): + self.data[key] = value + def __getitem__(self,key): + return self.data[key] -def create_h_apply(s): + def __repr__(self): buffer = template - for key in s: - buffer = buffer.replace('$'+key, s[key]) - print buffer - + for key,value in self.data.items(): + buffer = buffer.replace('$'+key, value) + return buffer + + def set_perturbation(self,pert): + self.perturbation = pert + if pert is not None: + self.data["parameters"] = ",sum_e_2_pert_in,sum_norm_pert_in,sum_H_pert_diag_in,N_st,Nint" + self.data["declarations"] = """ + integer, intent(in) :: N_st,Nint + double precision, intent(inout) :: sum_e_2_pert_in(N_st) + double precision, intent(inout) :: sum_norm_pert_in(N_st) + double precision, intent(inout) :: sum_H_pert_diag_in + double precision :: sum_e_2_pert(N_st) + double precision :: sum_norm_pert(N_st) + double precision :: sum_H_pert_diag + """ + self.data["size_max"] = "256" + self.data["initialization"] = """ + E_ref = diag_H_mat_elem(key_in,N_int) + sum_e_2_pert = sum_e_2_pert_in + sum_norm_pert = sum_norm_pert_in + sum_H_pert_diag = sum_H_pert_diag_in + """ + self.data["keys_work"] += """ + call perturb_buffer_%s(keys_out,key_idx,sum_e_2_pert, & + sum_norm_pert,sum_H_pert_diag,N_st,Nint) + """%(pert,) + self.data["finalization"] = """ + sum_e_2_pert_in = sum_e_2_pert + sum_norm_pert_in = sum_norm_pert + sum_H_pert_diag_in = sum_H_pert_diag + """ + if self.openmp: + self.data["omp_test_lock"] = ".False." + self.data["omp_parallel"] += """& + !$OMP SHARED(N_st) & + !$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)""" diff --git a/scripts/perturbation.py b/scripts/perturbation.py new file mode 100755 index 00000000..8b0f8b01 --- /dev/null +++ b/scripts/perturbation.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +import os + +Pert_dir = os.environ["QPACKAGE_ROOT"]+"/src/Perturbation/" + +perturbations = [] + +for filename in filter(lambda x: x.endswith(".irp.f"), os.listdir(Pert_dir)): + + filename = Pert_dir+filename + file = open(filename,'r') + lines = file.readlines() + file.close() + for line in lines: + buffer = line.lower().lstrip().split() + if len(buffer) > 1: + if buffer[0] == "subroutine" and buffer[1].startswith("pt2_"): + p = (buffer[1].split('(')[0])[4:] + perturbations.append( p ) + + +if __name__ == '__main__': + print 'Perturbations:' + for k in perturbations: + print '* ', k diff --git a/src/CISD/H_apply.irp.f b/src/CISD/H_apply.irp.f index b7f27dfc..be80e6dd 100644 --- a/src/CISD/H_apply.irp.f +++ b/src/CISD/H_apply.irp.f @@ -54,14 +54,14 @@ subroutine H_apply_cisd particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1)) particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2)) - call H_apply_cisd_monoexc(HF_bitmask, & + + call H_apply_cisd_OpenMP_monoexc(HF_bitmask, & hole_mask, particle_mask) - call H_apply_cisd_diexc(HF_bitmask, & + call H_apply_cisd_OpenMP_diexc(HF_bitmask, & hole_mask, particle_mask, & hole_mask, particle_mask ) - - call copy_H_apply_buffer_to_wf - + + call copy_h_apply_buffer_to_wf end diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index d2f6de94..ba84b8c6 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -3,6 +3,9 @@ program cisd integer :: i,k double precision, allocatable :: eigvalues(:),eigvectors(:,:) PROVIDE ref_bitmask_energy + + double precision :: pt2(10), norm_pert(10), H_pert_diag + call H_apply_cisd allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det diff --git a/src/CISD/h_apply.py b/src/CISD/h_apply.py index 25a1e328..138415d0 100755 --- a/src/CISD/h_apply.py +++ b/src/CISD/h_apply.py @@ -1,11 +1,11 @@ #!/usr/bin/env python -import generate_h_apply +from generate_h_apply import * -# H_apply -s = generate_h_apply.new_dict(openmp=True) -s["subroutine"] = "H_apply_cisd" -s["keys_work"] = "call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int)" -generate_h_apply.create_h_apply(s) +s = H_apply("cisd",openmp=True) +s["keys_work"] += """ +call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int) +""" +print s diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 3c13d31f..6a7efa8d 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -2,6 +2,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame use omp_lib use bitmasks implicit none + BEGIN_DOC + ! Generate all double excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -24,7 +29,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame integer,parameter :: size_max = $size_max double precision :: hij_elec, mo_bielec_integral, thresh integer, allocatable :: ia_ja_pairs(:,:,:) - double precision :: diag_H_mat_elem, E_ref + double precision :: diag_H_mat_elem PROVIDE mo_integrals_map ref_bitmask_energy PROVIDE mo_bielec_integrals_in_map @@ -34,8 +39,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame $omp_init_lock - E_ref = diag_H_mat_elem(key_in,N_int) - $initialization $omp_parallel @@ -233,6 +236,11 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) use omp_lib use bitmasks implicit none + BEGIN_DOC + ! Generate all single excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -255,7 +263,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer,parameter :: size_max = $size_max double precision :: hij_elec, thresh integer, allocatable :: ia_ja_pairs(:,:,:) - double precision :: diag_H_mat_elem, E_ref + double precision :: diag_H_mat_elem PROVIDE mo_integrals_map ref_bitmask_energy PROVIDE mo_bielec_integrals_in_map @@ -265,8 +273,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) $omp_init_lock - E_ref = diag_H_mat_elem(key_in,N_int) - $initialization $omp_parallel diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 69c9cd33..ab562b17 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -135,10 +135,10 @@ Documentation `get_s2 `_ Returns -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem `decode_exc `_ @@ -148,15 +148,9 @@ Documentation s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`filter_connected `_ - Filters out the determinants that are not connected by H - -`filter_connected_i_h_psi0 `_ - Undocumented - `get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase @@ -172,10 +166,10 @@ Documentation `get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -185,7 +179,7 @@ Documentation `i_h_j `_ Returns where i and j are determinants -`i_h_psim `_ +`i_h_psi `_ Undocumented `h_matrix_all_dets `_ diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f new file mode 100644 index 00000000..81e4e12a --- /dev/null +++ b/src/Dets/connected_to_ref.irp.f @@ -0,0 +1,254 @@ +integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh) + use bitmasks + implicit none + integer, intent(in) :: Nint, N_past_in, Ndet + integer(bit_kind), intent(in) :: keys(ishft(Nint,-1),2,Ndet) + integer(bit_kind), intent(in) :: key(ishft(Nint,-1),2) + double precision, intent(in) :: thresh + + integer :: N_past + integer :: i, l + integer :: degree_x2 + logical :: det_is_not_or_may_be_in_ref, t + double precision :: hij_elec + + ! output : 0 : not connected + ! i : connected to determinant i of the past + ! -i : is the ith determinant of the refernce wf keys + + ASSERT (Nint == N_int) + + connected_to_ref = 0 + N_past = max(1,N_past_in) + if (Nint == 1) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + do i=N_past,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)) ) then + cycle + endif + connected_to_ref = -i + return + enddo + return + + else if (Nint==2) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + !DIR$ LOOP COUNT (1000) + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)).or. & + (key(2,1) /= keys(2,1,i)).or. & + (key(2,2) /= keys(2,2,i)) ) then + cycle + endif + connected_to_ref = -i + return + enddo + return + + else if (Nint==3) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + & + popcnt(xor( key(3,1), keys(3,1,i))) + & + popcnt(xor( key(3,2), keys(3,2,i))) + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)).or. & + (key(2,1) /= keys(2,1,i)).or. & + (key(2,2) /= keys(2,2,i)).or. & + (key(3,1) /= keys(3,1,i)).or. & + (key(3,2) /= keys(3,2,i)) ) then + cycle + endif + connected_to_ref = -i + return + enddo + return + + else + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) + & + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if(degree_x2 == 0)then + connected_to_ref = -i + return + endif + if (degree_x2 > 5) then + cycle + else + call i_H_j(keys(1,1,i),key,N_int,hij_elec) + if(dabs(hij_elec).lt.thresh)cycle + connected_to_ref = i + return + endif + enddo + + !DIR$ FORCEINLINE + t = det_is_not_or_may_be_in_ref(key,N_int) + if ( t ) then + return + endif + + do i=N_past+1,Ndet + if ( (key(1,1) /= keys(1,1,i)).or. & + (key(1,2) /= keys(1,2,i)) ) then + cycle + else + connected_to_ref = -1 + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + if ( (key(l,1) /= keys(l,1,i)).or. & + (key(l,2) /= keys(l,2,i)) ) then + connected_to_ref = 0 + exit + endif + enddo + if (connected_to_ref /= 0) then + return + endif + endif + enddo + + endif + +end + +logical function det_is_not_or_may_be_in_ref(key,Nint) + implicit none + ! If true, det is not in ref + ! If false, det may be in ref + + integer, intent(in) :: key(Nint,2), Nint + integer :: key_int + integer*1 :: key_short(4) + !DIR$ ATTRIBUTES ALIGN : 32 :: key_short + equivalence (key_int,key_short) + + integer :: i, ispin + + det_is_not_or_may_be_in_ref = .False. + do ispin=1,2 + do i=1,Nint + key_int = key(i,ispin) + if ( & + key_pattern_not_in_ref(key_short(1), i, ispin) & + .or.key_pattern_not_in_ref(key_short(2), i, ispin) & + .or.key_pattern_not_in_ref(key_short(3), i, ispin) & + .or.key_pattern_not_in_ref(key_short(4), i, ispin) & + ) then + det_is_not_or_may_be_in_ref = .True. + return + endif + enddo + enddo + +end + + +BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ] + implicit none + BEGIN_DOC +! Min and max values of the integers of the keys of the reference + END_DOC + + integer :: i, j, ispin + integer :: key + integer*1 :: key_short(4) + equivalence (key,key_short) + integer :: idx + + key_pattern_not_in_ref = .True. + + do j=1,N_det + do ispin=1,2 + do i=1,N_int + key = psi_det(i,ispin,j) + key_pattern_not_in_ref( key_short(1), i, ispin ) = .False. + key_pattern_not_in_ref( key_short(2), i, ispin ) = .False. + key_pattern_not_in_ref( key_short(3), i, ispin ) = .False. + key_pattern_not_in_ref( key_short(4), i, ispin ) = .False. + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index 780d54b5..f0dc2437 100644 --- a/src/Dets/davidson.irp.f +++ b/src/Dets/davidson.irp.f @@ -133,7 +133,6 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) !$OMP END PARALLEL do iter=1,davidson_sze_max-1 - print *, 'iter = ',iter ! print *, '***************' ! do i=1,iter @@ -174,7 +173,6 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ! ------------- call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) - print *, lambda(1:4) ! Express eigenvectors of h in the determinant basis ! -------------------------------------------------- @@ -203,13 +201,11 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) enddo residual_norm(k) = u_dot_u(R(1,k),sze) enddo - print *, 'Lambda' - print *, lambda(1:N_st) + nuclear_repulsion - print *, 'Residual_norm' - print *, residual_norm(1:N_st) - print *, '' + + print '(I3,15(F16.8,x))', iter, lambda(1:N_st) + nuclear_repulsion + print '(3x,15(E16.5,x))', residual_norm(1:N_st) - converged = maxval(residual_norm) < 1.d-5 + converged = maxval(residual_norm) < 1.d-10 if (converged) then exit endif diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f new file mode 100644 index 00000000..18acd748 --- /dev/null +++ b/src/Dets/filter_connected.irp.f @@ -0,0 +1,174 @@ + +subroutine filter_connected(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC +! Filters out the determinants that are not connected by H + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + +subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,l + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (sze > 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + endif + enddo + + endif + idx(0) = l-1 +end + diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 59dc68ad..c160b0c2 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -488,7 +488,7 @@ end -subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) +subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate @@ -600,180 +600,6 @@ end -subroutine filter_connected(key1,key2,Nint,sze,idx) - use bitmasks - implicit none - BEGIN_DOC -! Filters out the determinants that are not connected by H - END_DOC - integer, intent(in) :: Nint, sze - integer(bit_kind), intent(in) :: key1(Nint,2,sze) - integer(bit_kind), intent(in) :: key2(Nint,2) - integer, intent(out) :: idx(0:sze) - - integer :: i,j,l - integer :: degree_x2 - - ASSERT (Nint > 0) - ASSERT (sze > 0) - - l=1 - - if (Nint==1) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) - if (degree_x2 < 5) then - idx(l) = i - l = l+1 - endif - enddo - - else if (Nint==2) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) + & - popcnt(xor( key1(2,1,i), key2(2,1))) + & - popcnt(xor( key1(2,2,i), key2(2,2))) - if (degree_x2 < 5) then - idx(l) = i - l = l+1 - endif - enddo - - else if (Nint==3) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) + & - popcnt(xor( key1(2,1,i), key2(2,1))) + & - popcnt(xor( key1(2,2,i), key2(2,2))) + & - popcnt(xor( key1(3,1,i), key2(3,1))) + & - popcnt(xor( key1(3,2,i), key2(3,2))) - if (degree_x2 < 5) then - idx(l) = i - l = l+1 - endif - enddo - - else - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = 0 - !DEC$ LOOP COUNT MIN(4) - do j=1,Nint - degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& - popcnt(xor( key1(j,2,i), key2(j,2))) - if (degree_x2 > 4) then - exit - endif - enddo - if (degree_x2 <= 5) then - idx(l) = i - l = l+1 - endif - enddo - - endif - idx(0) = l-1 -end - -subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) - use bitmasks - implicit none - integer, intent(in) :: Nint, sze - integer(bit_kind), intent(in) :: key1(Nint,2,sze) - integer(bit_kind), intent(in) :: key2(Nint,2) - integer, intent(out) :: idx(0:sze) - - integer :: i,l - integer :: degree_x2 - - ASSERT (Nint > 0) - ASSERT (sze > 0) - - l=1 - - if (Nint==1) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) - if (degree_x2 < 5) then - if(degree_x2 .ne. 0)then - idx(l) = i - l = l+1 - endif - endif - enddo - - else if (Nint==2) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) + & - popcnt(xor( key1(2,1,i), key2(2,1))) + & - popcnt(xor( key1(2,2,i), key2(2,2))) - if (degree_x2 < 5) then - if(degree_x2 .ne. 0)then - idx(l) = i - l = l+1 - endif - endif - enddo - - else if (Nint==3) then - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & - popcnt(xor( key1(1,2,i), key2(1,2))) + & - popcnt(xor( key1(2,1,i), key2(2,1))) + & - popcnt(xor( key1(2,2,i), key2(2,2))) + & - popcnt(xor( key1(3,1,i), key2(3,1))) + & - popcnt(xor( key1(3,2,i), key2(3,2))) - if (degree_x2 < 5) then - if(degree_x2 .ne. 0)then - idx(l) = i - l = l+1 - endif - endif - enddo - - else - - !DIR$ LOOP COUNT (1000) - do i=1,sze - degree_x2 = 0 - !DEC$ LOOP COUNT MIN(4) - do l=1,Nint - degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& - popcnt(xor( key1(l,2,i), key2(l,2))) - if (degree_x2 > 4) then - exit - endif - enddo - if (degree_x2 <= 5) then - if(degree_x2 .ne. 0)then - idx(l) = i - l = l+1 - endif - endif - enddo - - endif - idx(0) = l-1 -end - - double precision function diag_H_mat_elem(det_in,Nint) implicit none @@ -963,6 +789,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) integer, allocatable :: idx(:) double precision :: hij + double precision, allocatable :: vt(:) integer :: i,j,k,l, jj integer :: i0, j0 ASSERT (Nint > 0) @@ -971,32 +798,34 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) PROVIDE ref_bitmask_energy integer, parameter :: block_size = 157 !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,idx,jj) SHARED(n,H_jj,u_0,keys_tmp,Nint)& + !$OMP PRIVATE(i,hij,j,k,idx,jj,vt) SHARED(n,H_jj,u_0,keys_tmp,Nint)& !$OMP SHARED(v_0) - allocate(idx(0:block_size)) !$OMP DO SCHEDULE(static) do i=1,n v_0(i) = H_jj(i) * u_0(i) enddo !$OMP END DO + allocate(idx(0:n), vt(n)) + Vt = 0.d0 !$OMP DO SCHEDULE(guided) - do i0=1,n,block_size - do j0=1,n,block_size - do i=i0,min(i0+block_size-1,n) - call filter_connected(keys_tmp(1,1,j0),keys_tmp(1,1,i),Nint,min(block_size,i-j0+1),idx) + do i=1,n + call filter_connected(keys_tmp(1,1,1),keys_tmp(1,1,i),Nint,i-1,idx) do jj=1,idx(0) - j = idx(jj)+j0-1 - if ( (j 1.d-8)) then - call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) - v_0(i) = v_0(i) + hij*u_0(j) - v_0(j) = v_0(j) + hij*u_0(i) - endif + j = idx(jj) + call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) + vt (i) = vt (i) + hij*u_0(j) + vt (j) = vt (j) + hij*u_0(i) enddo - enddo - enddo enddo - !$OMP END DO NOWAIT - deallocate(idx) + !$OMP END DO + !$OMP CRITICAL + do i=1,n + v_0(i) = v_0(i) + vt(i) + enddo + !$OMP END CRITICAL + deallocate(idx,vt) !$OMP END PARALLEL end + + diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 1a75cb06..d85c28b1 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD +AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst new file mode 100644 index 00000000..41a8d6a9 --- /dev/null +++ b/src/Perturbation/README.rst @@ -0,0 +1,105 @@ +=================== +Perturbation Module +=================== + + +All subroutines in `*.irp.f` starting with ``pt2_`` in the current directory are +perturbation computed using the routine ``i_H_psi``. Other cases are not allowed. +The arguments of the ``pt2_`` are always: + + subroutine pt2_...( & + psi_ref, & + psi_ref_coefs, & + E_refs, & + det_pert, & + c_pert, & + e_2_pert, & + H_pert_diag, & + Nint, & + ndet, & + n_st ) + + + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: psi_ref(Nint,2,ndet) + double precision , intent(in) :: psi_ref_coefs(ndet,n_st) + double precision , intent(in) :: E_refs(n_st) + 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 + + +psi_ref + bitstring of the determinants present in the various n_st states + +psi_ref_coefs + coefficients of the determinants on the various n_st states + +E_refs + Energy of the various n_st states + +det_pert + Perturber determinant + +c_pert + Pertrubative coefficients for the various states + +e_2_pert + Perturbative energetic contribution for the various states + +H_pert_diag + Diagonal H matrix element of the perturber + +Nint + Should be equal to N_int + +Ndet + Number of determinants `i` in Psi on which we apply + +N_st + Number of states + + + + + +Assumptions +=========== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* This is not allowed: + + subroutine & + pt2_.... + + + + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ + diff --git a/src/Perturbation/perturbation.irp.f b/src/Perturbation/perturbation.irp.f new file mode 100644 index 00000000..8ce91f37 --- /dev/null +++ b/src/Perturbation/perturbation.irp.f @@ -0,0 +1,13 @@ +BEGIN_SHELL [ /usr/bin/env python ] +from perturbation import perturbations +import os + +filename = os.environ["QPACKAGE_ROOT"]+"/src/Perturbation/perturbation_template.f" +file = open(filename,'r') +template = file.read() +file.close() + +for p in perturbations: + print template.replace("$PERT",p) + +END_SHELL