From a5ff15f4592d4a36d2140b8c73cd3b0bc3057fcf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 13 Sep 2017 16:50:45 +0200 Subject: [PATCH] Introduced Phi_S --- plugins/analyze_wf/attachment.irp.f | 64 +++++++++++ plugins/analyze_wf/dump_nto.irp.f | 19 +++ plugins/analyze_wf/dump_one_body_mos.irp.f | 38 ++++++ plugins/analyze_wf/phi_s.irp.f | 128 +++++++++++++++++++++ src/AO_Basis/ao_overlap.irp.f | 38 ++++++ src/Determinants/density_matrix.irp.f | 26 +++++ src/Utils/LinearAlgebra.irp.f | 26 ++--- 7 files changed, 326 insertions(+), 13 deletions(-) create mode 100644 plugins/analyze_wf/attachment.irp.f create mode 100644 plugins/analyze_wf/dump_nto.irp.f create mode 100644 plugins/analyze_wf/dump_one_body_mos.irp.f create mode 100644 plugins/analyze_wf/phi_s.irp.f diff --git a/plugins/analyze_wf/attachment.irp.f b/plugins/analyze_wf/attachment.irp.f new file mode 100644 index 00000000..d086aa21 --- /dev/null +++ b/plugins/analyze_wf/attachment.irp.f @@ -0,0 +1,64 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_detachment, (mo_tot_num,mo_tot_num,2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_attachment, (mo_tot_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Detachment and attachment density matrices in MO basis + END_DOC + integer :: i,j, k, istate + double precision :: km(mo_tot_num), kp(mo_tot_num) + + one_body_dm_mo_detachment = 0.d0 + one_body_dm_mo_attachment = 0.d0 + + do istate=2,N_states + + km(:) = 0.d0 + kp(:) = 0.d0 + do i=1,mo_tot_num + if (one_body_dm_mo_diff_eigvalues(i,istate) < 0) then + km(i) = -one_body_dm_mo_diff_eigvalues(i,istate) + else + kp(i) = one_body_dm_mo_diff_eigvalues(i,istate) + endif + enddo + + ! Attachment + do k=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + one_body_dm_mo_detachment(i,j,istate) = one_body_dm_mo_detachment(i,j,istate) + & + one_body_dm_mo_diff_eigvectors(i,k,istate)*km(k)* & + one_body_dm_mo_diff_eigvectors(j,k,istate) + one_body_dm_mo_attachment(i,j,istate) = one_body_dm_mo_attachment(i,j,istate) + & + one_body_dm_mo_diff_eigvectors(i,k,istate)*kp(k)* & + one_body_dm_mo_diff_eigvectors(j,k,istate) + enddo + enddo + enddo + + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_ao_detachment, (ao_num,ao_num,2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_ao_attachment, (ao_num,ao_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Detachment and attachment density matrices in AO basis + END_DOC + integer :: istate + do istate=2,N_states + call mo_to_ao_no_overlap( & + one_body_dm_mo_attachment(1,1,istate), & + size(one_body_dm_mo_attachment,1), & + one_body_dm_ao_attachment(1,1,istate), & + size(one_body_dm_ao_attachment,1) ) + call mo_to_ao_no_overlap( & + one_body_dm_mo_detachment(1,1,istate), & + size(one_body_dm_mo_detachment,1), & + one_body_dm_ao_detachment(1,1,istate), & + size(one_body_dm_ao_detachment,1) ) + enddo + +END_PROVIDER + diff --git a/plugins/analyze_wf/dump_nto.irp.f b/plugins/analyze_wf/dump_nto.irp.f new file mode 100644 index 00000000..8d19c3eb --- /dev/null +++ b/plugins/analyze_wf/dump_nto.irp.f @@ -0,0 +1,19 @@ +program dump_nto + implicit none + integer :: i,j, istate + + print *, 'Phi_S' + do i=2,N_states + print *, i, Phi_S(i) + enddo + + do istate=2,N_states + do j=1,mo_tot_num + print *, 'MO: ', j, 'State:', istate, 'Eig:', one_body_dm_mo_diff_eigvalues(j,istate) + do i=1,ao_num + print *, i, transition_natorb(i,j,istate) + enddo + enddo + enddo + +end diff --git a/plugins/analyze_wf/dump_one_body_mos.irp.f b/plugins/analyze_wf/dump_one_body_mos.irp.f new file mode 100644 index 00000000..7ab841ef --- /dev/null +++ b/plugins/analyze_wf/dump_one_body_mos.irp.f @@ -0,0 +1,38 @@ +program dump_one_body_mos + implicit none + BEGIN_DOC +! Output density matrices of all the states + END_DOC + read_wf = .True. + TOUCH read_wf + call run() +end + +subroutine run + implicit none + integer :: istate + integer, parameter :: iunit = 66 + character*(64) :: filename, fmt + integer :: i,j,k + + write(fmt,'(''('',I4.4,''(X,E20.14))'')') mo_tot_num + do istate=1,N_states + write(filename,'(''state.'',I2.2)') istate + open(unit=iunit, form='formatted', file=filename) + write(iunit,*) mo_tot_num + do j=1,mo_tot_num + write(iunit,fmt) one_body_dm_mo_alpha(1:mo_tot_num,j,istate) + one_body_dm_mo_beta(1:mo_tot_num,j,istate) + enddo + enddo + + + call run2() +end +subroutine run2 + integer :: i,j, istate + print *, 'Phi_S' + do i=2,N_states + print *, i, Phi_S(i) + enddo +end + diff --git a/plugins/analyze_wf/phi_s.irp.f b/plugins/analyze_wf/phi_s.irp.f new file mode 100644 index 00000000..12bdb970 --- /dev/null +++ b/plugins/analyze_wf/phi_s.irp.f @@ -0,0 +1,128 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvalues, (mo_tot_num, 2:N_states) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff_eigvectors, (mo_tot_num, mo_tot_num, 2:N_states) ] + implicit none + BEGIN_DOC + ! Eigenvalues and eigenvectors of one_body_dm_mo_diff + END_DOC + integer :: i,j,istate + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:) + + + one_body_dm_mo_diff_eigvectors(1:mo_tot_num, 1:mo_tot_num, 2:N_states) =& + one_body_dm_mo_diff(1:mo_tot_num, 1:mo_tot_num, 2:N_states) + + n = mo_tot_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork), iwork(liwork)) + + lwork=-1 + liwork=-1 + istate=2 + + call dsyevd( 'V', 'U', mo_tot_num, & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), & + one_body_dm_mo_diff_eigvalues(1,istate), & + work, lwork, iwork, liwork, info) + + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(iwork,work) + + allocate(work(lwork), iwork(liwork)) + + do istate=2,N_states + call dsyevd( 'V', 'U', mo_tot_num, & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), & + one_body_dm_mo_diff_eigvalues(1,istate), & + work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + + enddo + + deallocate(iwork,work) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, transition_natorb, (ao_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Natural transition molecular orbitals + END_DOC + + integer :: istate + + do istate=2,N_states + call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, & + mo_coef, size(mo_coef,1), & + one_body_dm_mo_diff_eigvectors(1,1,istate), & + size(one_body_dm_mo_diff_eigvectors,1), 0.d0, & + transition_natorb(1,1,istate), size(transition_natorb,1)) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, phi_s, (2:N_states) ] + implicit none + BEGIN_DOC + ! + END_DOC + + integer :: i,istate + double precision, allocatable :: T(:,:), A(:,:), D(:,:) + double precision :: trace, norm + allocate(T(ao_num,ao_num), A(ao_num,ao_num), D(ao_num,ao_num)) + + do istate=2,N_states + + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + S_half, size(S_half,1), & + one_body_dm_ao_attachment(1,1,istate), size(one_body_dm_ao_attachment,1), 0.d0,& + T, size(T,1)) + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + T, size(T,1), & + S_half, size(S_half,1), 0.d0, & + A, size(A,1)) +! + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + S_half, size(S_half,1), & + one_body_dm_ao_detachment(1,1,istate), size(one_body_dm_ao_detachment,1), 0.d0,& + T, size(T,1)) + call dgemm('N','N',ao_num,ao_num,ao_num,1.d0, & + T, size(T,1), & + S_half, size(S_half,1), 0.d0, & + D, size(D,1)) + + trace = 0.d0 + do i=1,ao_num + trace = trace + A(i,i) + enddo + norm = 0.d0 + do i=1,ao_num + norm = norm + D(i,i) + enddo + norm = 0.5d0*(norm + trace) + + trace = 0.d0 + do i=1,mo_tot_num + trace = trace + dsqrt(A(i,i)*D(i,i)) + enddo + phi_s(istate) = trace/norm + enddo + +END_PROVIDER + diff --git a/src/AO_Basis/ao_overlap.irp.f b/src/AO_Basis/ao_overlap.irp.f index edf48b25..15d8a810 100644 --- a/src/AO_Basis/ao_overlap.irp.f +++ b/src/AO_Basis/ao_overlap.irp.f @@ -129,3 +129,41 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ] !$OMP END PARALLEL DO END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, S_half, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! S^{1/2} + END_DOC + + integer :: i,j,k + double precision, allocatable :: U(:,:) + double precision, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + + allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num)) + + call svd(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num) + + do i=1,ao_num + D(i) = dsqrt(D(i)) + do j=1,ao_num + S_half(j,i) = 0.d0 + enddo + enddo + + do k=1,ao_num + do j=1,ao_num + do i=1,ao_num + S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + enddo + + deallocate(U,Vt,D) + +END_PROVIDER + + diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index e0764d96..4fa33daa 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -15,6 +15,32 @@ enddo END_PROVIDER +BEGIN_PROVIDER [ double precision, one_body_dm_mo_diff, (mo_tot_num,mo_tot_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Difference of the one-body density matrix with respect to the ground state + END_DOC + integer :: i,j, istate + + do istate=2,N_states + do j=1,mo_tot_num + do i=1,mo_tot_num + one_body_dm_mo_diff(i,j,istate) = & + one_body_dm_mo_alpha(i,j,istate) - one_body_dm_mo_alpha(i,j,1) + & + one_body_dm_mo_beta (i,j,istate) - one_body_dm_mo_beta (i,j,1) + enddo + enddo + double precision :: trace + trace = 0.d0 + do i=1,mo_tot_num + trace += one_body_dm_mo_diff(i,i,istate) + enddo + print *, irp_here, trace + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_dm_mo_spin_index, (mo_tot_num_align,mo_tot_num,N_states,2) ] implicit none integer :: i,j,ispin,istate diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 1a3d1d36..7e81756c 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -69,7 +69,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) - double precision, allocatable :: S_half(:,:) + double precision, allocatable :: S(:,:) !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D integer :: info, i, j @@ -77,7 +77,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) return endif - allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n)) + allocate (U(ldc,n), Vt(lda,n), D(n), S(lda,n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) @@ -103,13 +103,13 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S_half,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j) !$OMP DO do j=1,n do i=1,n - S_half(i,j) = U(i,j)*D(j) + S(i,j) = U(i,j)*D(j) enddo do i=1,n U(i,j) = C(i,j) @@ -119,8 +119,8 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) !$OMP END PARALLEL - call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) - deallocate (U, Vt, D, S_half) + call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) + deallocate (U, Vt, D, S) end @@ -210,19 +210,19 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) double precision, allocatable :: U(:,:) double precision, allocatable :: Vt(:,:) double precision, allocatable :: D(:) - double precision, allocatable :: S_half(:,:) + double precision, allocatable :: S(:,:) integer :: info, i, j, k if (n < 2) then return endif - allocate(U(ldc,n),Vt(lda,n),S_half(lda,n),D(n)) + allocate(U(ldc,n),Vt(lda,n),S(lda,n),D(n)) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(S_half,U,D,Vt,n,C,m) & + !$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j,k) !$OMP DO @@ -233,7 +233,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) D(i) = 1.d0/dsqrt(D(i)) endif do j=1,n - S_half(j,i) = 0.d0 + S(j,i) = 0.d0 enddo enddo !$OMP END DO @@ -243,7 +243,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do j=1,n do i=1,n - S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j) + S(i,j) = S(i,j) + U(i,k)*D(k)*Vt(k,j) enddo enddo !$OMP END DO NOWAIT @@ -261,9 +261,9 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP END PARALLEL - call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) + call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S,size(S,1),0.d0,C,size(C,1)) - deallocate(U,Vt,S_half,D) + deallocate(U,Vt,S,D) end