From 82772b96c7675b263e1b0a0b13e3678d38f17ae9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 29 Dec 2016 22:00:41 +0100 Subject: [PATCH] MRCC_selected --- plugins/Psiref_threshold/psi_ref.irp.f | 71 ++++-- plugins/mrcc_selected/dressing.irp.f | 80 ++----- plugins/mrcc_selected/mrcc_selected.irp.f | 1 - plugins/mrcc_selected/mrcepa0_general.irp.f | 15 +- src/Davidson/u0Hu0.irp.f | 234 ++++++++++++++------ src/Determinants/Fock_diag.irp.f | 9 + 6 files changed, 239 insertions(+), 171 deletions(-) diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f index ee69ef5c..62321140 100644 --- a/plugins/Psiref_threshold/psi_ref.irp.f +++ b/plugins/Psiref_threshold/psi_ref.irp.f @@ -1,5 +1,44 @@ use bitmasks +! BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] +!&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] +!&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] +!&BEGIN_PROVIDER [ integer, N_det_ref ] +! implicit none +! BEGIN_DOC +! ! Reference wave function, defined as determinants with amplitudes > 0.05 +! ! idx_ref gives the indice of the ref determinant in psi_det. +! END_DOC +! integer :: i, k, l +! logical :: good +! double precision, parameter :: threshold=0.01d0 +! double precision :: t(N_states) +! N_det_ref = 0 +! do l = 1, N_states +! t(l) = threshold * abs_psi_coef_max(l) +! enddo +! do i=1,N_det +! good = .False. +! do l=1, N_states +! psi_ref_coef(i,l) = 0.d0 +! good = good.or.(dabs(psi_coef(i,l)) > t(l)) +! enddo +! if (good) then +! N_det_ref = N_det_ref+1 +! do k=1,N_int +! psi_ref(k,1,N_det_ref) = psi_det(k,1,i) +! psi_ref(k,2,N_det_ref) = psi_det(k,2,i) +! enddo +! idx_ref(N_det_ref) = i +! do k=1,N_states +! psi_ref_coef(N_det_ref,k) = psi_coef(i,k) +! enddo +! endif +! enddo +! call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') +! +!END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] &BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] @@ -10,30 +49,16 @@ use bitmasks ! idx_ref gives the indice of the ref determinant in psi_det. END_DOC integer :: i, k, l - logical :: good - double precision, parameter :: threshold=0.05d0 - double precision :: t(N_states) - N_det_ref = 0 - do l = 1, N_states - t(l) = threshold * abs_psi_coef_max(l) - enddo - do i=1,N_det - good = .False. - do l=1, N_states - psi_ref_coef(i,l) = 0.d0 - good = good.or.(dabs(psi_coef(i,l)) > t(l)) + double precision, parameter :: threshold=0.01d0 + + call find_reference(threshold, N_det_ref, idx_ref) + do l=1,N_states + do i=1,N_det_ref + psi_ref_coef(i,l) = psi_coef(idx_ref(i), l) enddo - if (good) then - N_det_ref = N_det_ref+1 - do k=1,N_int - psi_ref(k,1,N_det_ref) = psi_det(k,1,i) - psi_ref(k,2,N_det_ref) = psi_det(k,2,i) - enddo - idx_ref(N_det_ref) = i - do k=1,N_states - psi_ref_coef(N_det_ref,k) = psi_coef(i,k) - enddo - endif + enddo + do i=1,N_det_ref + psi_ref(:,:,i) = psi_det(:,:,idx_ref(i)) enddo call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') diff --git a/plugins/mrcc_selected/dressing.irp.f b/plugins/mrcc_selected/dressing.irp.f index c772e2aa..23fedcee 100644 --- a/plugins/mrcc_selected/dressing.irp.f +++ b/plugins/mrcc_selected/dressing.irp.f @@ -534,63 +534,9 @@ END_PROVIDER END_PROVIDER -! BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] -! use bitmasks -! implicit none -! integer :: i,j,k -! double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall -! integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref) -! -! ! provide lambda_mrcc -! npres = 0 -! delta_cas = 0d0 -! call wall_time(wall) -! print *, "dcas ", wall -! do i_state = 1, N_states -! !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref) -! do k=1,N_det_non_ref -! if(lambda_mrcc(i_state, k) == 0d0) cycle -! npre = 0 -! do i=1,N_det_ref -! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) -! if(Hki /= 0d0) then -! !!$OMP ATOMIC -! npres(i) += 1 -! npre += 1 -! ipre(npre) = i -! pre(npre) = Hki -! end if -! end do -! -! -! do i=1,npre -! do j=1,i -! !!$OMP ATOMIC -! delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k) -! end do -! end do -! end do -! !!$OMP END PARALLEL DO -! npre=0 -! do i=1,N_det_ref -! npre += npres(i) -! end do -! !stop -! do i=1,N_det_ref -! do j=1,i -! delta_cas(j,i,i_state) = delta_cas(i,j,i_state) -! end do -! end do -! end do -! -! call wall_time(wall) -! print *, "dcas", wall -! ! stop -! END_PROVIDER - - BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] -&BEGIN_PROVIDER [ double precision, delta_cas_s2, (N_det_ref, N_det_ref, N_states) ] + BEGIN_PROVIDER [ double precision, delta_ref, (N_det_ref, N_det_ref, N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ref_s2, (N_det_ref, N_det_ref, N_states) ] use bitmasks implicit none integer :: i,j,k @@ -600,22 +546,22 @@ END_PROVIDER provide lambda_mrcc dIj do i_state = 1, N_states - !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,delta_cas_s2,N_det_ref,dij) + !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_ref,delta_ref_s2,N_det_ref,dij) do i=1,N_det_ref do j=1,i call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int) - delta_cas(i,j,i_state) = 0d0 - delta_cas_s2(i,j,i_state) = 0d0 + delta_ref(i,j,i_state) = 0d0 + delta_ref_s2(i,j,i_state) = 0d0 do k=1,N_det_non_ref call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) call get_s2(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Sjk) - delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) - delta_cas_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k) + delta_ref(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) + delta_ref_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k) end do - delta_cas(j,i,i_state) = delta_cas(i,j,i_state) - delta_cas_s2(j,i,i_state) = delta_cas_s2(i,j,i_state) + delta_ref(j,i,i_state) = delta_ref(i,j,i_state) + delta_ref_s2(j,i,i_state) = delta_ref_s2(i,j,i_state) end do end do !$OMP END PARALLEL DO @@ -739,7 +685,7 @@ end function !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) & !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) & !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & - !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas, delta_cas_s2) & + !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_ref, delta_ref_s2) & !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) do blok=1,cepa0_shortcut(0) do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 @@ -781,8 +727,8 @@ end function notf = notf+1 ! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) - contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) - contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + contrib = delta_ref(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + contrib_s2 = delta_ref_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) @@ -828,7 +774,7 @@ END_PROVIDER integer :: II, blok - provide delta_cas lambda_mrcc + provide delta_ref lambda_mrcc allocate(idx_sorted_bit(N_det)) idx_sorted_bit(:) = -1 do i=1,N_det_non_ref diff --git a/plugins/mrcc_selected/mrcc_selected.irp.f b/plugins/mrcc_selected/mrcc_selected.irp.f index 91592e62..b64f968d 100644 --- a/plugins/mrcc_selected/mrcc_selected.irp.f +++ b/plugins/mrcc_selected/mrcc_selected.irp.f @@ -8,7 +8,6 @@ program mrsc2sub read_wf = .True. SOFT_TOUCH read_wf - call print_cas_coefs call set_generators_bitmasks_as_holes_and_particles call run(N_states,energy) if(do_pt2_end)then diff --git a/plugins/mrcc_selected/mrcepa0_general.irp.f b/plugins/mrcc_selected/mrcepa0_general.irp.f index e3a2d1f5..812aeef0 100644 --- a/plugins/mrcc_selected/mrcepa0_general.irp.f +++ b/plugins/mrcc_selected/mrcepa0_general.irp.f @@ -60,16 +60,17 @@ subroutine run(N_st,energy) end -subroutine print_cas_coefs +subroutine print_ref_coefs implicit none integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, (psi_cas_coef(i,j), j=1,N_states) - call debug_det(psi_cas(1,1,i),N_int) + print *, 'Reference' + print *, '=========' + do i=1,N_det_ref + print *, (psi_ref_coef(i,j), j=1,N_states) + call debug_det(psi_ref(1,1,i),N_int) enddo + print *, '' call write_double(6,ci_energy(1),"Initial CI energy") end @@ -202,7 +203,7 @@ subroutine run_pt2(N_st,energy) print*,'Last iteration only to compute the PT2' - N_det_generators = N_det_cas + N_det_generators = N_det_ref N_det_selectors = N_det_non_ref do i=1,N_det_generators diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 6e20f0d0..9e76bc92 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -32,20 +32,20 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> + ! Computes v_0 = H|u_0> ! ! n : number of determinants ! ! H_jj : array of + ! END_DOC integer, intent(in) :: N_st,n,Nint, sze_8 double precision, intent(out) :: v_0(sze_8,N_st) double precision, intent(in) :: u_0(sze_8,N_st) double precision, intent(in) :: H_jj(n) integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - double precision :: hij - double precision, allocatable :: vt(:,:) - double precision, allocatable :: ut(:,:) + double precision :: hij,s2 + double precision, allocatable :: vt(:,:), ut(:,:), st(:,:) integer :: i,j,k,l, jj,ii integer :: i0, j0 @@ -57,77 +57,41 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) integer :: N_st_8 integer, external :: align_double - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut + integer :: blockb, blockb2, istep + double precision :: ave_workload, workload, target_workload_inv + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut, st N_st_8 = align_double(N_st) ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) - PROVIDE ref_bitmask_energy + PROVIDE ref_bitmask_energy allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - allocate(ut(N_st_8,n)) + allocate( ut(N_st_8,n)) v_0 = 0.d0 - do i=1,n - do istate=1,N_st - ut(istate,i) = u_0(i,istate) - enddo - enddo - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) - + !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,H_jj,keys_tmp,ut,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) - allocate(vt(N_st_8,n)) + !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& + !$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) + allocate(vt(N_st_8,n),st(N_st_8,n)) Vt = 0.d0 - - !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0,1) - do sh2=1,shortcut(0,1) - exa = popcnt(xor(version(1,sh,1), version(1,sh2,1))) - if(exa > 2) then - cycle - end if - do ni=2,Nint - exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) - end do - if(exa > 2) then - cycle - end if - - do i=shortcut(sh,1),shortcut(sh+1,1)-1 - org_i = sort_idx(i,1) - do ni=1,Nint - sorted_i(ni) = sorted(ni,i,1) - enddo - - jloop: do j=shortcut(sh2,1),shortcut(sh2+1,1)-1 - org_j = sort_idx(j,1) - ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) - if(ext > 4) then - cycle jloop - endif - do ni=2,Nint - ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) - if(ext > 4) then - cycle jloop - endif - end do - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - enddo - enddo jloop - enddo + St = 0.d0 + + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,2),istate) enddo enddo - !$OMP END DO NOWAIT - + !$OMP END DO + !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0,2) do i=shortcut(sh,2),shortcut(sh+1,2)-1 @@ -135,40 +99,164 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) do j=shortcut(sh,2),shortcut(sh+1,2)-1 org_j = sort_idx(j,2) ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2))) + if (ext > 4) cycle do ni=2,Nint ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) + if (ext > 4) exit end do - if(ext /= 4) then - cycle - endif - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - enddo + if(ext == 4) then + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + end if end do end do enddo - !$OMP END DO NOWAIT - - !$OMP CRITICAL + !$OMP END DO + + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,1),istate) + enddo + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,1) + do sh2=1,shortcut(0,1) + if (sh==sh2) cycle + + exa = 0 + do ni=1,Nint + exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) + end do + if(exa > 2) then + cycle + end if + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh2,1),shortcut(sh2+1,1)-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif + enddo + + enddo + enddo + + exa = 0 + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh,1),i-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif + enddo + + do j=i+1,shortcut(sh+1,1)-1 + if (i==j) cycle + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif + enddo + enddo + enddo + !$OMP END DO + + !$OMP CRITICAL (u0Hu0) do istate=1,N_st - do i=n,1,-1 + do i=1,n v_0(i,istate) = v_0(i,istate) + vt(istate,i) enddo enddo - !$OMP END CRITICAL + !$OMP END CRITICAL (u0Hu0) - deallocate(vt) + deallocate(vt,st) !$OMP END PARALLEL - + do istate=1,N_st do i=1,n - v_0(i,istate) += H_jj(i) * u_0(i,istate) + v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate) enddo enddo deallocate (shortcut, sort_idx, sorted, version, ut) end + BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ] implicit none BEGIN_DOC diff --git a/src/Determinants/Fock_diag.irp.f b/src/Determinants/Fock_diag.irp.f index a99bbcad..01393fe1 100644 --- a/src/Determinants/Fock_diag.irp.f +++ b/src/Determinants/Fock_diag.irp.f @@ -19,6 +19,15 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint) fock_diag_tmp = 0.d0 E0 = 0.d0 + if (Ne(1) /= elec_alpha_num) then + print *, 'Error in build_fock_tmp (alpha)', Ne(1), Ne(2) + stop -1 + endif + if (Ne(2) /= elec_beta_num) then + print *, 'Error in build_fock_tmp (beta)', Ne(1), Ne(2) + stop -1 + endif + ! Occupied MOs do ii=1,elec_alpha_num i = occ(ii,1)