LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-22 10:47:33 +02:00
Code Issues Releases Wiki Activity

print_integrals_ao

Browse Source
This commit is contained in:
Anthony Scemama 2017-09-25 17:58:23 +02:00
parent 9776667a3a
commit 6b6ca9e7b6
5 changed files with 245 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
plugins/Hartree_Fock/SCF_old.irp.f Normal file
View File

@ -0,0 +1,61 @@
program scf
BEGIN_DOC
! Produce `Hartree_Fock` MO orbital
! output: mo_basis.mo_tot_num mo_basis.mo_label mo_basis.ao_md5 mo_basis.mo_coef mo_basis.mo_occ
! output: hartree_fock.energy
! optional: mo_basis.mo_coef
END_DOC
call create_guess
call orthonormalize_mos
call run
end
subroutine create_guess
implicit none
BEGIN_DOC
! Create a MO guess if no MOs are present in the EZFIO directory
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_mo_basis_mo_coef(exists)
if (.not.exists) then
if (mo_guess_type == "HCore") then
mo_coef = ao_ortho_lowdin_coef
TOUCH mo_coef
mo_label = 'Guess'
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
SOFT_TOUCH mo_coef mo_label
else if (mo_guess_type == "Huckel") then
call huckel_guess
else
print *, 'Unrecognized MO guess type : '//mo_guess_type
stop 1
endif
endif
end
subroutine run
BEGIN_DOC
! Run SCF calculation
END_DOC
use bitmasks
implicit none
double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem
double precision :: EHF
integer :: i_it, i, j, k
EHF = HF_energy
mo_label = "Canonical"
! Choose SCF algorithm
call damping_SCF ! Deprecated routine
! call Roothaan_Hall_SCF
end

94
plugins/read_integral/print_integrals_ao.irp.f Normal file
View File

@ -0,0 +1,94 @@
program print_integrals
PROVIDE ezfio_filename
call ezfio_set_integrals_monoelec_disk_access_ao_one_integrals('None')
call ezfio_set_integrals_bielec_disk_access_ao_integrals('None')
call run
end
subroutine run
implicit none
integer :: iunit
integer :: getunitandopen
integer ::i,j,k,l
double precision :: integral
iunit = getunitandopen('kinetic_ao','w')
do i=1,ao_num
do j=1,ao_num
integral = ao_kinetic_integral(i,j)
if (dabs(integral) > ao_integrals_threshold) then
write(iunit,*) i,j, integral
endif
enddo
enddo
close(iunit)
iunit = getunitandopen('overlap_ao','w')
do i=1,ao_num
do j=1,ao_num
integral = ao_overlap(i,j)
if (dabs(integral) > ao_integrals_threshold) then
write(iunit,*) i,j, integral
endif
enddo
enddo
close(iunit)
iunit = getunitandopen('nuclear_ao','w')
do i=1,ao_num
do j=1,ao_num
integral = ao_nucl_elec_integral(i,j)
if (dabs(integral) > ao_integrals_threshold) then
write(iunit,*) i,j, integral
endif
enddo
enddo
close(iunit)
! iunit = getunitandopen('pseudo_ao','w')
! do i=1,ao_num
! do j=1,ao_num
! write(iunit,*) i,j, ao_pseudo_integral(i,j)
! enddo
! enddo
! close(iunit)
PROVIDE ao_bielec_integrals_in_map
iunit = getunitandopen('bielec_ao','w')
integer*8 :: i8
integer :: i_idx, n_elements_max, k1, n_elements
integer :: ii(8), jj(8), kk(8), ll(8)
double precision, external :: ao_bielec_integral
integer(key_kind), allocatable :: keys(:)
double precision, allocatable :: values(:)
call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
allocate(keys(n_elements_max), values(n_elements_max))
do i8=0_8,ao_integrals_map%map_size
n_elements = n_elements_max
call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
do k1=1,n_elements
call bielec_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
if ( (kk(1)>ao_num).or. &
(ii(1)>ao_num).or. &
(jj(1)>ao_num).or. &
(ll(1)>ao_num) ) then
cycle
endif
k = kk(1)
i = ii(1)
l = ll(1)
j = jj(1)
integral = values(k1)
write (iunit,'(4(I5,X),D22.15)') k,i,l,j, integral
enddo
enddo
close(iunit)
end

75
plugins/read_integral/read_integrals_ao.irp.f Normal file
View File

@ -0,0 +1,75 @@
program read_integrals
PROVIDE ezfio_filename
call ezfio_set_integrals_monoelec_disk_access_ao_one_integrals("None")
call run
end
subroutine run
use map_module
implicit none
integer :: iunit
integer :: getunitandopen
integer ::i,j,k,l
double precision :: integral
double precision, allocatable :: A(:,:)
integer :: n_integrals
integer(key_kind), allocatable :: buffer_i(:)
real(integral_kind), allocatable :: buffer_values(:)
integer(key_kind) :: key
allocate (A(ao_num,ao_num))
A = 0.d0
iunit = getunitandopen('kinetic_ao','r')
do
read (iunit,*,end=10) i,j, integral
A(i,j) = integral
A(j,i) = integral
enddo
10 continue
close(iunit)
call write_one_e_integrals('ao_kinetic_integral', A, size(A,1), size(A,2))
A = 0.d0
iunit = getunitandopen('nuclear_ao','r')
do
read (iunit,*,end=12) i,j, integral
A(i,j) = integral
A(j,i) = integral
enddo
12 continue
close(iunit)
call write_one_e_integrals('ao_ne_integral', A, size(A,1), size(A,2))
call write_one_e_integrals('ao_pseudo_integral', ao_pseudo_integral,&
size(ao_pseudo_integral,1), size(ao_pseudo_integral,2))
call ezfio_set_integrals_monoelec_disk_access_ao_one_integrals("Read")
allocate(buffer_i(ao_num**4), buffer_values(ao_num**4))
iunit = getunitandopen('bielec_ao','r')
n_integrals=0
do
read (iunit,*,end=13) i,j,k,l, integral
n_integrals += 1
call bielec_integrals_index(i, j, k, l, buffer_i(n_integrals) )
buffer_values(n_integrals) = integral
enddo
13 continue
close(iunit)
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_values)
call map_sort(ao_integrals_map)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
call ezfio_set_integrals_bielec_disk_access_ao_integrals('Read')
end

25
src/Integrals_Monoelec/kin_ao_ints.irp.f
View File

@ -1,6 +1,6 @@
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num_align,ao_num) ]
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num,ao_num) ]
implicit none
integer :: i,j,n,l
double precision :: f
@ -45,8 +45,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0
@ -59,7 +57,6 @@
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
@ -122,7 +119,7 @@
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! array of the priminitve basis kinetic integrals
@ -131,27 +128,23 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
integer :: i,j,k,l
if (read_ao_one_integrals) then
call ezfio_get_ao_basis_integral_kinetic(ao_kinetic_integral(1:ao_num, 1:ao_num))
call ezfio_set_ao_basis_integral_kinetic(ao_kinetic_integral(1:ao_num, 1:ao_num))
call read_one_e_integrals('ao_kinetic_integral', ao_kinetic_integral,&
size(ao_kinetic_integral,1), size(ao_kinetic_integral,2))
print *, 'AO kinetic integrals read from disk'
else
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
!$OMP SHARED(ao_num, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
do j = 1, ao_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do i = 1, ao_num
ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) )
enddo
do i = ao_num +1,ao_num_align
ao_kinetic_integral(i,j) = 0.d0
enddo
enddo
!$OMP END PARALLEL DO
endif
if (write_ao_one_integrals) then
call ezfio_set_ao_basis_integral_kinetic(ao_kinetic_integral(1:ao_num, 1:ao_num))
call write_one_e_integrals('ao_kinetic_integral', ao_kinetic_integral,&
size(ao_kinetic_integral,1), size(ao_kinetic_integral,2))
print *, 'AO kinetic integrals written to disk'
endif
END_PROVIDER

10
src/Integrals_Monoelec/pot_ao_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)]
BEGIN_DOC
! interaction nuclear electron
END_DOC
@ -11,7 +11,8 @@ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
if (read_ao_one_integrals) then
call ezfio_get_ao_basis_integral_nuclear(ao_nucl_elec_integral(1:ao_num, 1:ao_num))
call read_one_e_integrals('ao_ne_integral', ao_nucl_elec_integral, &
size(ao_nucl_elec_integral,1), size(ao_nucl_elec_integral,2))
print *, 'AO N-e integrals read from disk'
else
@ -73,14 +74,15 @@ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
!$OMP END PARALLEL
endif
if (write_ao_one_integrals) then
call ezfio_set_ao_basis_integral_nuclear(ao_nucl_elec_integral(1:ao_num, 1:ao_num))
call write_one_e_integrals('ao_ne_integral', ao_nucl_elec_integral, &
size(ao_nucl_elec_integral,1), size(ao_nucl_elec_integral,2))
print *, 'AO N-e integrals written to disk'
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num,ao_num,nucl_num)]
BEGIN_DOC
! ao_nucl_elec_integral_per_atom(i,j,k) = -<AO(i)|1/|r-Rk|AO(j)>
! where Rk is the geometry of the kth atom