diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index c78a3826..7c1f43c2 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -386,39 +386,52 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun ! ============== - 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 + if (N_st > 1) then - call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) + 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) + + + ! 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) + + else + + overlap(1,1) = u_dot_u(U_in(1,1),sze) + double precision :: f + f = 1.d0 / dsqrt(overlap(1,1)) + do i=1,sze + U_in(i,1) = U_in(i,1) * f + enddo + + endif ! Davidson iterations ! =================== diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 13138499..11de1270 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -73,6 +73,10 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j + if (n < 2) then + return + endif + allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) @@ -151,6 +155,10 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j, k + if (n < 2) then + return + endif + call svd(overlap,lda,U,ldc,D,Vt,lda,m,n) !$OMP PARALLEL DEFAULT(NONE) &