diff --git a/src/ao_one_e_ints_periodic/EZFIO.cfg b/src/ao_one_e_ints_periodic/EZFIO.cfg new file mode 100644 index 00000000..5dac91f7 --- /dev/null +++ b/src/ao_one_e_ints_periodic/EZFIO.cfg @@ -0,0 +1,72 @@ +[ao_integrals_e_n] +type: double precision +doc: Nucleus-electron integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_integrals_e_n] +type: Disk_access +doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_integrals_kinetic] +type: double precision +doc: Kinetic energy integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_integrals_kinetic] +type: Disk_access +doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_integrals_pseudo] +type: double precision +doc: Pseudopotential integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_integrals_pseudo] +type: Disk_access +doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_integrals_overlap] +type: double precision +doc: Overlap integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_integrals_overlap] +type: Disk_access +doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_one_e_integrals] +type: double precision +doc: Combined integrals in |AO| basis set +size: (ao_basis.ao_num,ao_basis.ao_num) +interface: ezfio + +[io_ao_one_e_integrals] +type: Disk_access +doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: Read + + +[num_kpts] +type: integer +doc: Number of k-points +interface: ezfio,provider,ocaml +default: 1 + + diff --git a/src/ao_one_e_ints_periodic/NEED b/src/ao_one_e_ints_periodic/NEED new file mode 100644 index 00000000..61d23b1e --- /dev/null +++ b/src/ao_one_e_ints_periodic/NEED @@ -0,0 +1,2 @@ +ao_basis +pseudo diff --git a/src/ao_one_e_ints_periodic/README.rst b/src/ao_one_e_ints_periodic/README.rst new file mode 100644 index 00000000..48ec8d01 --- /dev/null +++ b/src/ao_one_e_ints_periodic/README.rst @@ -0,0 +1,15 @@ +================== +ao_one_e_integrals +================== + +All the one-electron integrals in the periodic |AO| basis are here. + +Warning: this is incompatible with non-periodic |AOs|. + +The most important providers for usual quantum-chemistry calculation are: + +* `ao_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis +* `ao_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis +* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis + + diff --git a/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f b/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f new file mode 100644 index 00000000..27e17754 --- /dev/null +++ b/src/ao_one_e_ints_periodic/ao_one_e_ints.irp.f @@ -0,0 +1,29 @@ + BEGIN_PROVIDER [ complex*16, ao_one_e_integrals,(ao_num,ao_num)] +&BEGIN_PROVIDER [ double precision, ao_one_e_integrals_diag,(ao_num)] + implicit none + integer :: i,j,n,l + BEGIN_DOC + ! One-electron Hamiltonian in the |AO| basis. + END_DOC + + IF (read_ao_one_e_integrals) THEN + call ezfio_get_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals) + ELSE + ao_one_e_integrals = ao_integrals_n_e + ao_kinetic_integrals + + IF (DO_PSEUDO) THEN + ao_one_e_integrals += ao_pseudo_integrals + ENDIF + ENDIF + + DO j = 1, ao_num + ao_one_e_integrals_diag(j) = real(ao_one_e_integrals(j,j)) + ENDDO + + IF (write_ao_one_e_integrals) THEN + call ezfio_set_ao_one_e_ints_ao_one_e_integrals(ao_one_e_integrals) + print *, 'AO one-e integrals written to disk' + ENDIF + +END_PROVIDER + diff --git a/src/ao_one_e_ints_periodic/ao_overlap.irp.f b/src/ao_one_e_ints_periodic/ao_overlap.irp.f new file mode 100644 index 00000000..9c74091f --- /dev/null +++ b/src/ao_one_e_ints_periodic/ao_overlap.irp.f @@ -0,0 +1,136 @@ +BEGIN_PROVIDER [ complex*16, ao_overlap,(ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Overlap between atomic basis functions: +! :math:`\int \chi_i(r) \chi_j(r) dr)` + END_DOC + if (read_ao_integrals_overlap) then + call read_one_e_integrals_complex('ao_overlap', ao_overlap,& + size(ao_overlap,1), size(ao_overlap,2)) + print *, 'AO overlap integrals read from disk' + else + print *, 'complex AO overlap integrals must be provided' + endif +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Overlap between absolute value of atomic basis functions: +! :math:`\int |\chi_i(r)| |\chi_j(r)| dr)` + END_DOC + integer :: i,j + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i,j) & + !$OMP SHARED(ao_overlap_abs,ao_overlap,ao_num) + do j=1,ao_num + do i= 1,ao_num + ao_overlap_abs(i,j)= cdabs(ao_overlap(i,j)) + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, S_inv,(ao_num,ao_num) ] + implicit none + BEGIN_DOC +! Inverse of the overlap matrix + END_DOC + call get_pseudo_inverse_complex(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, S_half_inv, (AO_num,AO_num) ] + + BEGIN_DOC +! :math:`X = S^{-1/2}` obtained by SVD + END_DOC + + implicit none + + integer :: num_linear_dependencies + integer :: LDA, LDC + double precision, allocatable :: D(:) + complex*16, allocatable :: U(:,:),Vt(:,:) + integer :: info, i, j, k + double precision, parameter :: threshold_overlap_AO_eigenvalues = 1.d-6 + + LDA = size(AO_overlap,1) + LDC = size(S_half_inv,1) + + allocate( & + U(LDC,AO_num), & + Vt(LDA,AO_num), & + D(AO_num)) + + call svd_complex( & + AO_overlap,LDA, & + U,LDC, & + D, & + Vt,LDA, & + AO_num,AO_num) + + num_linear_dependencies = 0 + do i=1,AO_num + print*,D(i) + if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then + D(i) = 0.d0 + num_linear_dependencies += 1 + else + ASSERT (D(i) > 0.d0) + D(i) = 1.d0/sqrt(D(i)) + endif + do j=1,AO_num + S_half_inv(j,i) = (0.d0,0.d0) + enddo + enddo + write(*,*) 'linear dependencies',num_linear_dependencies + + do k=1,AO_num + if(D(k) /= 0.d0) then + do j=1,AO_num + do i=1,AO_num + S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + endif + enddo + + +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, S_half, (ao_num,ao_num) ] + implicit none + BEGIN_DOC + ! :math:`S^{1/2}` + END_DOC + + integer :: i,j,k + complex*16, allocatable :: U(:,:) + complex*16, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + + allocate(U(ao_num,ao_num),Vt(ao_num,ao_num),D(ao_num)) + + call svd_complex(ao_overlap,size(ao_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_num,ao_num) + + do i=1,ao_num + D(i) = dsqrt(D(i)) + do j=1,ao_num + S_half(j,i) = (0.d0,0.d0) + enddo + enddo + + do k=1,ao_num + do j=1,ao_num + do i=1,ao_num + S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + enddo + + deallocate(U,Vt,D) + +END_PROVIDER + diff --git a/src/ao_one_e_ints_periodic/aos_kpts.irp.f b/src/ao_one_e_ints_periodic/aos_kpts.irp.f new file mode 100644 index 00000000..e44d3a23 --- /dev/null +++ b/src/ao_one_e_ints_periodic/aos_kpts.irp.f @@ -0,0 +1,5 @@ + +BEGIN_PROVIDER [integer, ao_num_per_kpt ] + implicit none + ao_num_per_kpt = ao_num / num_kpts +END_PROVIDER diff --git a/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f b/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f new file mode 100644 index 00000000..ea84bd4f --- /dev/null +++ b/src/ao_one_e_ints_periodic/kin_ao_ints.irp.f @@ -0,0 +1,18 @@ +BEGIN_PROVIDER [complex*16, ao_kinetic_integrals, (ao_num,ao_num)] + implicit none + BEGIN_DOC + ! array of the priminitve basis kinetic integrals + ! \langle \chi_i |\hat{T}| \chi_j \rangle + END_DOC + + if (read_ao_integrals_kinetic) then + call read_one_e_integrals_complex('ao_kinetic_integrals', ao_kinetic_integrals,& + size(ao_kinetic_integrals,1), size(ao_kinetic_integrals,2)) + print *, 'AO kinetic integrals read from disk' + else + print *, 'complex AO kinetic integrals must be provided' + stop + endif +END_PROVIDER + + diff --git a/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f b/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f new file mode 100644 index 00000000..36b7f7ae --- /dev/null +++ b/src/ao_one_e_ints_periodic/pot_ao_ints.irp.f @@ -0,0 +1,28 @@ +BEGIN_PROVIDER [ complex*16, ao_integrals_n_e, (ao_num,ao_num)] + BEGIN_DOC + ! Nucleus-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + END_DOC + + if (read_ao_integrals_e_n) then + call read_one_e_integrals_complex('ao_ne_integral', ao_integrals_n_e, & + size(ao_integrals_n_e,1), size(ao_integrals_n_e,2)) + print *, 'AO N-e integrals read from disk' + else + print *, 'complex AO N-e integrals must be provided' + stop + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_integrals_n_e_per_atom, (ao_num,ao_num,nucl_num)] + BEGIN_DOC +! Nucleus-electron interaction in the |AO| basis set, per atom A. +! +! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle` + END_DOC + print *, 'ao_nucl_elec_integral_per_atom not implemented for k-points' + stop + +END_PROVIDER + diff --git a/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f new file mode 100644 index 00000000..243efdcb --- /dev/null +++ b/src/ao_one_e_ints_periodic/pot_ao_pseudo_ints.irp.f @@ -0,0 +1,25 @@ +BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)] + implicit none + BEGIN_DOC + ! Pseudo-potential integrals in the |AO| basis set. + END_DOC + + if (read_ao_integrals_pseudo) then + call ezfio_get_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals) + print *, 'AO pseudopotential integrals read from disk' + else + + if (do_pseudo) then + print *, irp_here, 'Not yet implemented' + stop -1 + endif + endif + + if (write_ao_integrals_pseudo) then + call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(ao_pseudo_integrals) + print *, 'AO pseudopotential integrals written to disk' + endif + +END_PROVIDER + + diff --git a/src/ao_one_e_ints_periodic/test.irp.f b/src/ao_one_e_ints_periodic/test.irp.f new file mode 100644 index 00000000..3519c17f --- /dev/null +++ b/src/ao_one_e_ints_periodic/test.irp.f @@ -0,0 +1,3 @@ +program test + print *, ' hello' +end