LCPQ
/
quantum_package
mirror of https://github.com/LCPQ/quantum_package
10
0
Fork 0
Code Issues Releases Wiki Activity

added the mulliken and hyperfine coupling constants analysis

This commit is contained in:
Emmanuel Giner 2016-02-16 18:32:53 +01:00
parent 6bbc24c1ec
commit 19e276fc0d
11 changed files with 407 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/Hartree_Fock/.gitignore vendored
View File

@ -5,7 +5,6 @@ AO_Basis
Bitmask
Electrons
Ezfio_files
Huckel_guess
IRPF90_man
IRPF90_temp
Integrals_Bielec
@ -16,7 +15,6 @@ Makefile
Makefile.depend
Nuclei
Pseudo
SCF
Utils
ZMQ
ezfio_interface.irp.f

135
plugins/Properties/hyperfine_constants.irp.f Normal file
View File

@ -0,0 +1,135 @@
BEGIN_PROVIDER [double precision, spin_density_at_nucleous, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
do k = 1, ao_num
accu += one_body_spin_density_ao(k,j) * aos_array(k) * aos_array(j)
enddo
enddo
spin_density_at_nucleous(i) = accu
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_from_mo, (nucl_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_per_mo, (nucl_num,mo_tot_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: contrib
double precision :: mo_values(mo_tot_num)
accu = 0.d0
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
spin_density_at_nucleous_from_mo(i) = 0.d0
do k = 1, mo_tot_num
mo_values(k) = 0.d0
do j = 1, ao_num
mo_values(k) += mo_coef(j,k) * aos_array(j)
enddo
enddo
do k = 1, mo_tot_num
spin_density_at_nucleous_contrib_per_mo(i,k) = 0.d0
do m = 1, mo_tot_num
contrib = one_body_spin_density_mo(k,m) * mo_values(k) * mo_values(m)
spin_density_at_nucleous_from_mo(i) += contrib
spin_density_at_nucleous_contrib_per_mo(i,k) += contrib
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo, (nucl_num,mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, spin_density_at_nucleous_contrib_mo_test, (nucl_num)]
implicit none
BEGIN_DOC
! value of the spin density at each nucleus
END_DOC
integer :: i,j,k,l,m
spin_density_at_nucleous_contrib_mo_test = 0.d0
do i = 1, nucl_num
double precision :: r(3),accu,aos_array(ao_num)
double precision :: c_i1,c_j1
r(1:3) = nucl_coord(i,1:3)
call give_all_aos_at_r(r,aos_array)
do k = 1, mo_tot_num
do m = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
c_i1 = mo_coef(j,k)
do l = 1, ao_num
c_j1 = c_i1*mo_coef(l,m)
accu += one_body_spin_density_mo(k,m) * aos_array(l) * aos_array(j) * c_j1
enddo
enddo
spin_density_at_nucleous_contrib_mo(i,k,m) = accu
spin_density_at_nucleous_contrib_mo_test(i) += accu
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, conversion_factor_mhz_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_gauss_hcc, (100)]
&BEGIN_PROVIDER [double precision, conversion_factor_cm_1_hcc, (100)]
BEGIN_DOC
! Conversion factor for the calculation of the hcc, according to the nuclear charge
END_DOC
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
conversion_factor_mhz_hcc =0.d0
! hydrogen
conversion_factor_mhz_hcc(1) = 4469.84692227102460d0
conversion_factor_gauss_hcc(1) = 1594.95296390862904d0
conversion_factor_cm_1_hcc(1) = 1490.98044430157870d0
! Li
conversion_factor_mhz_hcc(3) = 1737.2746512855997d0
conversion_factor_gauss_hcc(3) = 619.9027742370165d0
conversion_factor_cm_1_hcc(3) = 579.4924475562677d0
! carbon
conversion_factor_mhz_hcc(6) = 1124.18303629792945d0
conversion_factor_gauss_hcc(6) = 401.136570647523058d0
conversion_factor_cm_1_hcc(6) = 374.987097339830086d0
! nitrogen
conversion_factor_mhz_hcc(7) = 323.102093833793390d0
conversion_factor_gauss_hcc(7) = 115.290892768082614d0
conversion_factor_cm_1_hcc(7) = 107.775257586297698d0
! Oxygen
conversion_factor_mhz_hcc(8) = -606.1958551736545d0
conversion_factor_gauss_hcc(8) = -216.30574771560407d0
conversion_factor_cm_1_hcc(8) = -202.20517197179822d0
END_PROVIDER
BEGIN_PROVIDER [double precision, iso_hcc_mhz, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_gauss, (nucl_num)]
&BEGIN_PROVIDER [double precision, iso_hcc_cm_1, (nucl_num)]
BEGIN_DOC
! isotropic hyperfine coupling constants among the various atoms
END_DOC
integer :: i
do i = 1, nucl_num
iso_hcc_mhz(i) = conversion_factor_mhz_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !* 0.5d0
iso_hcc_gauss(i) = conversion_factor_gauss_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i)!* 0.5d0
iso_hcc_cm_1(i) = conversion_factor_cm_1_hcc(nint(nucl_charge(i))) * spin_density_at_nucleous(i) !*0.5d0
enddo
END_PROVIDER

107
plugins/Properties/mulliken.irp.f Normal file
View File

@ -0,0 +1,107 @@
BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
spin_population = 0.d0
do i = 1, ao_num
do j = 1, ao_num
spin_population(j,i) = one_body_spin_density_ao(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)]
implicit none
integer :: i
double precision :: accu
spin_population_angular_momentum = 0.d0
do i = 1, ao_num
spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, spin_gross_orbital_product, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product for the spin population
END_DOC
do i = 1, ao_num
do j = 1, ao_num
spin_gross_orbital_product(i) += spin_population(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_spin_densities, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!ATOMIC SPIN POPULATION (ALPHA MINUS BETA)
END_DOC
mulliken_spin_densities = 0.d0
do i = 1, ao_num
mulliken_spin_densities(ao_nucl(i)) += spin_gross_orbital_product(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, electronic_population_alpha, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [double precision, electronic_population_beta, (ao_num_align,ao_num)]
implicit none
integer :: i,j
BEGIN_DOC
! spin population on the ao basis :
! spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
END_DOC
electronic_population_alpha = 0.d0
electronic_population_beta = 0.d0
do i = 1, ao_num
do j = 1, ao_num
electronic_population_alpha(j,i) = one_body_dm_ao_alpha(i,j) * ao_overlap(i,j)
electronic_population_beta(j,i) = one_body_dm_ao_beta(i,j) * ao_overlap(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, gross_orbital_product_alpha, (ao_num)]
&BEGIN_PROVIDER [double precision, gross_orbital_product_beta, (ao_num)]
implicit none
spin_gross_orbital_product = 0.d0
integer :: i,j
BEGIN_DOC
! gross orbital product
END_DOC
do i = 1, ao_num
do j = 1, ao_num
gross_orbital_product_alpha(i) += electronic_population_alpha(j,i)
gross_orbital_product_beta(i) += electronic_population_beta(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, mulliken_densities_alpha, (nucl_num)]
&BEGIN_PROVIDER [double precision, mulliken_densities_beta, (nucl_num)]
implicit none
integer :: i,j
BEGIN_DOC
!
END_DOC
mulliken_densities_alpha = 0.d0
mulliken_densities_beta = 0.d0
do i = 1, ao_num
mulliken_densities_alpha(ao_nucl(i)) += gross_orbital_product_alpha(i)
mulliken_densities_beta(ao_nucl(i)) += gross_orbital_product_beta(i)
enddo
END_PROVIDER

35
plugins/Properties/print_mulliken.irp.f Normal file
View File

@ -0,0 +1,35 @@
program print_mulliken
implicit none
read_wf = .True.
touch read_wf
print*,'Mulliken spin densities'
call test
end
subroutine test
double precision :: accu
integer :: i
integer :: j
accu= 0.d0
do i = 1, nucl_num
print*,i,nucl_charge(i),mulliken_spin_densities(i)
accu += mulliken_spin_densities(i)
enddo
print*,'Sum of Mulliken SD = ',accu
print*,'AO SPIN POPULATIONS'
accu = 0.d0
do i = 1, ao_num
accu += spin_gross_orbital_product(i)
write(*,'(X,I3,X,A4,X,I2,X,A4,X,F10.7)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),trim(l_to_charater(ao_l(i))),spin_gross_orbital_product(i)
enddo
print*,'sum = ',accu
accu = 0.d0
print*,'Angular momentum analysis'
do i = 0, ao_l_max
accu += spin_population_angular_momentum(i)
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu
enddo
end

4
scripts/generate_h_apply.py
View File

@ -156,11 +156,11 @@ class H_apply(object):
def filter_only_1h1p(self):
self["filter_only_1h1p_single"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(hole).eq..False.) cycle
if (is_a_1h1p(hole).eqv..False.) cycle
"""
self["filter_only_1h1p_double"] = """
! ! DIR$ FORCEINLINE
if (is_a_1h1p(key).eq..False.) cycle
if (is_a_1h1p(key).eqv..False.) cycle
"""

27
src/AO_Basis/aos.irp.f
View File

@ -146,3 +146,30 @@ integer function ao_power_index(nx,ny,nz)
ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
end
BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
&BEGIN_PROVIDER [ integer, ao_l_max ]
&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ]
implicit none
BEGIN_DOC
! ao_l = l value of the AO: a+b+c in x^a y^b z^c
END_DOC
integer :: i
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
ao_l_char(i) = l_to_charater(ao_l(i))
enddo
ao_l_max = maxval(ao_l)
END_PROVIDER
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
BEGIN_DOC
! character corresponding to the "L" value of an AO orbital
END_DOC
implicit none
l_to_charater(0)='S'
l_to_charater(1)='P'
l_to_charater(2)='D'
l_to_charater(3)='F'
l_to_charater(4)='G'
END_PROVIDER

48
src/AO_Basis/aos_value.irp.f Normal file
View File

@ -0,0 +1,48 @@
double precision function ao_value(i,r)
implicit none
BEGIN_DOC
! return the value of the ith ao at point r
END_DOC
double precision, intent(in) :: r(3)
integer, intent(in) :: i
integer :: m,num_ao
double precision :: center_ao(3)
double precision :: beta
integer :: power_ao(3)
num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3)
double precision :: accu,dx,dy,dz,r2
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2))
dz = (r(3) - center_ao(3))
r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1)
dy = dy**power_ao(2)
dz = dz**power_ao(3)
accu = 0.d0
do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
enddo
ao_value = accu * dx * dy * dz
end
subroutine give_all_aos_at_r(r,aos_array)
implicit none
BEGIN_dOC
! gives the values of aos at a given point r
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num)
integer :: i
double precision :: ao_value
do i = 1, ao_num
aos_array(i) = ao_value(i,r)
enddo
end

51
src/Determinants/density_matrix.irp.f
View File

@ -206,3 +206,54 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
state_average_weight = 1.d0/dble(N_states)
END_PROVIDER
BEGIN_PROVIDER [ double precision, one_body_spin_density_ao, (ao_num_align,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_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_spin_density_mo(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_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_body_dm_ao_alpha, (ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, one_body_dm_ao_beta, (ao_num_align,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 :: dm_mo
one_body_spin_density_ao = 0.d0
do k = 1, ao_num
do l = 1, ao_num
do i = 1, mo_tot_num
do j = 1, mo_tot_num
dm_mo = one_body_dm_mo_alpha(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle
one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo
enddo
enddo
enddo
enddo
END_PROVIDER

3
src/Integrals_Bielec/.gitignore vendored
View File

@ -17,5 +17,4 @@ ZMQ
ezfio_interface.irp.f
irpf90.make
irpf90_entities
tags
test_integrals
tags

2
src/Integrals_Monoelec/.gitignore vendored
View File

@ -12,9 +12,7 @@ Makefile.depend
Nuclei
Pseudo
Utils
check_orthonormality
ezfio_interface.irp.f
irpf90.make
irpf90_entities
save_ortho_mos
tags

5
src/MOGuess/.gitignore vendored
View File

@ -4,7 +4,6 @@
AO_Basis
Electrons
Ezfio_files
H_CORE_guess
IRPF90_man
IRPF90_temp
Integrals_Monoelec
@ -15,8 +14,6 @@ Nuclei
Pseudo
Utils
ezfio_interface.irp.f
guess_overlap
irpf90.make
irpf90_entities
tags
truncate_mos
tags