BEGIN_PROVIDER [ double precision, one_e_dm_mo_alpha_average, (mo_num,mo_num) ] &BEGIN_PROVIDER [ double precision, one_e_dm_mo_beta_average, (mo_num,mo_num) ] implicit none BEGIN_DOC ! $\alpha$ and $\beta$ one-body density matrix for each state END_DOC integer :: i one_e_dm_mo_alpha_average = 0.d0 one_e_dm_mo_beta_average = 0.d0 do i = 1,N_states one_e_dm_mo_alpha_average(:,:) += one_e_dm_mo_alpha(:,:,i) * state_average_weight(i) one_e_dm_mo_beta_average(:,:) += one_e_dm_mo_beta(:,:,i) * state_average_weight(i) enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_mo_diff, (mo_num,mo_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_num do i=1,mo_num one_e_dm_mo_diff(i,j,istate) = & one_e_dm_mo_alpha(i,j,istate) - one_e_dm_mo_alpha(i,j,1) +& one_e_dm_mo_beta (i,j,istate) - one_e_dm_mo_beta (i,j,1) enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_mo_spin_index, (mo_num,mo_num,N_states,2) ] implicit none integer :: i,j,ispin,istate ispin = 1 do istate = 1, N_states do j = 1, mo_num do i = 1, mo_num one_e_dm_mo_spin_index(i,j,istate,ispin) = one_e_dm_mo_alpha(i,j,istate) enddo enddo enddo ispin = 2 do istate = 1, N_states do j = 1, mo_num do i = 1, mo_num one_e_dm_mo_spin_index(i,j,istate,ispin) = one_e_dm_mo_beta(i,j,istate) enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_dagger_mo_spin_index, (mo_num,mo_num,N_states,2) ] implicit none integer :: i,j,ispin,istate ispin = 1 do istate = 1, N_states do j = 1, mo_num one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_alpha(j,j,istate) do i = j+1, mo_num one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) enddo enddo enddo ispin = 2 do istate = 1, N_states do j = 1, mo_num one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_beta(j,j,istate) do i = j+1, mo_num one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_mo_alpha, (mo_num,mo_num,N_states) ] &BEGIN_PROVIDER [ double precision, one_e_dm_mo_beta, (mo_num,mo_num,N_states) ] implicit none BEGIN_DOC ! $\alpha$ and $\beta$ one-body density matrix for each state END_DOC integer :: j,k,l,m,k_a,k_b 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(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int) integer :: exc(0:2,2),n_occ(2) double precision, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:) integer :: krow, kcol, lrow, lcol PROVIDE psi_det one_e_dm_mo_alpha = 0.d0 one_e_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,& !$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)& !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num, & !$OMP elec_beta_num,one_e_dm_mo_alpha,one_e_dm_mo_beta,N_det,& !$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,& !$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,& !$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,& !$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,& !$OMP N_det_alpha_unique,N_det_beta_unique,irp_here) allocate(tmp_a(mo_num,mo_num,N_states), tmp_b(mo_num,mo_num,N_states) ) tmp_a = 0.d0 !$OMP DO SCHEDULE(guided) do k_a=1,N_det krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) ! Diagonal part ! ------------- call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) do m=1,N_states ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,m) do l=1,elec_alpha_num j = occ(l,1) tmp_a(j,j,m) += ck enddo enddo if (k_a == N_det) cycle l = k_a+1 lrow = psi_bilinear_matrix_rows(l) lcol = psi_bilinear_matrix_columns(l) ! Fix beta determinant, loop over alphas do while ( lcol == kcol ) tmp_det2(:) = psi_det_alpha_unique(:, lrow) call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int) if (degree == 1) then exc = 0 call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int) call decode_exc_spin(exc,h1,p1,h2,p2) do m=1,N_states ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,m) * phase tmp_a(h1,p1,m) += ckl tmp_a(p1,h1,m) += ckl enddo endif l = l+1 if (l>N_det) exit lrow = psi_bilinear_matrix_rows(l) lcol = psi_bilinear_matrix_columns(l) enddo enddo !$OMP END DO NOWAIT !$OMP CRITICAL one_e_dm_mo_alpha(:,:,:) = one_e_dm_mo_alpha(:,:,:) + tmp_a(:,:,:) !$OMP END CRITICAL deallocate(tmp_a) tmp_b = 0.d0 !$OMP DO SCHEDULE(guided) do k_b=1,N_det krow = psi_bilinear_matrix_transp_rows(k_b) ASSERT (krow <= N_det_alpha_unique) kcol = psi_bilinear_matrix_transp_columns(k_b) ASSERT (kcol <= N_det_beta_unique) tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) ! Diagonal part ! ------------- call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) do m=1,N_states ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,m) do l=1,elec_beta_num j = occ(l,2) tmp_b(j,j,m) += ck enddo enddo if (k_b == N_det) cycle l = k_b+1 lrow = psi_bilinear_matrix_transp_rows(l) lcol = psi_bilinear_matrix_transp_columns(l) ! Fix beta determinant, loop over alphas do while ( lrow == krow ) tmp_det2(:) = psi_det_beta_unique(:, lcol) call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int) if (degree == 1) then exc = 0 call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int) call decode_exc_spin(exc,h1,p1,h2,p2) do m=1,N_states ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,m) * phase tmp_b(h1,p1,m) += ckl tmp_b(p1,h1,m) += ckl enddo endif l = l+1 if (l>N_det) exit lrow = psi_bilinear_matrix_transp_rows(l) lcol = psi_bilinear_matrix_transp_columns(l) enddo enddo !$OMP END DO NOWAIT !$OMP CRITICAL one_e_dm_mo_beta(:,:,:) = one_e_dm_mo_beta(:,:,:) + tmp_b(:,:,:) !$OMP END CRITICAL deallocate(tmp_b) !$OMP END PARALLEL END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_mo, (mo_num,mo_num) ] implicit none BEGIN_DOC ! One-body density matrix END_DOC one_e_dm_mo = one_e_dm_mo_alpha_average + one_e_dm_mo_beta_average END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_spin_density_mo, (mo_num,mo_num) ] implicit none BEGIN_DOC ! $\rho(\alpha) - \rho(\beta)$ END_DOC one_e_spin_density_mo = one_e_dm_mo_alpha_average - one_e_dm_mo_beta_average END_PROVIDER subroutine set_natural_mos implicit none BEGIN_DOC ! Set natural orbitals, obtained by diagonalization of the one-body density matrix ! in the |MO| basis END_DOC character*(64) :: label double precision, allocatable :: tmp(:,:) label = "Natural" integer :: i,j,iorb,jorb do i = 1, n_virt_orb iorb = list_virt(i) do j = 1, n_core_inact_act_orb jorb = list_core_inact_act(j) enddo enddo call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) soft_touch mo_occ end subroutine save_natural_mos_canon_label implicit none BEGIN_DOC ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in ! the |MO| basis END_DOC call set_natural_mos_canon_label call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) call orthonormalize_mos call save_mos end subroutine set_natural_mos_canon_label implicit none BEGIN_DOC ! Set natural orbitals, obtained by diagonalization of the one-body density matrix ! in the |MO| basis END_DOC character*(64) :: label double precision, allocatable :: tmp(:,:) label = "Canonical" integer :: i,j,iorb,jorb do i = 1, n_virt_orb iorb = list_virt(i) do j = 1, n_core_inact_act_orb jorb = list_core_inact_act(j) enddo enddo call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) soft_touch mo_occ end subroutine set_natorb_no_ov_rot implicit none BEGIN_DOC ! Set natural orbitals, obtained by diagonalization of the one-body density matrix ! in the |MO| basis END_DOC character*(64) :: label double precision, allocatable :: tmp(:,:) allocate(tmp(mo_num, mo_num)) label = "Natural" tmp = one_e_dm_mo integer :: i,j,iorb,jorb do i = 1, n_virt_orb iorb = list_virt(i) do j = 1, n_core_inact_act_orb jorb = list_core_inact_act(j) tmp(iorb, jorb) = 0.d0 tmp(jorb, iorb) = 0.d0 enddo enddo call mo_as_svd_vectors_of_mo_matrix_eig(tmp,size(tmp,1),mo_num,mo_num,mo_occ,label) soft_touch mo_occ end subroutine save_natural_mos_no_ov_rot implicit none BEGIN_DOC ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in ! the |MO| basis END_DOC call set_natorb_no_ov_rot call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) call orthonormalize_mos call save_mos end subroutine save_natural_mos implicit none BEGIN_DOC ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in ! the |MO| basis END_DOC call set_natural_mos call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) call orthonormalize_mos call save_mos end BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ] implicit none BEGIN_DOC ! Weight of the states in the selection : $\frac{1}{c_0^2}$. END_DOC if (N_states > 1) then integer :: i double precision :: c do i=1,N_states c0_weight(i) = 1.d-31 c = maxval(psi_coef(:,i) * psi_coef(:,i)) c0_weight(i) = 1.d0/(c+1.d-20) enddo c = 1.d0/sum(c0_weight(:)) do i=1,N_states c0_weight(i) = c0_weight(i) * c enddo else c0_weight(:) = 1.d0 endif END_PROVIDER BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] implicit none BEGIN_DOC ! Weights in the state-average calculation of the density matrix END_DOC logical :: exists state_average_weight(:) = 1.d0 if (weight_one_e_dm == 0) then state_average_weight(:) = c0_weight(:) else if (weight_one_e_dm == 1) then state_average_weight(:) = 1.d0/N_states else call ezfio_has_determinants_state_average_weight(exists) if (exists) then call ezfio_get_determinants_state_average_weight(state_average_weight) endif endif state_average_weight(:) = state_average_weight(:)+1.d-31 state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:))) END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_spin_density_ao, (ao_num,ao_num) ] BEGIN_DOC ! One body spin density matrix on the |AO| basis : $\rho_{AO}(\alpha) - \rho_{AO}(\beta)$ END_DOC implicit none integer :: i,j,k,l double precision :: dm_mo one_e_spin_density_ao = 0.d0 do k = 1, ao_num do l = 1, ao_num do i = 1, mo_num do j = 1, mo_num dm_mo = one_e_spin_density_mo(j,i) ! if(dabs(dm_mo).le.1.d-10)cycle one_e_spin_density_ao(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo enddo enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_ao_alpha, (ao_num,ao_num) ] &BEGIN_PROVIDER [ double precision, one_e_dm_ao_beta, (ao_num,ao_num) ] BEGIN_DOC ! One body density matrix on the |AO| basis : $\rho_{AO}(\alpha), \rho_{AO}(\beta)$. END_DOC implicit none integer :: i,j,k,l double precision :: mo_alpha,mo_beta one_e_dm_ao_alpha = 0.d0 one_e_dm_ao_beta = 0.d0 do k = 1, ao_num do l = 1, ao_num do i = 1, mo_num do j = 1, mo_num mo_alpha = one_e_dm_mo_alpha_average(j,i) mo_beta = one_e_dm_mo_beta_average(j,i) ! if(dabs(dm_mo).le.1.d-10)cycle one_e_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha one_e_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_beta enddo enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_ao_alpha_nstates, (ao_num,ao_num,N_states) ] &BEGIN_PROVIDER [ double precision, one_e_dm_ao_beta_nstates, (ao_num,ao_num,N_states) ] BEGIN_DOC ! One body density matrix on the |AO| basis : $\rho_{AO}(\alpha), \rho_{AO}(\beta)$. END_DOC implicit none integer :: i,j,k,l,istate double precision :: mo_alpha,mo_beta one_e_dm_ao_alpha_nstates = 0.d0 one_e_dm_ao_beta_nstates = 0.d0 do istate = 1, N_states do k = 1, ao_num do l = 1, ao_num do i = 1, mo_num do j = 1, mo_num mo_alpha = one_e_dm_mo_alpha(j,i,istate) mo_beta = one_e_dm_mo_beta(j,i,istate) ! if(dabs(dm_mo).le.1.d-10)cycle one_e_dm_ao_alpha_nstates(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha one_e_dm_ao_beta_nstates(l,k,istate) += mo_coef(k,i) * mo_coef(l,j) * mo_beta enddo enddo enddo enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)] implicit none BEGIN_DOC ! one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta END_DOC one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta END_PROVIDER subroutine get_occupation_from_dets(istate,occupation) implicit none double precision, intent(out) :: occupation(mo_num) integer, intent(in) :: istate BEGIN_DOC ! Returns the average occupation of the MOs END_DOC integer :: i,j, ispin integer :: list(N_int*bit_kind_size,2) integer :: n_elements(2) double precision :: c, norm_2 ASSERT (istate > 0) ASSERT (istate <= N_states) occupation = 0.d0 double precision, external :: u_dot_u norm_2 = 1.d0/u_dot_u(psi_coef(1,istate),N_det) do i=1,N_det c = psi_coef(i,istate)*psi_coef(i,istate)*norm_2 call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int) do ispin=1,2 do j=1,n_elements(ispin) ASSERT ( list(j,ispin) < mo_num ) occupation( list(j,ispin) ) += c enddo enddo enddo end BEGIN_PROVIDER [double precision, difference_dm, (mo_num, mo_num, N_states)] implicit none BEGIN_DOC ! difference_dm(i,j,istate) = dm(i,j,1) - dm(i,j,istate) END_DOC integer :: istate do istate = 1, N_states difference_dm(:,:,istate) = one_e_dm_mo_alpha(:,:,1) + one_e_dm_mo_beta(:,:,1) & - (one_e_dm_mo_alpha(:,:,istate) + one_e_dm_mo_beta(:,:,istate)) enddo END_PROVIDER BEGIN_PROVIDER [double precision, difference_dm_eigvect, (mo_num, mo_num, N_states) ] &BEGIN_PROVIDER [double precision, difference_dm_eigval, (mo_num, N_states) ] implicit none BEGIN_DOC ! eigenvalues and eigevenctors of the difference_dm END_DOC integer :: istate,i do istate = 2, N_states call lapack_diag(difference_dm_eigval(1,istate),difference_dm_eigvect(1,1,istate)& ,difference_dm(1,1,istate),mo_num,mo_num) print*,'Eigenvalues of difference_dm for state ',istate do i = 1, mo_num print*,i,difference_dm_eigval(i,istate) enddo enddo END_PROVIDER BEGIN_PROVIDER [ integer , n_attachment, (N_states)] &BEGIN_PROVIDER [ integer , n_dettachment, (N_states)] &BEGIN_PROVIDER [ integer , list_attachment, (mo_num,N_states)] &BEGIN_PROVIDER [ integer , list_dettachment, (mo_num,N_states)] implicit none integer :: i,istate integer :: list_attachment_tmp(mo_num) n_attachment = 0 n_dettachment = 0 do istate = 2, N_states do i = 1, mo_num if(difference_dm_eigval(i,istate).lt.0.d0)then ! dettachment_orbitals n_dettachment(istate) += 1 list_dettachment(n_dettachment(istate),istate) = i ! they are already sorted else n_attachment(istate) += 1 list_attachment_tmp(n_attachment(istate)) = i ! they are not sorted endif enddo ! sorting the attachment do i = 0, n_attachment(istate) - 1 list_attachment(i+1,istate) = list_attachment_tmp(n_attachment(istate) - i) enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, attachment_numbers_sorted, (mo_num, N_states)] &BEGIN_PROVIDER [ double precision, dettachment_numbers_sorted, (mo_num, N_states)] implicit none integer :: i,istate do istate = 2, N_states print*,'dettachment' do i = 1, n_dettachment(istate) dettachment_numbers_sorted(i,istate) = difference_dm_eigval(list_dettachment(i,istate),istate) print*,i,list_dettachment(i,istate),dettachment_numbers_sorted(i,istate) enddo print*,'attachment' do i = 1, n_attachment(istate) attachment_numbers_sorted(i,istate) = difference_dm_eigval(list_attachment(i,istate),istate) print*,i,list_attachment(i,istate),attachment_numbers_sorted(i,istate) enddo enddo END_PROVIDER BEGIN_PROVIDER [ double precision, attachment_orbitals, (ao_num, mo_num, N_states)] &BEGIN_PROVIDER [ double precision, dettachment_orbitals, (ao_num, mo_num, N_states)] implicit none integer :: i,j,k,istate attachment_orbitals = 0.d0 dettachment_orbitals = 0.d0 do istate = 2, N_states do i = 1, n_dettachment(istate) do j = 1, mo_num do k = 1, ao_num dettachment_orbitals(k,list_dettachment(i,istate),istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_dettachment(i,istate),istate) enddo enddo enddo do i = 1, n_attachment(istate) do j = 1, mo_num do k = 1, ao_num attachment_orbitals(k,i,istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_attachment(i,istate),istate) enddo enddo enddo enddo END_PROVIDER