PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 12:23:43 +01:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
Anthony Scemama 2023-10-05 15:05:07 +02:00
commit b5f1116fca
5 changed files with 290 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

98
src/determinants/density_matrix.irp.f
View File

@ -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

16
src/determinants/dipole_moments.irp.f
View File

@ -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

15
src/tc_bi_ortho/tc_bi_ortho_prop.irp.f
View File

@ -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

1
src/tc_bi_ortho/tc_prop.irp.f
View File

@ -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

168
src/tools/attachement_orb.irp.f Normal file
View File

@ -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