diff --git a/src/tc_bi_ortho/spin_mulliken.irp.f b/src/tc_bi_ortho/spin_mulliken.irp.f new file mode 100644 index 00000000..922cc1b9 --- /dev/null +++ b/src/tc_bi_ortho/spin_mulliken.irp.f @@ -0,0 +1,106 @@ + +BEGIN_PROVIDER [double precision, tc_spin_population, (ao_num,ao_num,N_states)] + implicit none + integer :: i,j,istate + BEGIN_DOC +! spin population on the ao basis : +! tc_spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * + END_DOC + tc_spin_population = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + do j = 1, ao_num + tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i) + enddo + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum, (0:ao_l_max,N_states)] +&BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum_per_atom, (0:ao_l_max,nucl_num,N_states)] + implicit none + integer :: i,istate + double precision :: accu + tc_spin_population_angular_momentum = 0.d0 + tc_spin_population_angular_momentum_per_atom = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + tc_spin_population_angular_momentum(ao_l(i),istate) += tc_spin_gross_orbital_product(i,istate) + tc_spin_population_angular_momentum_per_atom(ao_l(i),ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate) + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, tc_spin_gross_orbital_product, (ao_num,N_states)] + implicit none + tc_spin_gross_orbital_product = 0.d0 + integer :: i,j,istate + BEGIN_DOC +! gross orbital product for the spin population + END_DOC + do istate = 1, N_states + do i = 1, ao_num + do j = 1, ao_num + tc_spin_gross_orbital_product(i,istate) += tc_spin_population(j,i,istate) + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [double precision, tc_mulliken_spin_densities, (nucl_num,N_states)] + implicit none + integer :: i,j,istate + BEGIN_DOC +!ATOMIC SPIN POPULATION (ALPHA MINUS BETA) + END_DOC + tc_mulliken_spin_densities = 0.d0 + do istate = 1, N_states + do i = 1, ao_num + tc_mulliken_spin_densities(ao_nucl(i),istate) += tc_spin_gross_orbital_product(i,istate) + enddo + enddo + +END_PROVIDER + +subroutine tc_print_mulliken_sd + implicit none + double precision :: accu + integer :: i + integer :: j + print*,'Mulliken spin densities' + accu= 0.d0 + do i = 1, nucl_num + print*,i,nucl_charge(i),tc_mulliken_spin_densities(i,1) + accu += tc_mulliken_spin_densities(i,1) + enddo + print*,'Sum of Mulliken SD = ',accu + print*,'AO SPIN POPULATIONS' + accu = 0.d0 + do i = 1, ao_num + accu += tc_spin_gross_orbital_product(i,1) + write(*,'(1X,I3,1X,A4,1X,I2,1X,A4,1X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_character(ao_l(i))),tc_spin_gross_orbital_product(i,1) + enddo + print*,'sum = ',accu + accu = 0.d0 + print*,'Angular momentum analysis' + do i = 0, ao_l_max + accu += tc_spin_population_angular_momentum(i,1) + print*,' ',trim(l_to_character(i)),tc_spin_population_angular_momentum(i,1) + print*,'sum = ',accu + enddo + print*,'Angular momentum analysis per atom' + print*,'Angular momentum analysis' + do j = 1,nucl_num + accu = 0.d0 + do i = 0, ao_l_max + accu += tc_spin_population_angular_momentum_per_atom(i,j,1) + write(*,'(1X,I3,1X,A4,1X,A4,1X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_character(i)),tc_spin_population_angular_momentum_per_atom(i,j,1) + print*,'sum = ',accu + enddo + enddo + +end + diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index 33410570..b7e5ae81 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -21,7 +21,7 @@ allocate(dm_tmp(mo_num,mo_num), fock_diag(mo_num)) - dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) + dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1) print *, ' dm_tmp' do i = 1, mo_num diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f index c7f6c986..5bb0e2c0 100644 --- a/src/tc_bi_ortho/tc_prop.irp.f +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -1,8 +1,11 @@ -BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ] + BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_beta, (mo_num, mo_num,N_states,N_states) ] +&BEGIN_PROVIDER [ double precision, tc_transition_matrix_mo_alpha, (mo_num, mo_num,N_states,N_states) ] implicit none BEGIN_DOC - ! tc_transition_matrix(p,h,istate,jstate) = + ! tc_transition_matrix_mo_alpha(p,h,istate,jstate) = + ! + ! tc_transition_matrix_mo_beta(p,h,istate,jstate) = ! ! where are the left/right eigenvectors on a bi-ortho basis END_DOC @@ -11,43 +14,65 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state integer, allocatable :: occ(:,:) integer :: n_occ_ab(2),degree,exc(0:2,2,2) allocate(occ(N_int*bit_kind_size,2)) - tc_transition_matrix = 0.d0 - do istate = 1, N_states - do jstate = 1, N_states + tc_transition_matrix_mo_alpha = 0.d0 + tc_transition_matrix_mo_beta = 0.d0 do i = 1, N_det do j = 1, N_det call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) - if(degree.gt.1)then - cycle - else if (degree == 0)then - call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) - do p = 1, n_occ_ab(1) ! browsing the alpha electrons - m = occ(p,1) - tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - enddo - do p = 1, n_occ_ab(2) ! browsing the beta electrons - m = occ(p,1) - tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - enddo - else - call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int) - if (exc(0,1,1) == 1) then - ! Single alpha - h = exc(1,1,1) ! hole in psi_det(1,1,j) - p = exc(1,2,1) ! particle in psi_det(1,1,j) - else - ! Single beta - h = exc(1,1,2) ! hole in psi_det(1,1,j) - p = exc(1,2,2) ! particle in psi_det(1,1,j) - endif - tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) - endif + if(degree.gt.1)cycle + do istate = 1, N_states + do jstate = 1, N_states + if (degree == 0)then + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + do p = 1, n_occ_ab(1) ! browsing the alpha electrons + m = occ(p,1) + tc_transition_matrix_mo_alpha(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + do p = 1, n_occ_ab(2) ! browsing the beta electrons + m = occ(p,1) + tc_transition_matrix_mo_beta(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + else + call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Single alpha + h = exc(1,1,1) ! hole in psi_det(1,1,j) + p = exc(1,2,1) ! particle in psi_det(1,1,j) + tc_transition_matrix_mo_alpha(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + else + ! Single beta + h = exc(1,1,2) ! hole in psi_det(1,1,j) + p = exc(1,2,2) ! particle in psi_det(1,1,j) + tc_transition_matrix_mo_beta(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + endif + endif + enddo enddo - enddo enddo enddo END_PROVIDER + BEGIN_PROVIDER [double precision, tc_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ] + implicit none + BEGIN_DOC + ! tc_transition_matrix_mo(p,h,istate,jstate) = \sum_{sigma=alpha,beta} + ! + ! where are the left/right eigenvectors on a bi-ortho basis + END_DOC + tc_transition_matrix_mo = tc_transition_matrix_mo_beta + tc_transition_matrix_mo_alpha + END_PROVIDER + + + BEGIN_PROVIDER [double precision, tc_spin_transition_matrix_mo, (mo_num, mo_num,N_states,N_states) ] + implicit none + BEGIN_DOC + ! tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta + ! + ! where are the left/right eigenvectors on a bi-ortho basis + END_DOC + tc_spin_transition_matrix_mo = tc_transition_matrix_mo_alpha - tc_transition_matrix_mo_beta + END_PROVIDER + BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)] implicit none @@ -57,9 +82,9 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state do istate = 1, N_states do i = 1, mo_num do j = 1, mo_num - tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i) - tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i) - tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i) + tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i) + tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i) + tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix_mo(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i) enddo enddo enddo @@ -78,3 +103,55 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state enddo END_PROVIDER + + BEGIN_PROVIDER [ double precision, tc_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ] + implicit none + BEGIN_DOC +! tc_transition_matrix(p,h,istate,jstate) in the AO basis + END_DOC + integer :: i,j,k,l + double precision :: dm_mo + tc_transition_matrix_ao = 0.d0 + integer :: istate,jstate + do istate = 1, N_states + do jstate = 1, N_states + do i = 1, mo_num + do j = 1, mo_num + dm_mo = tc_transition_matrix_mo(j,i,jstate,istate) + do k = 1, ao_num + do l = 1, ao_num + tc_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo + enddo + enddo + enddo + enddo + enddo + enddo + + END_PROVIDER + + BEGIN_PROVIDER [ double precision, tc_spin_transition_matrix_ao, (ao_num, ao_num,N_states,N_states) ] + implicit none + BEGIN_DOC +! tc_spin_transition_matrix_ao(p,h,istate,jstate) in the AO basis + END_DOC + integer :: i,j,k,l + double precision :: dm_mo + tc_spin_transition_matrix_ao = 0.d0 + integer :: istate,jstate + do istate = 1, N_states + do jstate = 1, N_states + do i = 1, mo_num + do j = 1, mo_num + dm_mo = tc_spin_transition_matrix_mo(j,i,jstate,istate) + do k = 1, ao_num + do l = 1, ao_num + tc_spin_transition_matrix_ao(l,k,jstate,istate) += mo_l_coef(l,j) * mo_r_coef(k,i) * dm_mo + enddo + enddo + enddo + enddo + enddo + enddo + + END_PROVIDER