diff --git a/scripts/check_dependencies.sh b/scripts/check_dependencies.sh index 943c5211..4e7e4bf1 100755 --- a/scripts/check_dependencies.sh +++ b/scripts/check_dependencies.sh @@ -51,7 +51,7 @@ DEPS=$(unique_list $DEPS_LONG) if [[ ! "$COMMAND_LINE" == "$DEPS" ]] then - DEPS=$(check_dependencies.sh $DEPS) + DEPS=$(${QPACKAGE_ROOT}/scripts/check_dependencies.sh ${DEPS}) fi echo $DEPS diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 5fa3e436..f40833d7 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -10,14 +10,25 @@ subroutine parameters initialization declarations -keys_work +keys_work_locked +keys_work_unlocked 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.selection_pt2 = None + 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, & @@ -27,9 +38,9 @@ def new_dict(openmp=True): !$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,& !$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, & !$OMP N_elec_in_key_hole_2,ia_ja_pairs) & - !$OMP SHARED(key_in,N_int,elec_num_tab, & + !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & !$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 +61,76 @@ 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): + if self.perturbation is not None: + raise + 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 + double precision :: e_2_pert_buffer(N_st,size_max) + double precision :: coef_pert_buffer(N_st,size_max) + """ + self.data["size_max"] = "256" + self.data["initialization"] = """ + 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_unlocked"] += """ + call perturb_buffer_%s(keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,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_set_lock"] = "" + self.data["omp_unset_lock"] = "" + self.data["omp_test_lock"] = ".False." + self.data["omp_parallel"] += """& + !$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) & + !$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)""" + + def set_selection_pt2(self,pert): + if self.selection_pt2 is not None: + raise + self.set_perturbation(pert) + self.selection_pt2 = pert + if pert is not None: + self.data["size_max"] = str(1024*128) + self.data["keys_work_unlocked"] = """ + e_2_pert_buffer = 0.d0 + coef_pert_buffer = 0.d0 + """ + self.data["keys_work_unlocked"] + self.data["keys_work_locked"] = """ + call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer,coef_pert_buffer,N_st,N_int) + """ + self.data["omp_test_lock"] = "omp_test_lock(lck)" + self.data["omp_set_lock"] = "call omp_set_lock(lck)" + self.data["omp_unset_lock"] = "call omp_unset_lock(lck)" 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/README.rst b/src/CISD/README.rst index 1a9c86d9..fbf1807f 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -23,12 +23,15 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fill_h_apply_buffer_cisd `_ +`fill_h_apply_buffer_cisd `_ Fill the H_apply buffer with determiants for CISD -`h_apply_cisd `_ +`h_apply_cisd `_ Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). +`cisd `_ + Undocumented + diff --git a/src/CISD/SC2.irp.f b/src/CISD/SC2.irp.f new file mode 100644 index 00000000..330e3665 --- /dev/null +++ b/src/CISD/SC2.irp.f @@ -0,0 +1,421 @@ +subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + PROVIDE ref_bitmask_energy + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + integer :: iter + integer :: i,j,k,l,m + logical :: converged + double precision :: overlap(N_st,N_st) + double precision :: u_dot_v, u_dot_u + + integer :: degree,N_double,index_hf,index_double(sze) + double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0 + double precision :: e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),hij_double(sze) + double precision :: e_corr_double_before,accu,cpu_2,cpu_1 + integer :: i_ok + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(sze,N_st, & + !$OMP H_jj_ref,Nint,dets_in,u_in) & + !$OMP PRIVATE(i) + + !$OMP DO + do i=1,sze + H_jj_ref(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + N_double = 0 + e_corr = 0.d0 + e_corr_double = 0.d0 + do i = 1, sze + call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint) + if(degree==0)then + index_hf=i + else if (degree == 2)then + N_double += 1 + index_double(N_double) = i + call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) + hij_double(N_double) = hij_elec + e_corr_array(N_double) = u_in(i,1)* hij_elec + e_corr_double += e_corr_array(N_double) + e_corr += e_corr_array(N_double) + index_double(N_double) = i + else if (degree == 1)then + call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) + print*,hij_elec + e_corr += u_in(i,1)* hij_elec + endif + enddo + inv_c0 = 1.d0/u_in(index_hf,1) + do i = 1, N_double + e_corr_array(i) = e_corr_array(i) * inv_c0 + enddo + e_corr = e_corr * inv_c0 + e_corr_double = e_corr_double * inv_c0 + print*, 'E_corr = ',e_corr + print*, 'E_corr_double = ', e_corr_double + + converged = .False. + e_corr_double_before = e_corr_double + iter = 0 + do while (.not.converged) + + iter +=1 + print*,'SC2 iteration : ',iter + call cpu_time(cpu_1) + do i=1,sze + H_jj_dressed(i) = H_jj_ref(i) + if (i==index_hf)cycle + accu = 0.d0 + do j=1,N_double + call repeat_excitation(dets_in(1,1,i),ref_bitmask,dets_in(1,1,index_double(j)),i_ok,Nint) + if (i_ok==1)cycle! you check if the excitation is possible + accu += e_corr_array(j) + enddo + H_jj_dressed(i) += accu + enddo + + call cpu_time(cpu_2) + print*,'time for the excitations = ',cpu_2 - cpu_1 + print*,H_jj_ref(1),H_jj_ref(2) + print*,H_jj_dressed(1),H_jj_dressed(2) + print*,u_in(index_hf,1),u_in(index_double(1),1) + call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint) + print*,u_in(index_hf,1),u_in(index_double(1),1) + e_corr_double = 0.d0 + inv_c0 = 1.d0/u_in(index_hf,1) + do i = 1, N_double + e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i) + e_corr_double += e_corr_array(i) + enddo + print*,'E_corr = ',e_corr_double + print*,'delta E_corr =',e_corr_double - e_corr_double_before + converged = dabs(e_corr_double - e_corr_double_before) < 1.d-10 + if (converged) then + exit + endif + e_corr_double_before = e_corr_double + + enddo + + + +end + +subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Davidson diagonalization with specific diagonal elements of the H matrix + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(in) :: H_jj(dim_in) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + + integer :: iter + integer :: i,j,k,l,m + logical :: converged + + double precision :: overlap(N_st,N_st) + double precision :: u_dot_v, u_dot_u + + integer, allocatable :: kl_pairs(:,:) + integer :: k_pairs, kl + + integer :: iter2 + double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) + double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) + double precision :: diag_h_mat_elem + double precision :: residual_norm(N_st) + + PROVIDE ref_bitmask_energy + + allocate( & + kl_pairs(2,N_st*(N_st+1)/2), & + W(sze,N_st,davidson_sze_max), & + U(sze,N_st,davidson_sze_max), & + R(sze,N_st), & + h(N_st,davidson_sze_max,N_st,davidson_sze_max), & + y(N_st,davidson_sze_max,N_st,davidson_sze_max), & + lambda(N_st*davidson_sze_max)) + + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + ! Initialization + ! ============== + + k_pairs=0 + do l=1,N_st + do k=1,l + k_pairs+=1 + kl_pairs(1,k_pairs) = k + kl_pairs(2,k_pairs) = l + enddo + enddo + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & + !$OMP Nint,dets_in,u_in) & + !$OMP PRIVATE(k,l,kl,i) + + + ! Orthonormalize initial guess + ! ============================ + + !$OMP DO + do kl=1,k_pairs + k = kl_pairs(1,kl) + l = kl_pairs(2,kl) + if (k/=l) then + overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) + overlap(l,k) = overlap(k,l) + else + overlap(k,k) = u_dot_u(U_in(1,k),sze) + endif + enddo + !$OMP END DO + !$OMP END PARALLEL + + call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) + + ! Davidson iterations + ! =================== + + converged = .False. + + do while (.not.converged) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) + do k=1,N_st + !$OMP DO + do i=1,sze + U(i,k,1) = u_in(i,k) + enddo + !$OMP END DO + enddo + !$OMP END PARALLEL + + do iter=1,davidson_sze_max-1 + + ! Compute W_k = H |u_k> + ! ---------------------- + + do k=1,N_st + call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + enddo + + ! Compute h_kl = = + ! ------------------------------------------- + + do l=1,N_st + do k=1,N_st + do iter2=1,iter-1 + h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) + h(k,iter,l,iter2) = h(k,iter2,l,iter) + enddo + enddo + do k=1,l + h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) + h(l,iter,k,iter) = h(k,iter,l,iter) + enddo + enddo + + ! Diagonalize h + ! ------------- + call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + ! call dgemm ( 'N','N', sze, N_st*iter, N_st, & + ! 1.d0, U(1,1,1), size(U,1), y(1,1,1,1), size(y,1)*size(y,2), & + ! 0.d0, U(1,1,iter+1), size(U,1) ) + do k=1,N_st + do i=1,sze + U(i,k,iter+1) = 0.d0 + W(i,k,iter+1) = 0.d0 + do l=1,N_st + do iter2=1,iter + U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) + W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) + enddo + enddo + enddo + enddo + + ! Compute residual vector + ! ----------------------- + + do k=1,N_st + do i=1,sze + R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) + enddo + residual_norm(k) = u_dot_u(R(1,k),sze) + enddo + + 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-10 + if (converged) then + exit + endif + + ! Davidson step + ! ------------- + + do k=1,N_st + do i=1,sze + U(i,k,iter+1) = 1.d0/(lambda(k) - H_jj(i)) * R(i,k) + enddo + enddo + + ! Gram-Schmidt + ! ------------ + + double precision :: c + do k=1,N_st + do iter2=1,iter + do l=1,N_st + c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) + do i=1,sze + U(i,k,iter+1) -= c * U(i,l,iter2) + enddo + enddo + enddo + do l=1,k-1 + c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze) + do i=1,sze + U(i,k,iter+1) -= c * U(i,l,iter+1) + enddo + enddo + call normalize( U(1,k,iter+1), sze ) + enddo + enddo + + if (.not.converged) then + iter = davidson_sze_max-1 + endif + + ! Re-contract to u_in + ! ----------- + + do k=1,N_st + energies(k) = lambda(k) + do i=1,sze + u_in(i,k) = 0.d0 + do iter2=1,iter + do l=1,N_st + u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1) + enddo + enddo + enddo + enddo + + enddo + + deallocate ( & + kl_pairs, & + W, & + U, & + R, & + h, & + y, & + lambda & + ) +end + +subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint) + use bitmasks + implicit none + integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint + integer,intent(out):: i_ok + integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2) + double precision :: phase + i_ok = 1 + call get_excitation(key_1,key_2,exc,degree,phase,Nint) + integer :: h1,p1,h2,p2,s1,s2 + if(degree==2)then + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ! first hole + k_hole = ishft(h1-1,-5)+1 + j_hole = h1-ishft(k_hole-1,5)-1 + if(iand(key_in(k_hole,s1),ibset(0,j_hole)).eq.0)then + i_ok = 0 + return + endif + + ! second hole + k_hole = ishft(h2-1,-5)+1 + j_hole = h2-ishft(k_hole-1,5)-1 + if(iand(key_in(k_hole,s2),ibset(0,j_hole)).eq.0)then + i_ok = 0 + return + endif + + ! first particle + k_particl = ishft(p1-1,-5)+1 + j_particl = p1-ishft(k_particl-1,5)-1 + if(iand(key_in(k_particl,s1),ibset(0,j_particl)).ne.0)then + i_ok = 0 + return + endif + + ! second particle + k_particl = ishft(p2-1,-5)+1 + j_particl = p2-ishft(k_particl-1,5)-1 + if(iand(key_in(k_particl,s2),ibset(0,j_particl)).ne.0)then + i_ok = 0 + return + endif + return + endif +end + diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index d2f6de94..f3bb5488 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -3,6 +3,11 @@ program cisd integer :: i,k double precision, allocatable :: eigvalues(:),eigvectors(:,:) PROVIDE ref_bitmask_energy + + double precision :: pt2(10), norm_pert(10), H_pert_diag + + N_states = 3 + touch N_states call H_apply_cisd allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det @@ -17,7 +22,9 @@ program cisd print *, 'HF:', HF_energy print *, '---' do i = 1,1 - print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion + 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) deallocate(eigvalues,eigvectors) end 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/CIS_dressed/CIS_providers.irp.f b/src/CIS_dressed/CIS_providers.irp.f new file mode 100644 index 00000000..b05f87cb --- /dev/null +++ b/src/CIS_dressed/CIS_providers.irp.f @@ -0,0 +1,334 @@ + + use bitmasks + BEGIN_PROVIDER [integer(bit_kind), psi_CIS,(N_int,2,size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_holes,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_particl,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_spin,(size_psi_CIS)] + &BEGIN_PROVIDER [integer, psi_CIS_adress,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, H_CIS,(size_psi_CIS,size_psi_CIS)] + BEGIN_DOC + !key of the CIS-matrix + END_DOC + + + implicit none + integer :: a !control variable + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + integer :: key !key for CIS-matrix + integer :: i_hole,j_hole,ispin,l_particle,k_particle + double precision :: hij + + + do a=1,N_int + psi_CIS(a,1,1)=ref_bitmask(a,1) + psi_CIS(a,2,1)=ref_bitmask(a,2) + enddo + + psi_CIS_holes(1) = 0 + psi_CIS_particl(1) = 0 + psi_CIS_spin(1) = 0 + + !loop on particles: create a particle in k + do k=elec_alpha_num+1,n_act_cis + + !loop on holes: destroy a particle in i + do i=n_core_cis+1,elec_alpha_num + + !alpha spin + ispin=1 + + key=2*((k-elec_alpha_num-1)*(elec_alpha_num-n_core_cis) + i-n_core_cis) !index of such an excited determinant in the CIS WF + psi_CIS_adress(i,k)=key + + do a=1,N_int + psi_CIS(a,1,key)=ref_bitmask(a,1) + psi_CIS(a,2,key)=ref_bitmask(a,2) + enddo + + j_hole=ishft(i-1,-5)+1 + i_hole=i-ishft(j_hole-1,5)-1 + + psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole) + + l_particle=ishft(k-1,-5)+1 + k_particle=k-ishft(l_particle-1,5)-1 + + psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle) + + psi_CIS_holes(key) = i + psi_CIS_particl(key) = k + psi_CIS_spin(key) = 1 + + !beta spin + ispin=2 + + key=key+1 + + do a=1,N_int + psi_CIS(a,1,key)=ref_bitmask(a,1) + psi_CIS(a,2,key)=ref_bitmask(a,2) + enddo + + j_hole=ishft(i-1,-5)+1 + i_hole=i-ishft(j_hole-1,5)-1 + + psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole) + + l_particle=ishft(k-1,-5)+1 + k_particle=k-ishft(l_particle-1,5)-1 + psi_CIS_holes(key) = i + psi_CIS_particl(key) = k + psi_CIS_spin(key) = 2 + + psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle) + enddo + enddo + + !Building the CIS-matrix + double precision :: diag_H_mat_elem + + do key=1,size_psi_CIS + H_CIS(key,key)=diag_H_mat_elem(psi_CIS(1,1,key),N_int) + + do a=key+1,size_psi_CIS + call i_H_j(psi_CIS(1,1,a),psi_CIS(1,1,key),N_int,hij) + + H_CIS(key,a)=hij + H_CIS(a,key)=hij + enddo + enddo + + END_PROVIDER + + + BEGIN_PROVIDER[double precision, eigenvalues_CIS,(n_state_CIS)] + &BEGIN_PROVIDER[double precision, coefs_CIS, (size_psi_CIS,n_state_CIS)] + use bitmasks + + BEGIN_DOC + !the first states of the CIS matrix + END_DOC + + implicit none + integer :: i,j,k + double precision :: eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS) + double precision :: coefs_tmp(size_psi_CIS) + double precision :: s2 + + + !Diagonalisation of CIS-matrix + call lapack_diag(eigvalues,eigvectors,H_CIS,size_psi_CIS,size_psi_CIS) + + do i = 1,n_state_CIS + eigenvalues_CIS(i) = eigvalues(i) + + do k=1,size_psi_CIS + + if (dabs(eigvectors(k,i)).ge.10.d-2) then + write(11,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i) + write(11,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(11,*),'' + endif + coefs_tmp(k) = eigvectors(k,i) + coefs_CIS(k,i)=eigvectors(k,i) + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + enddo + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D,(size_psi_CIS,n_state_CIS)] +&BEGIN_PROVIDER [double precision, overlap_D ] + use bitmasks + + BEGIN_DOC + !The first states of the CIS matrix dressed by the doubles + END_DOC + implicit none + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision :: overlap,max_overlap,s2 + integer :: i_overlap,i,j,k + allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS)) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) + do i = 1,n_state_CIS + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix_doub(j,k) += H_CIS(j,k) + enddo + enddo + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + + ! state following + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo + ! print*,'overlap = ',max_overlap + overlap_D=max_overlap + do k = 1,size_psi_CIS + eigenvectors_CIS_dress_D(k,i) = eigvectors(k,i_overlap) + if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then + write(12,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap) + write(12,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(12,*),'' + endif + enddo + call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D(1,i),size_psi_CIS,size_psi_CIS,s2) + s_2_CIS_dress_D(i) = s2 + eigenvalues_CIS_dress_D(i) = eigvalues(i_overlap) + enddo + + END_PROVIDER + + + + BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D_dt,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D_dt,(n_state_CIS)] +&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D_dt,(size_psi_CIS,n_state_CIS) +&BEGIN_PROVIDER [double precision, overlap_Ddt] + use bitmasks + + BEGIN_DOC + !The first states of the CIS matrix dressed by the doubles and the disconnected triples + END_DOC + implicit none + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision :: overlap,max_overlap,s2 + integer :: i_overlap,i,j,k + allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS)) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) + do i = 1,n_state_CIS + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles + + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix_doub(j,k) += H_CIS(j,k) + enddo + delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j) + enddo + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + + ! state following + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo + ! print*,'overlap = ',max_overlap + overlap_Ddt=max_overlap + do k = 1,size_psi_CIS + eigenvectors_CIS_dress_D_dt(k,i) = eigvectors(k,i_overlap) + if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then + write(13,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap) + write(13,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k) + write(13,*),'' + endif + enddo + call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D_dt(1,i),size_psi_CIS,size_psi_CIS,s2) + s_2_CIS_dress_D_dt(i) = s2 + eigenvalues_CIS_dress_D_dt(i) = eigvalues(i_overlap) + enddo + + END_PROVIDER + +!BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_tot,(n_state_CIS)] +!BEGIN_PROVIDER [double precision, s_2_CIS_dress_tot,(n_state_CIS)] +!BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_tot,(size_psi_CIS,n_state_CIS)] +!BEGIN_PROVIDER [double precision, overlap_tot] + +!BEGIN_DOC +!!The first states of the CIS matrix dressed by the doubles +!END_DOC +!implicit none +!double precision,allocatable :: delta_H_matrix_doub(:,:) +!double precision,allocatable :: eigvalues(:),eigvectors(:,:) +!double precision,allocatable :: delta_H_trip(:,:) +!double precision :: overlap,max_overlap,s2,average_eigvalue +!integer :: i_overlap,i,j,k,m +!allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS),delta_H_trip(size_psi_CIS,size_psi_CIS) ) +!allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) +! do i = 1,n_state_CIS +! call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles +! call dress_T_con(eigenvalues_CIS(i),delta_H_trip,size_psi_CIS) +! do j = 1,size_psi_CIS +! do k = 1,size_psi_CIS +! delta_H_matrix_doub(j,k) += H_CIS(j,k) +! delta_H_matrix_doub(j,k) += delta_H_trip(j,k) +! enddo +! delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j) +! enddo +! call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) +! do m=1,n_state_CIS +! write(12,*),'m,eigvalues(m)',m,eigvalues(m) +! enddo +! ! state following +! max_overlap = 0.d0 +! do k = 1, size_psi_CIS +! overlap = 0.d0 +! do j = 1,size_psi_CIS +! overlap += eigvectors(j,k)*coefs_CIS(j,i) +! enddo +! if(dabs(overlap).gt.max_overlap)then +! max_overlap = dabs(overlap) +! i_overlap = k +! endif +! ! +! enddo +!! print*,'overlap = ',max_overlap +! overlap_tot=max_overlap +! do k = 1,size_psi_CIS +! eigenvectors_CIS_dress_tot(k,i) = eigvectors(k,i_overlap) +! enddo +! call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_tot(1,i),size_psi_CIS,size_psi_CIS,s2) +! s_2_CIS_dress_tot(i) = s2 +! eigenvalues_CIS_dress_tot(i) = eigvalues(i_overlap) +! enddo + +!END_PROVIDER + + + + + + + + + BEGIN_PROVIDER [double precision, diag_elements, (size_psi_CIS)] + use bitmasks + + + BEGIN_DOC + !Array of the energy of the CIS determinants ordered in the CIS matrix + END_DOC + + implicit none + double precision :: hij + integer :: i + + do i = 1, size_psi_CIS + call i_H_j(psi_CIS(1,1,i),psi_CIS(1,1,i),N_int,hij) + diag_elements(i) = hij + enddo + + END_PROVIDER diff --git a/src/CIS_dressed/MP2.irp.f b/src/CIS_dressed/MP2.irp.f new file mode 100644 index 00000000..cfcb3369 --- /dev/null +++ b/src/CIS_dressed/MP2.irp.f @@ -0,0 +1,1611 @@ + BEGIN_PROVIDER [double precision, MP2_corr_energy] + &BEGIN_PROVIDER [double precision, p_imp,(elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, h_imp,(n_core_cis+1:elec_alpha_num)] + &BEGIN_PROVIDER [double precision, hp_imp,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + + BEGIN_DOC + !Calculation of the MP2 correlation energy (MP2_corr_energy) + !and calculation of the contribution of the disconnected Triples on the + !Singles, via the impossible (p_imp, h_imp, hp_imp) + END_DOC + + implicit none + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + double precision :: direct,exchg,hij !calculating direct, exchange and total contribution + double precision :: get_mo_bielec_integral + double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l + double precision :: delta_e_ik,delta_e_ikj + double precision :: delta_e !delta epsilons + + + + MP2_corr_energy=0.d0 + + do k =elec_beta_num + 1, n_act_cis + p_imp(k)=0.d0 + enddo + + + do i=n_core_cis+1,elec_alpha_num + h_imp(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis + hp_imp(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for MP2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for MP2_corr_energy and a part of the contribution to p_imp and h_imp + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + MP2_corr_energy=MP2_corr_energy+2*hij + p_imp(k)=p_imp(k)+hij + h_imp(i)=h_imp(i)+hij + + enddo + + !remaining same spin contribution for p_imp + do l=elec_alpha_num+1,k-1 + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + p_imp(k)=p_imp(k)+hij + + enddo + enddo + + !remaining same spin contribution for h_imp + do j=n_core_cis+1,i-1 + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + h_imp(i)=h_imp(i)+hij + + enddo + enddo + + !same spin contribution for hp_imp + do j=n_core_cis+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg =get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg)/delta_e + + hp_imp(i,k)=hp_imp(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct/delta_e + + MP2_corr_energy=MP2_corr_energy+hij + p_imp(k)=p_imp(k)+hij + h_imp(i)=h_imp(i)+hij + hp_imp(i,k)=hp_imp(i,k)+hij + + enddo + enddo + enddo + enddo + + print*,'MP2 correlation energy=',MP2_corr_energy + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, EN2_corr_energy] + &BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)] + &BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + + BEGIN_DOC + !Calculation of the EN2 correlation energy (EN2_corr_energy) + !and calculation of the contribution of the disconnected Triples on the + !Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN) + END_DOC + + implicit none + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + double precision :: direct,exchg,hij !calculating direct, exchange and total contribution + double precision :: get_mo_bielec_integral + double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l + double precision :: delta_e_ik,delta_e_ikj + double precision :: delta_e !delta epsilons + double precision :: delta_e_tmp,H_jj_total + integer, allocatable :: key_in(:,:) + integer, allocatable :: key_out(:,:) + integer :: ispin1,ispin2,i_ok + allocate(key_in(N_int,2)) + allocate(key_out(N_int,2)) + + print*,'EN2_corr_energy' + + EN2_corr_energy=0.d0 + +do i=n_core_cis+1,elec_alpha_num + h_imp_EN(i)=0.d0 + do k =elec_beta_num + 1, n_act_cis + p_imp_EN(k)=0.d0 + hp_imp_EN(i,k)=0.d0 + enddo +enddo + + do i=n_core_cis+1,elec_alpha_num +! h_imp_EN(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis +! hp_imp_EN(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for EN2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + ! ispin1 = 1 + ! ispin2 = 1 + ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + ! delta_e_tmp = HF_energy - H_jj_total(key_out) + ! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then + ! print*,'same spin first' + ! print*,'delta_e_SSS = ',delta_e_tmp-delta_e + ! stop + ! endif + ! + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + ! ispin1 = 1 + ! ispin2 = 1 + ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + ! delta_e_tmp = HF_energy - H_jj_total(key_out) + + EN2_corr_energy=EN2_corr_energy+2*hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + + p_imp_EN(l)=p_imp_EN(l)+hij + h_imp_EN(j)=h_imp_EN(j)+hij + + enddo +enddo +! !remaining same spin contribution for p_imp_EN + +! do l=elec_alpha_num+1,k-1 +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=e_i + e_j - e_k -e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) + + + + + +!! ispin1 = 1 +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! call print_key(key_out) +!! call print_key(ref_bitmask) +!! print*,i,k,j,l +!! print*,'same spin middle' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! print*,'delta_e_SSS = ',delta_e_tmp,delta_e +!! stop +!! endif + +! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + +! hij=(direct-exchg)*(direct-exchg) +! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + +! p_imp_EN(k)=p_imp_EN(k)+hij + +! enddo +! enddo + +! !remaining same spin contribution for h_imp_EN +! do j=n_core_cis+1,i-1 +! e_j=diagonal_Fock_matrix_mo(j) +! delta_e_ikj=delta_e_ik+e_j + +! do l=k+1,n_act_cis +! e_l=diagonal_Fock_matrix_mo(l) +! delta_e=delta_e_ikj-e_l +! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & +! mo_bielec_integral_jj_anti(k,l) +! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & +! mo_bielec_integral_jj_anti(i,l) + & +! mo_bielec_integral_jj_anti(j,k) + & +! mo_bielec_integral_jj_anti(j,l) +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! print*,'same spin third ' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! stop +!! endif + +! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + +! hij=(direct-exchg)*(direct-exchg) +! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + +! h_imp_EN(i)=h_imp_EN(i)+hij + +! enddo +! enddo + + !same spin contribution for hp_imp_EN + + do j=n_core_cis+1,elec_alpha_num + if(j==i)cycle + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + if(l==k)cycle + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) +!!! ispin1 = 1 +!!! ispin2 = 1 +!!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!!! print*,'same spin fourth' +!!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!!! call print_key(key_out) +!!! call print_key(ref_bitmask) +!!! print*,i,k,j,l +!!! stop +!!! endif + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij = 0.5d0 * (-delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij)) + + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj(i,j) - & + mo_bielec_integral_jj(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(j,l) + & + mo_bielec_integral_jj(i,l) + & + mo_bielec_integral_jj(j,k) +!! ispin1 = 2 +!! ispin2 = 1 +!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) +!! delta_e_tmp = HF_energy - H_jj_total(key_out) +!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then +!! print*,'different spin ' +!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e +!! stop +!! endif + + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + + EN2_corr_energy=EN2_corr_energy+hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + enddo + enddo + + print*,'EN correlation energy=',EN2_corr_energy + print*,'EN correlation energy=',EN2_corr_energy + + END_PROVIDER + + + + + BEGIN_PROVIDER [integer, size_psi_CIS] + + BEGIN_DOC + !Definition of the size of the CIS vector + END_DOC + + implicit none + size_psi_CIS=(elec_alpha_num-n_core_cis)*(n_act_cis-elec_alpha_num)*2+1 !!!Is this still correct for open shell systems??? + print*,'size_psi_CIS = ',size_psi_CIS + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, diag_elements_sorted, (size_psi_CIS)] + &BEGIN_PROVIDER [integer , diag_elements_sorted_index, (size_psi_CIS)] + BEGIN_DOC + !Array of the energy of the CIS determinants sorted by energy and + !Index in the CIS matrix + END_DOC + + implicit none + integer :: iorder(size_psi_CIS),i + + print*,'diag_elements_sorted' + + do i = 1,size_psi_CIS + diag_elements_sorted(i) = diag_elements(i) + iorder(i) = i + enddo + + call dsort(diag_elements_sorted,iorder,size_psi_CIS) + + do i = 1,size_psi_CIS + diag_elements_sorted(i) = diag_elements(iorder(i)) + diag_elements_sorted_index(i) = iorder(i) + enddo + + END_PROVIDER + + + + + + + BEGIN_PROVIDER [double precision, eigenvalues_dressed_CIS,(n_state_CIS)] + + BEGIN_DOC + !The first states of the dressed CIS matrix + END_DOC + + implicit none + integer :: i,j,k + double precision :: s2,e_corr,tmp + double precision,allocatable :: coefs_tmp(:) + double precision,allocatable :: eigvalues(:),eigvectors(:,:) + double precision,allocatable :: delta_H_trip(:,:) + double precision,allocatable :: delta_H_matrix_doub(:,:) + double precision,allocatable :: delta_H_matrix_dpdiscon(:,:) + double precision,allocatable :: delta_H_matrix(:,:) + double precision :: eigenvalues_dressed_CIS_D(n_state_CIS) + double precision :: eigenvalues_dressed_CIS_DT(n_state_CIS) + double precision,allocatable :: eigvald(:),eigvecd(:,:) + allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS),delta_H_trip(size_psi_CIS,size_psi_CIS),delta_H_matrix_doub(size_psi_CIS,size_psi_CIS),delta_H_matrix(size_psi_CIS,size_psi_CIS)) + allocate (eigvald(size_psi_CIS),eigvecd(size_psi_CIS,size_psi_CIS),coefs_tmp(size_psi_CIS),delta_H_matrix_dpdiscon(size_psi_CIS,size_psi_CIS) ) + + print*,'eigenvalues_dressed_CIS' + + provide coefs_CIS + do i = 1,n_state_CIS + + write(66,*),'i=',i + + call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix,size_psi_CIS) !dressing of the Doubles + do j = 1,size_psi_CIS + do k = 1,size_psi_CIS + delta_H_matrix(j,k)=delta_H_matrix(j,k)+H_CIS(j,k) + + delta_H_matrix_doub(j,k)=delta_H_matrix(j,k) + delta_H_matrix(j,k)=delta_H_matrix(j,k)+delta_H_trip(j,k) + delta_H_matrix_dpdiscon(j,k)=delta_H_matrix(j,k) + enddo + delta_H_matrix(j,j)=delta_H_matrix(j,j)+dress_T_discon_array_CIS(j) !dressing by the disconnected Triples + enddo + + + + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) + double precision :: overlap + double precision :: max_overlap + integer :: i_overlap + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + eigenvalues_dressed_CIS_D(i) = eigvalues(i_overlap) + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) +! print*,'s2(D)=',s2 +! print*,'' +! print*,'diagonal elements sorted(i) =',diag_elements_sorted(i) +! print*,'eigenvalues of the CIS(i) =',eigenvalues_CIS(i) +! print*,'eigenvalues dressed by Doubles (D) =',eigenvalues_dressed_CIS_D(i) + +! write(12,*),'' +! write(12,*),'i=',i +! write(12,*),'eigenvalues of the CIS(i) =',eigenvalues_CIS(i) +! write(12,*),'diagonal elements sorted(i)=',diag_elements_sorted(i) +! write(12,*),'' +! write(12,*),'s2(D)=',s2 +! write(12,*),'eigenvalues dressed by Doubles (D) =',eigenvalues_dressed_CIS_D(i) + + write(66,*),'CIS-HF ',eigenvalues_CIS(i)-eigenvalues_CIS(1) + write(66,*),'Doub-HF ',eigenvalues_dressed_CIS_D(i)-eigenvalues_dressed_CIS_D(1) + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix_dpdiscon,size_psi_CIS,size_psi_CIS) + max_overlap = 0.d0 + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + + + eigenvalues_dressed_CIS_DT(i) = eigvalues(i_overlap) + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + enddo + + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + +! write(12,*),'s2(D+cT)=',s2 +! write(12,*),'eigenvalues dressed by D and connected Triples=',eigenvalues_dressed_CIS_DT(i) +! print*,'s2(DcT)=',s2 +! print*,'eigenvalues dressed D and cT=',eigenvalues_dressed_CIS_DT(i) + + write(66,*),'D+cT-HF ',eigenvalues_dressed_CIS_DT(i)-eigenvalues_dressed_CIS_DT(1) + + + + call lapack_diag(eigvalues,eigvectors,delta_H_matrix,size_psi_CIS,size_psi_CIS) + do k = 1, size_psi_CIS + overlap = 0.d0 + do j = 1,size_psi_CIS + overlap += eigvectors(j,k)*coefs_CIS(j,i) + enddo + if(dabs(overlap).gt.max_overlap)then + max_overlap = dabs(overlap) + i_overlap = k + endif + ! + enddo +! print*,'i_overlap = ',i_overlap + + + eigenvalues_dressed_CIS(i) = eigvalues(i_overlap) + + + do j = 1,size_psi_CIS + coefs_tmp(j) = eigvectors(j,i_overlap) + enddo + call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2) + +! print*,'s2(Tot)=',s2 +! print*,'eigenvalues dressed =',eigenvalues_dressed_CIS(i) + +! write(12,*),'' +! write(12,*),'s2(tot)=',s2 +! write(12,*),'eigenvalues of dressed CIS(i)=',eigenvalues_dressed_CIS(i) + + write(66,*)'Eigval-HF ',eigenvalues_dressed_CIS(i)-eigenvalues_dressed_CIS(1) + enddo + + deallocate(eigvalues,eigvectors,delta_H_trip,delta_H_matrix_doub,delta_H_matrix_dpdiscon) + deallocate (eigvald,eigvecd,coefs_tmp) +!call system ( "mkdir results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp ../newcode/MP2.irp.f results_$(date +%d%m%y_%I%M)" ) +!call system ( "mv fort.10 Ecorr.out" ) +!call system ( "mv fort.11 CIS.out" ) +!call system ( "mv fort.12 Dressing.out" ) +!call system ( "mv fort.66 delta_e_states" ) +!call system ( "cp Ecorr.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp CIS.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp Dressing.out results_$(date +%d%m%y_%I%M)" ) +!call system ( "cp delta_e_states results_$(date +%d%m%y_%I%M)" ) +!call system ( "rm delta_e_states" ) +!call system ( "rm Ecorr.out" ) +!call system ( "rm CIS.out" ) +!call system ( "rm Dressing.out" ) + + END_PROVIDER + + + BEGIN_PROVIDER [double precision, dress_T_discon,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + &BEGIN_PROVIDER [double precision, dress_T_discon_array_CIS,(size_psi_CIS)] + + BEGIN_DOC + !Calculation of the dressing by the disconnected Triples, via the impossible + END_DOC + + implicit none + integer :: i,k !variables for going over the occupied (i) and virutal (k) + integer :: key !key for CIS-matrix + + print*,'dress_T_discon' + + if(mp2_dressing)then + dress_T_discon_array_CIS(1)=MP2_corr_energy + + do i=n_core_cis+1,elec_alpha_num + do k=elec_alpha_num+1,n_act_cis + key=psi_CIS_adress(i,k) + + dress_T_discon(i,k)=MP2_corr_energy-p_imp(k)-h_imp(i)+hp_imp(i,k) + + dress_T_discon_array_CIS(key) = dress_T_discon(i,k) + dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) + + enddo + enddo + + else !EN Dressing + dress_T_discon_array_CIS(1)=EN2_corr_energy + + do i=n_core_cis+1,elec_alpha_num + do k=elec_alpha_num+1,n_act_cis + key=psi_CIS_adress(i,k) + + dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k) + +!print*,'i,k,key',i,k,key +!print*,'EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k)',EN2_corr_energy,p_imp_EN(k),h_imp_EN(i),hp_imp_EN(i,k) +!print*,'dress_T_discon',dress_T_discon(i,k) + + dress_T_discon_array_CIS(key) = dress_T_discon(i,k) + dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) + + enddo + enddo + end if + + END_PROVIDER + + + subroutine dress_by_doubles(ref_energy,coefs,delta_H_matrix,ndet) + + use bitmasks + implicit none + integer, intent(in) :: ndet + double precision, intent(in) :: ref_energy + double precision, intent(in) ::coefs(size_psi_CIS) + double precision, intent(out) :: delta_H_matrix(ndet,ndet) + + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + integer :: key !key for CIS-matrix + integer :: a,b,i_ok !control variables + integer(bit_kind) :: key_out(N_int,2) !Doubles + integer :: ispin1,ispin2 !spin variables + integer :: degree,i_count_connected + double precision :: hmD,hij,hij_array(ndet),hij_index(ndet) + double precision :: delta_hij + double precision :: delta_e , e_double!delta epsilons + double precision :: accu_H_mD,diag_H_mat_elem + double precision :: epsilon_D,s2 + + print*,'dress_by_doubles' + + double precision :: accu,accu_2 + accu = 0.d0 + accu_2 = 0.d0 + do i =1,ndet + accu += coefs(i) * coefs(i) + enddo + if (dabs(accu-1.d0) > 1.d-10)then + print*,'coefficients not normalized !!' + print*,accu + !stop + endif + + accu = 0.d0 + if(standard_doubles)then + delta_H_matrix = 0.d0 + + do i=n_core_cis+1,elec_alpha_num + + do k=elec_alpha_num+1,n_act_cis + + !same spin contribution + do j=i+1,elec_alpha_num + + do l=k+1,n_act_cis + !Alpha + ispin1=1 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + !Beta + ispin1=2 + ispin2=2 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a)*delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b)*delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + + do l=elec_beta_num+1,n_act_cis + ispin1=2 + ispin2=1 + + hij_array=0.d0 + + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + delta_e = 1.d0 / (ref_energy - e_double) + + i_count_connected = 0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + enddo + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a)*delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b)*delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + enddo + enddo + + else !2x2 Matrix diagonalisation + + + delta_H_matrix = 0.d0 + + do i=n_core_cis+1,elec_alpha_num + + do k=elec_alpha_num+1,n_act_cis + + !same spin contribution + do j=i+1,elec_alpha_num + + do l=k+1,n_act_cis + !Alpha + ispin1=1 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + +! if (dabs(hmD).le.1.d-7) cycle + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + !Beta + ispin1=2 + ispin2=2 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + +! if (dabs(hmD).le.1.d-7) cycle + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + + do l=elec_beta_num+1,n_act_cis + ispin1=2 + ispin2=1 + + call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) + + e_double =diag_H_mat_elem(key_out,N_int) + + i_count_connected = 0 + accu_H_md=0.d0 + + do a=2, ndet + call get_excitation_degree(key_out,psi_CIS(1,1,a),degree,N_int) + + if (degree /= 1 .and.degree /= 2) cycle + + call i_H_j(key_out,psi_CIS(1,1,a),N_int,hij) + + if (dabs(hij).le.1.d-10) cycle + + i_count_connected = i_count_connected +1 + + hij_array(i_count_connected) = hij + hij_index(i_count_connected) = a + + accu_H_mD=accu_H_mD+(coefs(a))*hij + + enddo + + hmD=accu_H_mD + + + epsilon_D=0.5*((e_double-ref_energy)-sqrt((ref_energy -e_double)*(ref_energy -e_double)+4*hmD*hmD)) + if (dabs(hmD*hmD).le.1.d-10.or.dabs(epsilon_D).le.1.d-10) cycle + accu += epsilon_D + accu_2 += hmD*hmD + + delta_e = epsilon_D/(hmD*hmD) +! e_double =1.d0/(ref_energy - H_jj_total_diff(key_out,ref_bitmask,HF_energy,N_int)) +! if(dabs(e_double - delta_e)/(dabs(e_double)).gt.1.d-3)then +! print*,delta_e,e_double +! endif + + do a=1,i_count_connected + delta_hij = hij_array(a)*hij_array(a) * delta_e + + delta_H_matrix(hij_index(a),hij_index(a))=delta_H_matrix(hij_index(a),hij_index(a))+delta_hij + + do b=a+1,i_count_connected + delta_hij = hij_array(a)*hij_array(b) * delta_e + + delta_H_matrix(hij_index(a),hij_index(b))=delta_H_matrix(hij_index(a),hij_index(b))+delta_hij + delta_H_matrix(hij_index(b),hij_index(a))=delta_H_matrix(hij_index(b),hij_index(a))+delta_hij + + enddo + enddo + + enddo + enddo + enddo + enddo + + + end if + + end + + subroutine dress_T_con(ref_energy,delta_H_trip,ndet) + + BEGIN_DOC + !Generating all the Triples and, in the loops, the connected Singles + END_DOC + + implicit none + double precision :: e_i,e_r,e_j,e_s,e_k,e_t !epsilons of occupied and virtuals + double precision :: delta_e_ir,delta_e_irj,delta_e_irjs,delta_e_irjsk + double precision :: delta_e,energy !delta epsilons + double precision :: phi_tmp,phi_doub + double precision :: direct,exchg,get_mo_bielec_integral + double precision :: phi_aaa(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab_alpha(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab_beta(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis),phi_aab(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis) + integer :: i,j,k,r,s,t !holes and particles of the Triples + integer :: occ,vir,occ_tmp,vir_tmp !hole and particle of the Singles + integer :: a,b,c,d !variables of the conected Singels + integer :: m,key_ir,key_js,key_is,key_jr,key_kt,key_ks,key_kr,key_it,key_jt,key_tmp + + integer, intent(in) :: ndet + double precision, intent(in) :: ref_energy + double precision, intent(out) :: delta_H_trip(ndet,ndet) + + delta_H_trip=0.d0 +!do i=5,6 + do i=n_core_cis+1,elec_alpha_num + e_i=diagonal_Fock_matrix_mo(i) + !r=7 + do r=elec_alpha_num+1,n_act_cis + e_r=diagonal_Fock_matrix_mo(r) + delta_e_ir=e_i-e_r + key_ir=psi_CIS_adress(i,r) + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_irj=delta_e_ir+e_j + do s=r+1,n_act_cis + e_s=diagonal_Fock_matrix_mo(s) + delta_e_irjs=delta_e_irj-e_s + !alpha-alpha-alpha + do k=j+1,elec_alpha_num + e_k=diagonal_Fock_matrix_mo(k) + delta_e_irjsk=delta_e_irjs+e_k + do t=s+1,n_act_cis + e_t=diagonal_Fock_matrix_mo(t) + delta_e=delta_e_irjsk-e_t + energy=1.d0/(ref_energy-delta_e) + occ=i + a=j + b=k + vir=r + c=s + d=t + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + phi_tmp=direct-exchg + do occ=i,k + if (occ==i) then + a=j + b=k + else if (occ==j) then + a=i + b=k + !else if (occ==k) then + ! a=i + ! b=j + !else + ! cycle + end if + do vir=r,t + if (vir==r) then + c=s + d=t + else if (vir==s) then + c=r + d=t + ! else if (vir==t) then + ! c=r + ! d=t + ! else + ! cycle + end if + key_tmp=psi_CIS_adress(occ,vir) + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + if (occ==j .and. vir==s .or. occ==k .and. vir==t ) then + phi_aaa(occ,vir)=direct-exchg + else if (occ==i .and. vir==r) then + phi_aaa(occ,vir)=0.d0 + !else if (occ/=i .or. occ/=j .or. occ/=k .or. vir/=r .or. vir/=s .or. vir /=t) then + ! phi_aaa(occ,vir)=0.d0 + else + phi_aaa(occ,vir)=-direct+exchg + endif + phi_doub=phi_tmp*phi_aaa(occ,vir)*energy + delta_H_trip(key_ir,key_tmp)=delta_H_trip(key_ir,key_tmp)+phi_doub + delta_H_trip(key_ir+1,key_tmp+1)=delta_H_trip(key_ir+1,key_tmp+1)+phi_doub + + ! if (phi_doub.ge.1.d-6) then + ! print*,'alpha alpha' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + ! endif + + ! if (key_ir==key_tmp) then + ! print*,'ir=tmp' + ! print*,'i,r,j,s,k,t',i,r,j,s,k,t + ! print*,'occ,vir',occ,vir + ! !stop + ! end if + + enddo + enddo + enddo + enddo + !alpha-alpha-beta + do k=n_core_cis+1,elec_alpha_num + e_k=diagonal_Fock_matrix_mo(k) + delta_e_irjsk=delta_e_irjs+e_k + do t=elec_alpha_num+1,n_act_cis + e_t=diagonal_Fock_matrix_mo(t) + delta_e=delta_e_irjsk-e_t + energy=1.d0/(ref_energy-delta_e) + occ=i + a=j + b=k + vir=r + c=s + d=t + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + phi_tmp=direct-exchg + do occ=i,k + if (occ==i) then + a=j + b=k + else if (occ==j) then + a=i + b=k + else if (occ==k) then + a=i + b=j + else + cycle + end if + do vir=r,t + if (vir==r) then + c=s + d=t + else if (vir==s) then + c=r + d=t + else if (vir==t) then + c=r + d=s + else + cycle + end if + key_tmp=psi_CIS_adress(occ,vir) + direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + if (occ==j .and. vir==s .or. occ==k .and. vir==t ) then + phi_aab(occ,vir)=direct + else if (occ==i .and. vir==r .or. occ==i .and. vir==t .or. occ==j .and. vir==t .or. occ==k .and. vir==r .or. occ==k .and. vir==s) then + phi_aab(occ,vir)=0.d0 + !else if (occ/=i .or. occ/=j .or. occ/=k .or. vir/=r .or. vir/=s .or. vir /=t) then + ! phi_aab(occ,vir)=0.d0 + + else + phi_aab(occ,vir)=-direct + endif + + phi_doub=phi_tmp*phi_aab(occ,vir)*energy + delta_H_trip(key_ir,key_tmp+1)=delta_H_trip(key_ir,key_tmp+1)+phi_doub + delta_H_trip(key_ir+1,key_tmp)=delta_H_trip(key_ir,key_tmp+1)+phi_doub + ! if (phi_doub.ge.1.d-6) then + ! if (key_ir==key_tmp) then + ! print*,'ir=tmp beta ' + ! print*,'i,r,j,s,k,t',i,r,j,s,k,t + ! print*,'occ,vir',occ,vir + ! !stop + ! end if + ! end if + + ! if (phi_doub.ge.1.d-6) then + ! if (vir==r .or. vir==s) then + ! print*,'alpha' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + + ! else if (vir==t) then + ! print*,'beta' + ! print*,'occ,vir,key_tmp',occ,vir,key_tmp + + ! endif + ! print*,'vir=?',vir + ! endif + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + +!delta_H_trip=0.d0 + + + +!!generating the Singles included in the Triples +!! do m=1,n_state_CIS +!do i=n_core_cis+1,elec_alpha_num +! e_i=diagonal_Fock_matrix_mo(i) + + + +!!do r=elec_alpha_num+1,n_state_CIS +! do r=elec_alpha_num+1,n_act_cis +! e_r=diagonal_Fock_matrix_mo(r) +! delta_e_ir=e_i-e_r +!key_ir=psi_CIS_adress(i,r) + +! !alpha-alpha-x (=beta-beta-x) +! do j=i+1,elec_alpha_num +! e_j=diagonal_Fock_matrix_mo(j) +! delta_e_irj=delta_e_ir+e_j +!key_jr=psi_CIS_adress(j,r) +! do s=r+1,n_act_cis +! e_s=diagonal_Fock_matrix_mo(s) +! delta_e_irjs=delta_e_irj-e_s + +! key_is=psi_CIS_adress(i,s) +! key_js=psi_CIS_adress(j,s) + + +! !alpha-alpha-alpha (=beta-beta-beta) +! do k=j+1,elec_alpha_num +! e_k=diagonal_Fock_matrix_mo(k) +! delta_e_irjsk=delta_e_irjs+e_k +! key_kr=psi_CIS_adress(k,r) +! key_ks=psi_CIS_adress(k,s) + +! do t=s+1,n_act_cis +! e_t=diagonal_Fock_matrix_mo(t) +! delta_e=delta_e_irjsk-e_t + +! key_it=psi_CIS_adress(i,t) +! key_jt=psi_CIS_adress(j,t) +! key_kt=psi_CIS_adress(k,t) + + + +! !generating the connected singles +! do occ=i,k +! if (occ==i) then +! a=j +! b=k +! else if (occ==j) then +! a=i +! b=k +! else if (occ==k) then +! a=i +! b=j +! else +! cycle +! end if + +! do vir=r,t +! if (vir==r) then +! c=s +! d=t +! else if (vir==s) then +! c=r +! d=t +! else if (vir==t) then +! c=r +! d=s +! else +! cycle +! end if + +! if (occ==i .and. vir==r .or. occ==j .and. vir==s .or. occ==k .and. vir==t ) then +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aaa(occ,vir)=direct-exchg +! else +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aaa(occ,vir)=-direct+exchg +! endif + +! ! phi_aaa(occ,vir)=phi_aaa(occ,vir)+phi_tmp + +! enddo +! enddo + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_js)=delta_H_trip(key_ir,key_js)+phi_tmp +! delta_H_trip(key_ir+1,key_js+1)=delta_H_trip(key_ir+1,key_js+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_kt)=delta_H_trip(key_ir,key_kt)+phi_tmp +! delta_H_trip(key_ir+1,key_kt+1)=delta_H_trip(key_ir+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(j,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_jt)=delta_H_trip(key_ir,key_jt)+phi_tmp +! delta_H_trip(key_ir+1,key_jt+1)=delta_H_trip(key_ir+1,key_jt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,r)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_ks)=delta_H_trip(key_ir,key_ks)+phi_tmp +! delta_H_trip(key_ir+1,key_ks+1)=delta_H_trip(key_ir+1,key_ks+1)+phi_tmp + + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jr)=delta_H_trip(key_is,key_jr)+phi_tmp +! delta_H_trip(key_is+1,key_jr+1)=delta_H_trip(key_is+1,key_jr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kt)=delta_H_trip(key_is,key_kt)+phi_tmp +! delta_H_trip(key_is+1,key_kt+1)=delta_H_trip(key_is+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(j,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jt)=delta_H_trip(key_is,key_jt)+phi_tmp +! delta_H_trip(key_is+1,key_jt+1)=delta_H_trip(key_is+1,key_jt+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,s)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kr)=delta_H_trip(key_is,key_kr)+phi_tmp +! delta_H_trip(key_is+1,key_kr+1)=delta_H_trip(key_is+1,key_kr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_js)=delta_H_trip(key_it,key_js)+phi_tmp +! delta_H_trip(key_it+1,key_js+1)=delta_H_trip(key_it+1,key_js+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_jr)=delta_H_trip(key_it,key_jr)+phi_tmp +! delta_H_trip(key_it+1,key_jr+1)=delta_H_trip(key_it+1,key_jr+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_ks)=delta_H_trip(key_it,key_ks)+phi_tmp +! delta_H_trip(key_it+1,key_ks+1)=delta_H_trip(key_it+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(i,t)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_it,key_kr)=delta_H_trip(key_it,key_kr)+phi_tmp +! delta_H_trip(key_it+1,key_kr+1)=delta_H_trip(key_it+1,key_kr+1)+phi_tmp + + +! phi_tmp=(phi_aaa(j,r)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_kt)=delta_H_trip(key_jr,key_kt)+phi_tmp +! delta_H_trip(key_jr+1,key_kt+1)=delta_H_trip(key_jr+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,r)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_ks)=delta_H_trip(key_jr,key_ks)+phi_tmp +! delta_H_trip(key_jr+1,key_ks+1)=delta_H_trip(key_jr+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,s)*phi_aaa(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js,key_kt)=delta_H_trip(key_js,key_kt)+phi_tmp +! delta_H_trip(key_js+1,key_kt+1)=delta_H_trip(key_js+1,key_kt+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,s)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js,key_kr)=delta_H_trip(key_js,key_kr)+phi_tmp +! delta_H_trip(key_js+1,key_kr+1)=delta_H_trip(key_js+1,key_kr+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,t)*phi_aaa(k,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jt,key_ks)=delta_H_trip(key_jt,key_ks)+phi_tmp +! delta_H_trip(key_jt+1,key_ks+1)=delta_H_trip(key_jt+1,key_ks+1)+phi_tmp + +! phi_tmp=(phi_aaa(j,t)*phi_aaa(k,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jt,key_kr)=delta_H_trip(key_jt,key_kr)+phi_tmp +! delta_H_trip(key_jt+1,key_kr+1)=delta_H_trip(key_jt+1,key_kr+1)+phi_tmp + +! enddo +! enddo + +! !alpha-alpha-beta (=beta-beta-alpha) +! do k=n_core_cis+1,elec_beta_num +! e_k=diagonal_Fock_matrix_mo(k) +! delta_e_irjsk=delta_e_irjs+e_k + +! key_kr=psi_CIS_adress(k,r) +! key_ks=psi_CIS_adress(k,s) + + + +! do t=elec_beta_num+1,n_act_cis +! e_t=diagonal_Fock_matrix_mo(t) +! delta_e=delta_e_irjsk-e_t + +! key_it=psi_CIS_adress(i,t) +! key_jt=psi_CIS_adress(j,t) +! key_kt=psi_CIS_adress(k,t) + + +! !generating the connected singles alpha (beta) +! do occ=i,j +! if (occ==i) then +! a=j +! b=k +! else if (occ==j) then +! a=i +! b=k +! end if + +! do vir = r,s +! if (vir==r) then +! c=s +! d=t +! else if (vir==s) then +! c=r +! d=t +! end if + +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) + +! if (occ==i .and. vir==r .or. occ==j .and. vir==s) then +! phi_aab_alpha(occ,vir)=direct +! else if (occ==i .and. vir==s .or. occ==j .and. vir==r) then +! phi_aab_alpha(occ,vir)=-1.d0*direct + +! end if + +!! phi_aab_alpha(occ,vir)=phi_aab_alpha(occ,vir)+phi_tmp + +! enddo +! enddo + +! !alpha single alpha single connected +! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_alpha(j,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_js)=delta_H_trip(key_ir,key_js)+phi_tmp +! delta_H_trip(key_ir+1,key_js+1)=delta_H_trip(key_ir+1,key_js+1)+phi_tmp + + +! ! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_alpha(i,s))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! ! delta_H_trip(key_ir,key_is)=delta_H_trip(key_ir,key_is)+phi_tmp +! ! delta_H_trip(key_ir+1,key_is+1)=delta_H_trip(key_ir+1,key_is+1)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(i,s)*phi_aab_alpha(j,r))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_jr)=delta_H_trip(key_is,key_jr)+phi_tmp +! delta_H_trip(key_is+1,key_jr+1)=delta_H_trip(key_is+1,key_jr+1)+phi_tmp + + + + + +! !generating the conncected singles beta (alpha) +! ! do occ=i,k +! ! if (occ==i) then +! ! a=j +! ! b=k +! ! else if (occ==j) then +! ! a=i +! ! b=k +! ! else if (occ==k) then +! ! +! ! a=i +! ! b=j +! ! end if +! ! do vir=r,t +! ! if (vir==r) then +! ! c=s +! ! d=t +! ! else if (vir==s) then +! ! c=r +! ! d=t +! ! else if (vir==k) then +! ! +! ! c=r +! ! d=t +! ! end if + +! occ=k +! a=i +! b=j +! vir=t +! c=r +! d=s +! direct=get_mo_bielec_integral(a,b,c,d,mo_integrals_map) +! exchg =get_mo_bielec_integral(a,b,d,c,mo_integrals_map) + +! phi_aab_beta(occ,vir)=direct-exchg +! ! enddo +! ! enddo +!! phi_aab_beta(occ,vir)=phi_aab_beta(occ,vir)+phi_tmp +! + + + + +! !alpha single beta single connected +! phi_tmp=(phi_aab_alpha(i,r)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_ir,key_kt+1)=delta_H_trip(key_ir,key_kt+1)+phi_tmp +! delta_H_trip(key_ir+1,key_kt)=delta_H_trip(key_ir+1,key_kt)+phi_tmp + +! phi_tmp=(phi_aab_alpha(i,s)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_is,key_kt+1)=delta_H_trip(key_is,key_kt+1)+phi_tmp +! delta_H_trip(key_is+1,key_kt)=delta_H_trip(key_is+1,key_kt)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(j,r)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_jr,key_kt+1)=delta_H_trip(key_jr,key_kt+1)+phi_tmp +! delta_H_trip(key_jr+1,key_kt)=delta_H_trip(key_jr+1,key_kt)+phi_tmp + + +! phi_tmp=(phi_aab_alpha(j,s)*phi_aab_beta(k,t))/(ref_energy-delta_e-eigenvalues_CIS(1)) +! delta_H_trip(key_js+1,key_kt)=delta_H_trip(key_js+1,key_kt)+phi_tmp +! delta_H_trip(key_js,key_kt+1)=delta_H_trip(key_js,key_kt+1)+phi_tmp + + +! enddo +! enddo +! enddo +! enddo +! enddo +!enddo + + END + +subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint) + implicit none + use bitmasks + ! realize the double excitation i-->k (ispin1) + j-->l (ispin2) on key_in + ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible) + integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint + integer(bit_kind), intent(in) :: key_in(Nint,2) + integer, intent(out):: i_ok + integer(bit_kind), intent(out):: key_out(Nint,2) + integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta + integer :: i_test_hole,i_test_particl + key_out = key_in + + k_hole = ishft(i-1,-5)+1 + j_hole = i-ishft(k_hole-1,5)-1 + i_test_hole = ibset(0,j_hole) + if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole) + k_particl = ishft(k-1,-5)+1 + j_particl = k-ishft(k_particl-1,5)-1 + i_test_particl= ibset(0,j_particl) + if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0)then + i_ok = 0 + return + endif + key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl) + + k_hole = ishft(j-1,-5)+1 + j_hole = j-ishft(k_hole-1,5)-1 + i_test_hole = ibset(0,j_hole) + if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole) + k_particl = ishft(l-1,-5)+1 + j_particl = l-ishft(k_particl-1,5)-1 + i_test_particl = ibset(0,j_particl) + if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0)then + i_ok = 0 + return + endif + key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl) + i_ok = 1 + end + diff --git a/src/Hartree_Fock_AOs/README.rst b/src/CIS_dressed/README.rst similarity index 80% rename from src/Hartree_Fock_AOs/README.rst rename to src/CIS_dressed/README.rst index 0418aca5..bb6eeeea 100644 --- a/src/Hartree_Fock_AOs/README.rst +++ b/src/CIS_dressed/README.rst @@ -1,8 +1,3 @@ -=================== -Hartree-Fock Module -=================== - - Needed Modules ============== @@ -12,13 +7,16 @@ Needed Modules * `AOs `_ * `BiInts `_ * `Bitmask `_ +* `Dets `_ * `Electrons `_ * `Ezfio_files `_ +* `Hartree_Fock `_ * `MonoInts `_ * `MOs `_ * `Nuclei `_ * `Output `_ * `Utils `_ +* `DensityMatrix `_ Documentation ============= diff --git a/src/CIS_dressed/cis_dressed.ezfio_config b/src/CIS_dressed/cis_dressed.ezfio_config new file mode 100644 index 00000000..fd9f291a --- /dev/null +++ b/src/CIS_dressed/cis_dressed.ezfio_config @@ -0,0 +1,6 @@ +cis_dressed + n_state_cis integer + n_core_cis integer + n_act_cis integer + mp2_dressing logical + standard_doubles logical diff --git a/src/CIS_dressed/density_matrix_suroutine.irp.f b/src/CIS_dressed/density_matrix_suroutine.irp.f new file mode 100644 index 00000000..59e099a8 --- /dev/null +++ b/src/CIS_dressed/density_matrix_suroutine.irp.f @@ -0,0 +1,66 @@ + + subroutine get_dm_from_psi(dets_in,u_in,sze,dim_in,Nint,dm_alpha,dm_beta) + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix + ! + ! dets_in :: bitsrings corresponding to the determinants in the wave function + ! + ! u_in :: coefficients of the wave function + ! + ! sze :: number of determinants in the wave function + ! + ! dim_in :: physical dimension of the array u_in and dets_in + ! + ! Nint :: should be equal to N_int + ! + ! dm_alpha :: alpha one body density matrix + ! + ! dm_beta :: beta one body density matrix + END_DOC + use bitmasks + integer, intent(in) :: sze,dim_in,Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,dim_in) + double precision, intent(in) :: u_in(dim_in) + double precision, intent(out) :: dm_alpha(mo_tot_num,mo_tot_num) + double precision, intent(out) :: dm_beta(mo_tot_num,mo_tot_num) + + integer :: j,k,l + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + dm_alpha = 0.d0 + dm_beta = 0.d0 + + do k=1,sze + call bitstring_to_list(dets_in(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(dets_in(1,2,k), occ(1,2), n_occ_alpha, N_int) + ck = u_in(k) + do l=1,elec_alpha_num + j = occ(l,1) + dm_alpha(j,j) += ck*ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + dm_beta(j,j) += ck*ck + enddo + do l=1,k-1 + call get_excitation_degree(dets_in(1,1,k),dets_in(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(dets_in(1,1,k),dets_in(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + ckl = ck * u_in(l) * phase + if (s1==1) then + dm_alpha(h1,p1) += ckl + dm_alpha(p1,h1) += ckl + else + dm_beta(h1,p1) += ckl + dm_beta(p1,h1) += ckl + endif + enddo + enddo + end diff --git a/src/CIS_dressed/natural_particl_orbitals.irp.f b/src/CIS_dressed/natural_particl_orbitals.irp.f new file mode 100644 index 00000000..532f640f --- /dev/null +++ b/src/CIS_dressed/natural_particl_orbitals.irp.f @@ -0,0 +1,92 @@ + BEGIN_PROVIDER [double precision, particle_natural_orb_CIS_properties, (6,n_state_cis)] +&BEGIN_PROVIDER [double precision, CIS_states_properties, (6,n_state_cis)] + implicit none + BEGIN_DOC +! properties of the natural orbital of the particle of the various n_state_cis eigenvectors of the CIS matrix +! +! You first build the density matrix of the one eigenvector and you take off the Hartree Fock density matrix +! +! particl(i,j)(state = k) == dm(i,j)(Hartree Fock) - dm(i,j)(state = k) +! +! you diagonalize particl(i,j) and the first eigenvector is the natural orbital corresponding to the particl +! +! that is specific to the excitation in the CIS state +! +! particle_natural_orb_CIS_properties(i,1) = +! +! particle_natural_orb_CIS_properties(i,2) = +! +! particle_natural_orb_CIS_properties(i,3) = +! +! particle_natural_orb_CIS_properties(i,5) = +! +! particle_natural_orb_CIS_properties(i,6) = +! +! particle_natural_orb_CIS_properties(i,7) = +! +! CIS_states_properties(i,1:6) = the same but for the hole state i + END_DOC + + integer :: i,j,k,l + double precision :: dm_alpha(mo_tot_num,mo_tot_num) + double precision :: dm_beta(mo_tot_num,mo_tot_num) + double precision :: dm(mo_tot_num,mo_tot_num) + double precision :: eigvalues(mo_tot_num) + double precision :: eigvectors(mo_tot_num,mo_tot_num) + double precision :: accu_n_elec,c_k + + + do i = 1, n_state_cis + print*,' state cis = ',i + call get_dm_from_psi(psi_CIS,coefs_CIS(1,i),size_psi_CIS,size_psi_CIS,N_int,dm_alpha,dm_beta) + + dm = dm_alpha + dm_beta + call get_properties_from_density_matrix(dm,CIS_states_properties(1,i)) + dm = -dm + do k = 1, elec_alpha_num + dm(k,k) += 1.d0 + enddo + do k = 1, elec_beta_num + dm(k,k) += 1.d0 + enddo + call lapack_diag(eigvalues,eigvectors,dm,mo_tot_num,mo_tot_num) + accu_n_elec = 0.d0 + do k = 1, mo_tot_num + accu_n_elec += eigvalues(k) + enddo + do k = 1, mo_tot_num + do l = 1, mo_tot_num + c_k = eigvectors(k,j) * eigvectors(l,j) + particle_natural_orb_CIS_properties(1,i) += c_k * mo_dipole_x(k,l) + particle_natural_orb_CIS_properties(2,i) += c_k * mo_dipole_y(k,l) + particle_natural_orb_CIS_properties(3,i) += c_k * mo_dipole_z(k,l) + particle_natural_orb_CIS_properties(4,i) += c_k * mo_spread_x(k,l) + particle_natural_orb_CIS_properties(5,i) += c_k * mo_spread_y(k,l) + particle_natural_orb_CIS_properties(6,i) += c_k * mo_spread_z(k,l) + enddo + enddo + enddo + +END_PROVIDER +subroutine get_properties_from_density_matrix(dm,properties) + implicit none + double precision, intent(in) :: dm(mo_tot_num,mo_tot_num) + double precision, intent(out) :: properties(6) + integer :: k,l + double precision :: c_k + do k = 1, 6 + properties(k) = 0.d0 + enddo + do k = 1, mo_tot_num + do l = 1, mo_tot_num + c_k = dm(k,l) + properties(1) += c_k * mo_dipole_x(k,l) + properties(2) += c_k * mo_dipole_y(k,l) + properties(3) += c_k * mo_dipole_z(k,l) + properties(4) += c_k * mo_spread_x(k,l) + properties(5) += c_k * mo_spread_y(k,l) + properties(6) += c_k * mo_spread_z(k,l) + enddo + enddo + +end diff --git a/src/CIS_dressed/options.irp.f b/src/CIS_dressed/options.irp.f new file mode 100644 index 00000000..59b67c4f --- /dev/null +++ b/src/CIS_dressed/options.irp.f @@ -0,0 +1,90 @@ +BEGIN_PROVIDER [ integer , n_state_cis ] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_state_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_state_cis(n_state_cis) + else + n_state_cis = 5 + call ezfio_set_cis_dressed_n_state_cis(n_state_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer , n_core_cis] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_core_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_core_cis(n_core_cis) + else + n_core_cis = 0 + call ezfio_set_cis_dressed_n_core_cis(n_core_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer , n_act_cis] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_n_act_cis(has) + if (has) then + call ezfio_get_cis_dressed_n_act_cis(n_act_cis) + else + n_act_cis = mo_tot_num + call ezfio_set_cis_dressed_n_act_cis(n_act_cis) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical , mp2_dressing] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_mp2_dressing(has) + if (has) then + call ezfio_get_cis_dressed_mp2_dressing(mp2_dressing) + else + mp2_dressing = .False. + call ezfio_set_cis_dressed_mp2_dressing(mp2_dressing) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical , standard_doubles] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_standard_doubles(has) + if (has) then + call ezfio_get_cis_dressed_standard_doubles(standard_doubles) + else + standard_doubles = .True. + call ezfio_set_cis_dressed_standard_doubles(standard_doubles) + endif + +END_PROVIDER + diff --git a/src/DensityMatrix/README.rst b/src/DensityMatrix/README.rst index d0734ed9..af3c405c 100644 --- a/src/DensityMatrix/README.rst +++ b/src/DensityMatrix/README.rst @@ -8,6 +8,39 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`iunit_two_body_dm_aa `_ + Temporary files for 2-body dm calculation + +`iunit_two_body_dm_ab `_ + Temporary files for 2-body dm calculation + +`iunit_two_body_dm_bb `_ + Temporary files for 2-body dm calculation + +`one_body_dm_a `_ + Alpha and beta one-body density matrix + +`one_body_dm_b `_ + Alpha and beta one-body density matrix + +`two_body_dm_diag_aa `_ + diagonal part of the two body density matrix + +`two_body_dm_diag_ab `_ + diagonal part of the two body density matrix + +`two_body_dm_diag_bb `_ + diagonal part of the two body density matrix + +`det_coef_provider `_ + Undocumented + +`det_num `_ + Undocumented + +`det_provider `_ + Undocumented + Needed Modules @@ -28,4 +61,5 @@ Needed Modules * `Nuclei `_ * `Output `_ * `Utils `_ +* `DensityMatrix `_ diff --git a/src/DensityMatrix/det_num.irp.f b/src/DensityMatrix/det_num.irp.f index 7c73b0a4..6ba0c1b4 100644 --- a/src/DensityMatrix/det_num.irp.f +++ b/src/DensityMatrix/det_num.irp.f @@ -43,14 +43,14 @@ -0.316337207748718D-01 & /) - do i=1,10 - call write_bitstring( 6, det_provider(1,1,i), N_int ) - enddo - print *, '' - do i=1,10 - call write_bitstring( 6, det_provider(1,2,i), N_int ) - enddo - print *, '' +!do i=1,10 +! call write_bitstring( 6, det_provider(1,1,i), N_int ) +!enddo +!print *, '' +!do i=1,10 +! call write_bitstring( 6, det_provider(1,2,i), N_int ) +!enddo +!print *, '' END_PROVIDER diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index 40fe4b1b..974e123a 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -33,11 +33,12 @@ subroutine resize_H_apply_buffer_det(new_size) integer, intent(in) :: new_size integer(bit_kind), allocatable :: buffer_det(:,:,:) double precision, allocatable :: buffer_coef(:,:) + double precision, allocatable :: buffer_e2(:,:) integer :: i,j,k integer :: Ndet ASSERT (new_size > 0) - allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states) ) + allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states), buffer_e2(new_size,N_states) ) do i=1,min(new_size,H_apply_buffer_N_det) do k=1,N_int @@ -48,9 +49,10 @@ subroutine resize_H_apply_buffer_det(new_size) ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num ) enddo do k=1,N_states - do i=1,min(new_size,H_apply_buffer_N_det) - buffer_coef(i,k) = H_apply_buffer_coef(i,k) - enddo + do i=1,min(new_size,H_apply_buffer_N_det) + buffer_coef(i,k) = H_apply_buffer_coef(i,k) + buffer_e2(i,k) = H_apply_buffer_e2(i,k) + enddo enddo H_apply_buffer_size = new_size @@ -70,20 +72,23 @@ subroutine resize_H_apply_buffer_det(new_size) do k=1,N_states do i=1,H_apply_buffer_N_det H_apply_buffer_coef(i,k) = buffer_coef(i,k) + H_apply_buffer_e2(i,k) = buffer_e2(i,k) enddo enddo - deallocate (buffer_det, buffer_coef) - SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det + deallocate (buffer_det, buffer_coef, buffer_e2) + SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det H_apply_buffer_e2 end BEGIN_PROVIDER [ integer(bit_kind), H_apply_buffer_det,(N_int,2,H_apply_buffer_size) ] &BEGIN_PROVIDER [ double precision, H_apply_buffer_coef,(H_apply_buffer_size,N_states) ] +&BEGIN_PROVIDER [ double precision, H_apply_buffer_e2,(H_apply_buffer_size,N_states) ] &BEGIN_PROVIDER [ integer, H_apply_buffer_N_det ] implicit none BEGIN_DOC - ! Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. + ! Buffer of determinants/coefficients/perturbative energy for H_apply. + ! Uninitialized. Filled by H_apply subroutines. END_DOC H_apply_buffer_N_det = 0 @@ -148,8 +153,9 @@ subroutine copy_H_apply_buffer_to_wf psi_coef(i+N_det_old,k) = H_apply_buffer_coef(i,k) enddo enddo + H_apply_buffer_N_det = 0 - SOFT_TOUCH psi_det psi_coef + SOFT_TOUCH psi_det psi_coef H_apply_buffer_N_det H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_e2 end diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 3c13d31f..14fea836 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -2,6 +2,12 @@ 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 + integer,parameter :: size_max = $size_max $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -21,21 +27,18 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) - 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 + PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref + PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref $set_i_H_j_threshold $omp_init_lock - E_ref = diag_H_mat_elem(key_in,N_int) - $initialization $omp_parallel @@ -153,8 +156,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame hij_tab(key_idx) = hij_elec ASSERT (key_idx <= size_max) if (key_idx == size_max) then + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -162,7 +166,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo if (key_idx > ishft(size_max,-5)) then if ($omp_test_lock) then - $keys_work + $keys_work_unlocked + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -201,8 +206,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame hij_tab(key_idx) = hij_elec ASSERT (key_idx <= size_max) if (key_idx == size_max) then + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock key_idx = 0 endif @@ -210,7 +216,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo if (key_idx > ishft(size_max,-5)) then if ($omp_test_lock) then - $keys_work + $keys_work_locked + $keys_work_unlocked $omp_unset_lock key_idx = 0 endif @@ -219,8 +226,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame enddo ! ii $omp_enddo enddo ! ispin + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock deallocate (keys_out,hij_tab,ia_ja_pairs) $omp_end_parallel @@ -233,6 +241,12 @@ 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 + integer,parameter :: size_max = $size_max $declarations integer(omp_lock_kind) :: lck integer(bit_kind),intent(in) :: key_in(N_int,2) @@ -252,21 +266,18 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) - 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 + PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref + PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref $set_i_H_j_threshold $omp_init_lock - E_ref = diag_H_mat_elem(key_in,N_int) - $initialization $omp_parallel @@ -305,7 +316,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) integer(bit_kind) :: test(N_int,2) double precision :: accu accu = 0.d0 - hij_elec = 0.d0 do ispin=1,2 other_spin = iand(ispin,1)+1 $omp_do @@ -319,33 +329,33 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters ) k_a = ishft(j_a-1,-bit_kind_shift)+1 l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) - call i_H_j(hole,key_in,N_int,hij_elec) - if(dabs(hij_elec) .ge. thresh)then - key_idx += 1 - do k=1,N_int - keys_out(k,1,key_idx) = hole(k,1) - keys_out(k,2,key_idx) = hole(k,2) - enddo - hij_tab(key_idx) = hij_elec - if (key_idx > ishft(size_max,-5)) then - if ($omp_test_lock) then - $keys_work - $omp_unset_lock - key_idx = 0 - endif - endif - if (key_idx == size_max) then - $omp_set_lock - $keys_work + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = hole(k,1) + keys_out(k,2,key_idx) = hole(k,2) + enddo + hij_tab(key_idx) = hij_elec + if (key_idx > ishft(size_max,-5)) then + if ($omp_test_lock) then + $keys_work_unlocked + $keys_work_locked $omp_unset_lock key_idx = 0 endif endif + if (key_idx == size_max) then + $keys_work_unlocked + $omp_set_lock + $keys_work_locked + $omp_unset_lock + key_idx = 0 + endif enddo ! ii $omp_enddo enddo ! ispin + $keys_work_unlocked $omp_set_lock - $keys_work + $keys_work_locked $omp_unset_lock deallocate (keys_out,hij_tab,ia_ja_pairs) $omp_end_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..1fdaacf5 --- /dev/null +++ b/src/Dets/connected_to_ref.irp.f @@ -0,0 +1,256 @@ +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(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,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 > 0) + 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) + use bitmasks + implicit none + BEGIN_DOC + ! If true, det is not in ref + ! If false, det may be in ref + END_DOC + integer(bit_kind), intent(in) :: key(Nint,2), Nint + integer(bit_kind) :: key_int + integer*1 :: key_short(bit_kind) + !DIR$ ATTRIBUTES ALIGN : 32 :: key_short + equivalence (key_int,key_short) + + integer :: i, ispin, k + + det_is_not_or_may_be_in_ref = .False. + do ispin=1,2 + do i=1,Nint + key_int = key(i,ispin) + do k=1,bit_kind + det_is_not_or_may_be_in_ref = & + det_is_not_or_may_be_in_ref .or. & + key_pattern_not_in_ref(key_short(k), i, ispin) + enddo + if(det_is_not_or_may_be_in_ref) then + return + endif + enddo + enddo + +end + + +BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ] + use bitmasks + implicit none + BEGIN_DOC + ! Min and max values of the integers of the keys of the reference + END_DOC + + integer :: i, j, ispin + integer(bit_kind) :: key + integer*1 :: key_short(bit_kind) + equivalence (key,key_short) + integer :: idx, k + + 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) + do k=1,bit_kind + key_pattern_not_in_ref( key_short(k), i, ispin ) = .False. + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index 780d54b5..09fad255 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 @@ -279,3 +275,4 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint) ) end + diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f new file mode 100644 index 00000000..0353c9c6 --- /dev/null +++ b/src/Dets/filter_connected.irp.f @@ -0,0 +1,175 @@ + +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 (Nint == N_int) + 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(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 + 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/s2.irp.f b/src/Dets/s2.irp.f index c970d628..b218e7cb 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -1,11 +1,12 @@ subroutine get_s2(key_i,key_j,phase,Nint) implicit none + use bitmasks BEGIN_DOC ! Returns END_DOC integer, intent(in) :: Nint - integer, intent(in) :: key_i(Nint,2) - integer, intent(in) :: key_j(Nint,2) + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer(bit_kind), intent(in) :: key_j(Nint,2) double precision, intent(out) :: phase integer :: exc(0:2,2,2) integer :: degree @@ -32,3 +33,38 @@ subroutine get_s2(key_i,key_j,phase,Nint) end select end +BEGIN_PROVIDER [ double precision, S_z ] +&BEGIN_PROVIDER [ double precision, S_z2_Sz ] + implicit none + + S_z = 0.5d0*dble(elec_alpha_num-elec_beta_num) + S_z2_Sz = S_z*(S_z-1.d0) + +END_PROVIDER + + +subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) + integer, intent(in) :: n,nmax + double precision, intent(in) :: psi_coefs_tmp(nmax) + double precision, intent(out) :: s2 + integer :: i,j,l + double precision :: s2_tmp + s2 = S_z2_Sz + !$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 + 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 + enddo + !$OMP END PARALLEL DO +end + diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 59dc68ad..3f56eba3 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -488,13 +488,13 @@ 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 - integer, intent(in) :: keys(Nint,2,Ndet_max) + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) double precision, intent(in) :: coef(Ndet_max,Nstate) - integer, intent(in) :: key(Nint,2) double precision, intent(out) :: i_H_psi_array(Nstate) integer :: i, ii,j @@ -503,6 +503,11 @@ subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) double precision :: hij integer :: idx(0:Ndet) + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) do ii=1,idx(0) @@ -512,6 +517,7 @@ subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij enddo + print *, 'x', coef(i,1), hij, i_H_psi_array(1) enddo end @@ -600,180 +606,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 +795,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 +804,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/MOGuess/README.rst b/src/MOGuess/README.rst index d23ebe86..d8e72641 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -23,7 +23,8 @@ Documentation .. NEEDED_MODULES file. `h_core_guess `_ -None + Undocumented + `ao_ortho_lowdin_coef `_ matrix of the coefficients of the mos generated by the orthonormalization by the S^{-1/2} canonical transformation of the aos @@ -34,6 +35,7 @@ None supposed to be the Identity `ao_ortho_lowdin_nucl_elec_integral `_ -None + Undocumented + diff --git a/src/Makefile.config.example b/src/Makefile.config.example deleted file mode 100644 index 82e0ca27..00000000 --- a/src/Makefile.config.example +++ /dev/null @@ -1,29 +0,0 @@ -OPENMP =1 -PROFILE =0 -DEBUG = 0 - -IRPF90_FLAGS+= --align=32 -FC = ifort -g -FCFLAGS= -FCFLAGS+= -xHost -#FCFLAGS+= -xAVX -FCFLAGS+= -O2 -FCFLAGS+= -ip -FCFLAGS+= -opt-prefetch -FCFLAGS+= -ftz -MKL=-mkl=parallel - -ifeq ($(PROFILE),1) -FC += -p -g -CXX += -pg -endif - -ifeq ($(OPENMP),1) -FC += -openmp -CXX += -fopenmp -endif - -ifeq ($(DEBUG),1) -FC += -C -traceback -fpe0 -#FCFLAGS =-O0 -endif diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 03ec9939..052e2053 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -42,13 +42,17 @@ Documentation :math:`\int \chi_i(r) \chi_j(r) dr)` `check_ortho `_ -None + Undocumented + `do_print `_ -None + Undocumented + `n_pt_max_i_x `_ -None + Undocumented + `n_pt_max_integrals `_ -None + Undocumented + `ao_deriv2_x `_ second derivatives matrix elements in the ao basis .. math:: @@ -67,46 +71,169 @@ None .br {\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle -`ao_kinetic_integral `_ +`ao_kinetic_integral `_ array of the priminitve basis kinetic integrals \langle \chi_i |\hat{T}| \chi_j \rangle `mo_kinetic_integral `_ -None + Undocumented + `mo_mono_elec_integral `_ array of the mono electronic hamiltonian on the MOs basis : sum of the kinetic and nuclear electronic potential `mo_overlap `_ -None + Undocumented + `orthonormalize_mos `_ -None + Undocumented + `ao_nucl_elec_integral `_ interaction nuclear electron `give_polynom_mult_center_mono_elec `_ -None + Undocumented + `i_x1_pol_mult_mono_elec `_ -None + Undocumented + `i_x2_pol_mult_mono_elec `_ -None + Undocumented + `int_gaus_pol `_ -None + Undocumented + `nai_pol_mult `_ -None + Undocumented + `v_e_n `_ -None + Undocumented + `v_phi `_ -None + Undocumented + `v_r `_ -None + Undocumented + `v_theta `_ -None + Undocumented + `wallis `_ -None + Undocumented + `mo_nucl_elec_integral `_ -None + Undocumented + `save_ortho_mos `_ -None + Undocumented + +`ao_deriv_1_x `_ + array of the integrals of AO_i * d/dx AO_j + array of the integrals of AO_i * d/dy AO_j + array of the integrals of AO_i * d/dz AO_j + +`ao_deriv_1_y `_ + array of the integrals of AO_i * d/dx AO_j + array of the integrals of AO_i * d/dy AO_j + array of the integrals of AO_i * d/dz AO_j + +`ao_deriv_1_z `_ + array of the integrals of AO_i * d/dx AO_j + array of the integrals of AO_i * d/dy AO_j + array of the integrals of AO_i * d/dz AO_j + +`ao_dipole_x `_ + array of the integrals of AO_i * x AO_j + array of the integrals of AO_i * y AO_j + array of the integrals of AO_i * z AO_j + +`ao_dipole_y `_ + array of the integrals of AO_i * x AO_j + array of the integrals of AO_i * y AO_j + array of the integrals of AO_i * z AO_j + +`ao_dipole_z `_ + array of the integrals of AO_i * x AO_j + array of the integrals of AO_i * y AO_j + array of the integrals of AO_i * z AO_j + +`ao_spread_x `_ + array of the integrals of AO_i * x^2 AO_j + array of the integrals of AO_i * y^2 AO_j + array of the integrals of AO_i * z^2 AO_j + +`ao_spread_y `_ + array of the integrals of AO_i * x^2 AO_j + array of the integrals of AO_i * y^2 AO_j + array of the integrals of AO_i * z^2 AO_j + +`ao_spread_z `_ + array of the integrals of AO_i * x^2 AO_j + array of the integrals of AO_i * y^2 AO_j + array of the integrals of AO_i * z^2 AO_j + +`overlap_bourrin_deriv_x `_ + Undocumented + +`overlap_bourrin_dipole `_ + Undocumented + +`overlap_bourrin_spread `_ + Undocumented + +`overlap_bourrin_x `_ + Undocumented + +`overlap_bourrin_x_abs `_ + Undocumented + +`power `_ + Undocumented + +`mo_deriv_1_x `_ + array of the integrals of MO_i * d/dx MO_j + array of the integrals of MO_i * d/dy MO_j + array of the integrals of MO_i * d/dz MO_j + +`mo_deriv_1_y `_ + array of the integrals of MO_i * d/dx MO_j + array of the integrals of MO_i * d/dy MO_j + array of the integrals of MO_i * d/dz MO_j + +`mo_deriv_1_z `_ + array of the integrals of MO_i * d/dx MO_j + array of the integrals of MO_i * d/dy MO_j + array of the integrals of MO_i * d/dz MO_j + +`mo_dipole_x `_ + array of the integrals of MO_i * x MO_j + array of the integrals of MO_i * y MO_j + array of the integrals of MO_i * z MO_j + +`mo_dipole_y `_ + array of the integrals of MO_i * x MO_j + array of the integrals of MO_i * y MO_j + array of the integrals of MO_i * z MO_j + +`mo_dipole_z `_ + array of the integrals of MO_i * x MO_j + array of the integrals of MO_i * y MO_j + array of the integrals of MO_i * z MO_j + +`mo_spread_x `_ + array of the integrals of MO_i * x^2 MO_j + array of the integrals of MO_i * y^2 MO_j + array of the integrals of MO_i * z^2 MO_j + +`mo_spread_y `_ + array of the integrals of MO_i * x^2 MO_j + array of the integrals of MO_i * y^2 MO_j + array of the integrals of MO_i * z^2 MO_j + +`mo_spread_z `_ + array of the integrals of MO_i * x^2 MO_j + array of the integrals of MO_i * y^2 MO_j + array of the integrals of MO_i * z^2 MO_j + diff --git a/src/MonoInts/spread_dipole_ao.irp.f b/src/MonoInts/spread_dipole_ao.irp.f new file mode 100644 index 00000000..2628c6a0 --- /dev/null +++ b/src/MonoInts/spread_dipole_ao.irp.f @@ -0,0 +1,416 @@ + BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * x^2 AO_j + ! array of the integrals of AO_i * y^2 AO_j + ! array of the integrals of AO_i * z^2 AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + dim1=500 + lower_exp_val = 40.d0 + ao_spread_x= 0.d0 + ao_spread_y= 0.d0 + ao_spread_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + c = ao_coef_transp(n,j)*ao_coef_transp(l,i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x += c*(tmp*overlap_y*overlap_z) + call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y += c*(tmp*overlap_x*overlap_z) + call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z += c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_spread_x(i,j) = accu_x + ao_spread_y(i,j) = accu_y + ao_spread_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * x AO_j + ! array of the integrals of AO_i * y AO_j + ! array of the integrals of AO_i * z AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c + dim1=500 + lower_exp_val = 40.d0 + ao_dipole_x= 0.d0 + ao_dipole_y= 0.d0 + ao_dipole_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + c = ao_coef_transp(l,i)*ao_coef_transp(n,j) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x = accu_x + c*(tmp*overlap_y*overlap_z) + call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y = accu_y + c*(tmp*overlap_x*overlap_z) + call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z = accu_z + c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_dipole_x(i,j) = accu_x + ao_dipole_y(i,j) = accu_y + ao_dipole_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)] + &BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)] + BEGIN_DOC + ! array of the integrals of AO_i * d/dx AO_j + ! array of the integrals of AO_i * d/dy AO_j + ! array of the integrals of AO_i * d/dz AO_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + integer :: i_component + dim1=500 + lower_exp_val = 40.d0 + ao_deriv_1_x= 0.d0 + ao_deriv_1_y= 0.d0 + ao_deriv_1_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_power,ao_prim_num, & + !$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_transp,ao_nucl, & + !$OMP ao_expo_transp,dim1,lower_exp_val) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + !DEC$ VECTOR ALIGNED + !DEC$ VECTOR ALWAYS + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + !DEC$ VECTOR ALIGNED + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + c = ao_coef_transp(l,i) * ao_coef_transp(n,j) + i_component = 1 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_x += c*(tmp*overlap_y*overlap_z) + i_component = 2 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_y += c*(tmp*overlap_x*overlap_z) + i_component = 3 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_z += c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_deriv_1_x(i,j) = accu_x + ao_deriv_1_y(i,j) = accu_y + ao_deriv_1_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + subroutine overlap_bourrin_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy +!print*,'xmin = ',x_min +!print*,'xmax = ',x_max + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += (power(power_A,x-A_center) * power(power_B,x-B_center)) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + overlap_x *= factor * dx +end + + subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! needed for the dipole and those things + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then +! print*,'factor = ',factor + dx = 0.d0 + overlap_x = 0.d0 + return + endif + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x * x + enddo + overlap_x *= factor * dx + + end + + double precision function power(n,x) + implicit none + integer :: i,n + double precision :: x,accu + power = x**n + return + end + + subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! needed for the dipole and those things + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + double precision :: pouet_timy + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x + enddo + overlap_x *= factor * dx + + end + + subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx) + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A(3),power_B(3),i_component + double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: overlap_first, overlap_second +! computes : = (a_x_i - 2 alpha ) + + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx) + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx) + overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second) + end + + subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then + dx = 0.d0 + overlap_x = 0.d0 + return + endif + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + overlap_x *= factor * dx + end + diff --git a/src/MonoInts/spread_dipole_mo.irp.f b/src/MonoInts/spread_dipole_mo.irp.f new file mode 100644 index 00000000..d7306727 --- /dev/null +++ b/src/MonoInts/spread_dipole_mo.irp.f @@ -0,0 +1,101 @@ + BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * x MO_j + ! array of the integrals of MO_i * y MO_j + ! array of the integrals of MO_i * z MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_dipole_x = 0.d0 + mo_dipole_y = 0.d0 + mo_dipole_z = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_dipole_x,mo_dipole_y,mo_dipole_z,ao_dipole_x,ao_dipole_y,ao_dipole_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_dipole_x(j,i) = mo_dipole_x(j,i) + c_j1 * ao_dipole_x(j1,i1) + mo_dipole_y(j,i) = mo_dipole_y(j,i) + c_j1 * ao_dipole_y(j1,i1) + mo_dipole_z(j,i) = mo_dipole_z(j,i) + c_j1 * ao_dipole_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + + BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * x^2 MO_j + ! array of the integrals of MO_i * y^2 MO_j + ! array of the integrals of MO_i * z^2 MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_nucl_elec_integral = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_spread_x,mo_spread_y,mo_spread_z,ao_spread_x,ao_spread_y,ao_spread_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_spread_x(j,i) = mo_spread_x(j,i) + c_j1 * ao_spread_x(j1,i1) + mo_spread_y(j,i) = mo_spread_y(j,i) + c_j1 * ao_spread_y(j1,i1) + mo_spread_z(j,i) = mo_spread_z(j,i) + c_j1 * ao_spread_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + + BEGIN_PROVIDER [double precision, mo_deriv_1_x , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_deriv_1_y , (mo_tot_num_align,mo_tot_num)] +&BEGIN_PROVIDER [double precision, mo_deriv_1_z , (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC + ! array of the integrals of MO_i * d/dx MO_j + ! array of the integrals of MO_i * d/dy MO_j + ! array of the integrals of MO_i * d/dz MO_j + END_DOC + implicit none + integer :: i1,j1,i,j + double precision :: c_i1,c_j1 + + mo_nucl_elec_integral = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_deriv_1_x,mo_deriv_1_y,mo_deriv_1_z,ao_spread_x,ao_spread_y,ao_spread_z) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_deriv_1_x(j,i) = mo_deriv_1_x(j,i) + c_j1 * ao_spread_x(j1,i1) + mo_deriv_1_y(j,i) = mo_deriv_1_y(j,i) + c_j1 * ao_spread_y(j1,i1) + mo_deriv_1_z(j,i) = mo_deriv_1_z(j,i) + c_j1 * ao_spread_z(j1,i1) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 1a75cb06..d09721ef 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 Selection diff --git a/src/Perturbation/ASSUMPTIONS.rst b/src/Perturbation/ASSUMPTIONS.rst new file mode 100644 index 00000000..fd22e8a3 --- /dev/null +++ b/src/Perturbation/ASSUMPTIONS.rst @@ -0,0 +1,6 @@ +* This is not allowed: + + subroutine & + pt2_.... + + diff --git a/src/Perturbation/Makefile b/src/Perturbation/Makefile new file mode 100644 index 00000000..ea95cbf1 --- /dev/null +++ b/src/Perturbation/Makefile @@ -0,0 +1,8 @@ +default: all + +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC=perturbation_template.f +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Perturbation/NEEDED_MODULES b/src/Perturbation/NEEDED_MODULES new file mode 100644 index 00000000..f6551dd9 --- /dev/null +++ b/src/Perturbation/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst new file mode 100644 index 00000000..83392b83 --- /dev/null +++ b/src/Perturbation/README.rst @@ -0,0 +1,125 @@ +=================== +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. + +`pt2_epstein_nesbet `_ + compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states. + .br + c_pert(i) = /( E(i) - ) + .br + e_2_pert(i) = ^2/( E(i) - ) + .br + +`pt2_epstein_nesbet_2x2 `_ + compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + .br + for the various n_st states. + .br + e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + .br + c_pert(i) = e_2_pert(i)/ + .br + + + +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/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f new file mode 100644 index 00000000..4c1e7c45 --- /dev/null +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -0,0 +1,68 @@ +subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,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 + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states. + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem + ASSERT (Nint == N_int) + 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) + do i =1,n_st + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end + +subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,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 + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + ! + ! for the various n_st states. + ! + ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + ! + ! c_pert(i) = e_2_pert(i)/ + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,delta_e + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + print *, 'coefs' + print *, psi_ref_coef(1:N_det_ref,1) + print *, '-----' + call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,N_det_ref,psi_ref_size,n_st,i_H_psi_array) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + do i =1,n_st + delta_e = H_pert_diag - reference_energy(i) + e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + enddo + print *, e_2_pert, delta_e , i_H_psi_array + +end 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 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f new file mode 100644 index 00000000..3365e67c --- /dev/null +++ b/src/Perturbation/perturbation_template.f @@ -0,0 +1,47 @@ +BEGIN_SHELL [ /usr/bin/env python ] +import perturbation +END_SHELL + +subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint) + implicit none + BEGIN_DOC + ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply + ! routine. + END_DOC + + integer, intent(in) :: Nint, N_st, buffer_size + integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) + double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) + double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) + double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag + integer :: i,k, c_ref + integer :: connected_to_ref + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (buffer_size >= 0) + ASSERT (minval(sum_norm_pert) >= 0.d0) + ASSERT (N_st > 0) + do i = 1,buffer_size + + c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_ref,N_det,h_apply_threshold) + + if (c_ref /= 0) then + cycle + endif + + call pt2_$PERT(buffer(1,1,i), & + c_pert,e_2_pert,H_pert_diag,Nint,n_det_ref,n_st) + + do k = 1,N_st + e_2_pert_buffer(k,i) = e_2_pert(k) + coef_pert_buffer(k,i) = c_pert(k) + sum_norm_pert(k) += c_pert(k) * c_pert(k) + sum_e_2_pert(k) += e_2_pert(k) + sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag + enddo + + enddo + +end + diff --git a/src/README.rst b/src/README.rst index 9c6b5bcb..3aae2807 100644 --- a/src/README.rst +++ b/src/README.rst @@ -129,4 +129,13 @@ Needed Modules * `Dets `_ * `DensityMatrix `_ * `CISD `_ +* `Perturbation `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index ecae8cdc..22ca38a8 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -76,6 +76,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP END DO NOWAIT enddo + !$OMP BARRIER !$OMP DO do j=1,n do i=1,m @@ -181,6 +182,13 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! END_DOC integer, intent(in) :: n,nmax double precision, intent(out) :: eigvectors(nmax,n) @@ -189,7 +197,6 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) double precision,allocatable :: eigenvalues(:) double precision,allocatable :: work(:) double precision,allocatable :: A(:,:) - !eigvectors(i,j) = where d_i is the basis function and psi_j is the j th eigenvector allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax)) integer :: LWORK, info, i,j,l,k A=H