From 6b6ca9e7b6d55c355fee9e55d59e03b0f7c81e6b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 25 Sep 2017 17:58:23 +0200 Subject: [PATCH] print_integrals_ao --- plugins/Hartree_Fock/SCF_old.irp.f | 61 ++++++++++++ .../read_integral/print_integrals_ao.irp.f | 94 +++++++++++++++++++ plugins/read_integral/read_integrals_ao.irp.f | 75 +++++++++++++++ src/Integrals_Monoelec/kin_ao_ints.irp.f | 25 ++--- src/Integrals_Monoelec/pot_ao_ints.irp.f | 10 +- 5 files changed, 245 insertions(+), 20 deletions(-) create mode 100644 plugins/Hartree_Fock/SCF_old.irp.f create mode 100644 plugins/read_integral/print_integrals_ao.irp.f create mode 100644 plugins/read_integral/read_integrals_ao.irp.f diff --git a/plugins/Hartree_Fock/SCF_old.irp.f b/plugins/Hartree_Fock/SCF_old.irp.f new file mode 100644 index 00000000..03d9a91d --- /dev/null +++ b/plugins/Hartree_Fock/SCF_old.irp.f @@ -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 + + diff --git a/plugins/read_integral/print_integrals_ao.irp.f b/plugins/read_integral/print_integrals_ao.irp.f new file mode 100644 index 00000000..3f489ba8 --- /dev/null +++ b/plugins/read_integral/print_integrals_ao.irp.f @@ -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 diff --git a/plugins/read_integral/read_integrals_ao.irp.f b/plugins/read_integral/read_integrals_ao.irp.f new file mode 100644 index 00000000..c323c7e1 --- /dev/null +++ b/plugins/read_integral/read_integrals_ao.irp.f @@ -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 diff --git a/src/Integrals_Monoelec/kin_ao_ints.irp.f b/src/Integrals_Monoelec/kin_ao_ints.irp.f index 6cb2aa49..d6d09fbc 100644 --- a/src/Integrals_Monoelec/kin_ao_ints.irp.f +++ b/src/Integrals_Monoelec/kin_ao_ints.irp.f @@ -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 diff --git a/src/Integrals_Monoelec/pot_ao_ints.irp.f b/src/Integrals_Monoelec/pot_ao_ints.irp.f index 7116d2c7..22869c4c 100644 --- a/src/Integrals_Monoelec/pot_ao_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_ints.irp.f @@ -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) = - ! where Rk is the geometry of the kth atom