From 541d7f5ff91b246bd9454d4f14879ca54470e837 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 3 Oct 2023 20:04:34 +0200 Subject: [PATCH] added attachment orbitals --- src/determinants/density_matrix.irp.f | 98 +++++++++++++++ src/determinants/dipole_moments.irp.f | 16 +-- src/tc_bi_ortho/tc_bi_ortho_prop.irp.f | 15 +++ src/tc_bi_ortho/tc_prop.irp.f | 1 + src/tools/attachement_orb.irp.f | 168 +++++++++++++++++++++++++ 5 files changed, 290 insertions(+), 8 deletions(-) create mode 100644 src/tools/attachement_orb.irp.f diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index ce4d96c2..46726df0 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -493,3 +493,101 @@ subroutine get_occupation_from_dets(istate,occupation) 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 diff --git a/src/determinants/dipole_moments.irp.f b/src/determinants/dipole_moments.irp.f index e445c56b..dae04369 100644 --- a/src/determinants/dipole_moments.irp.f +++ b/src/determinants/dipole_moments.irp.f @@ -26,10 +26,10 @@ enddo enddo -! print*,'electron part for z_dipole = ',z_dipole_moment -! print*,'electron part for y_dipole = ',y_dipole_moment -! print*,'electron part for x_dipole = ',x_dipole_moment -! + print*,'electron part for z_dipole = ',z_dipole_moment + print*,'electron part for y_dipole = ',y_dipole_moment + print*,'electron part for x_dipole = ',x_dipole_moment + nuclei_part_z = 0.d0 nuclei_part_y = 0.d0 nuclei_part_x = 0.d0 @@ -38,10 +38,10 @@ nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) enddo -! print*,'nuclei part for z_dipole = ',nuclei_part_z -! print*,'nuclei part for y_dipole = ',nuclei_part_y -! print*,'nuclei part for x_dipole = ',nuclei_part_x -! + print*,'nuclei part for z_dipole = ',nuclei_part_z + print*,'nuclei part for y_dipole = ',nuclei_part_y + print*,'nuclei part for x_dipole = ',nuclei_part_x + do istate = 1, N_states z_dipole_moment(istate) += nuclei_part_z y_dipole_moment(istate) += nuclei_part_y diff --git a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f index 9168fb3d..a5fe9249 100644 --- a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f @@ -34,4 +34,19 @@ subroutine test do i= 1, 3 print*,tc_bi_ortho_dipole(i,1) enddo + integer, allocatable :: occ(:,:) + integer :: n_occ_ab(2) + allocate(occ(N_int*bit_kind_size,2)) + call bitstring_to_list_ab(ref_bitmask, occ, n_occ_ab, N_int) + integer :: ispin,j,jorb + double precision :: accu + accu = 0.d0 + do ispin=1, 2 + do i = 1, n_occ_ab(ispin) + jorb = occ(i,ispin) + accu += mo_bi_orth_bipole_z(jorb,jorb) + enddo + enddo + print*,'accu = ',accu + end diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f index a13dc9a2..3375fed6 100644 --- a/src/tc_bi_ortho/tc_prop.irp.f +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -90,6 +90,7 @@ enddo enddo enddo + print*,'tc_bi_ortho_dipole(3) elec = ',tc_bi_ortho_dipole(3,1) nuclei_part = 0.d0 do m = 1, 3 diff --git a/src/tools/attachement_orb.irp.f b/src/tools/attachement_orb.irp.f new file mode 100644 index 00000000..92a51ca8 --- /dev/null +++ b/src/tools/attachement_orb.irp.f @@ -0,0 +1,168 @@ +program molden_detachment_attachment + implicit none + read_wf=.True. + touch read_wf + call molden_attachment +end + +subroutine molden_attachment + implicit none + BEGIN_DOC + ! Produces a Molden file + END_DOC + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + integer :: i,j,k,l + double precision, parameter :: a0 = 0.529177249d0 + + PROVIDE ezfio_filename + + output=trim(ezfio_filename)//'.attachement.mol' + print*,'output = ',trim(output) + + i_unit_output = getUnitAndOpen(output,'w') + + write(i_unit_output,'(A)') '[Molden Format]' + + write(i_unit_output,'(A)') '[Atoms] Angs' + do i = 1, nucl_num + write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') & + trim(element_name(int(nucl_charge(i)))), & + i, & + int(nucl_charge(i)), & + nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0 + enddo + + write(i_unit_output,'(A)') '[GTO]' + + character*(1) :: character_shell + integer :: i_shell,i_prim,i_ao + integer :: iorder(ao_num) + integer :: nsort(ao_num) + + i_shell = 0 + i_prim = 0 + do i=1,nucl_num + write(i_unit_output,*) i, 0 + do j=1,nucl_num_shell_aos(i) + i_shell +=1 + i_ao = nucl_list_shell_aos(i,j) + character_shell = trim(ao_l_char(i_ao)) + write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00' + do k = 1, ao_prim_num(i_ao) + i_prim +=1 + write(i_unit_output,'(ES20.10,2X,ES20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k) + enddo + l = i_ao + do while ( ao_l(l) == ao_l(i_ao) ) + nsort(l) = i*10000 + j*100 + l += 1 + if (l > ao_num) exit + enddo + enddo + write(i_unit_output,*)'' + enddo + + + do i=1,ao_num + iorder(i) = i + ! p + if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 1 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 2 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 3 + ! d + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 1 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 2 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 3 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 4 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 5 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 6 + ! f + else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 1 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 2 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then + nsort(i) += 3 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 4 + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 5 + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 6 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 7 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 8 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 9 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 10 + ! g + else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 1 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 2 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then + nsort(i) += 3 + else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 4 + else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 5 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 6 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 7 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then + nsort(i) += 8 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then + nsort(i) += 9 + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then + nsort(i) += 10 + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 11 + else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 12 + else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 13 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then + nsort(i) += 14 + else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then + nsort(i) += 15 + endif + enddo + + call isort(nsort,iorder,ao_num) + write(i_unit_output,'(A)') '[MO]' + integer :: istate + istate = 2 + do i=1,n_dettachment(istate) + write (i_unit_output,*) 'Sym= 1' + write (i_unit_output,*) 'Ene=', dettachment_numbers_sorted(i,istate) + write (i_unit_output,*) 'Spin= Alpha' + write (i_unit_output,*) 'Occup=', dettachment_numbers_sorted(i,istate) + do j=1,ao_num + write(i_unit_output, '(I6,2X,ES20.10)') j, dettachment_orbitals(iorder(j),i,istate) + enddo + enddo + do i=1,n_attachment(istate) + write (i_unit_output,*) 'Sym= 1' + write (i_unit_output,*) 'Ene=', attachment_numbers_sorted(i,istate) + write (i_unit_output,*) 'Spin= Alpha' + write (i_unit_output,*) 'Occup=', attachment_numbers_sorted(i,istate) + do j=1,ao_num + write(i_unit_output, '(I6,2X,ES20.10)') j, attachment_orbitals(iorder(j),i,istate) + enddo + enddo + close(i_unit_output) +end +