From d4dc02363e1a89f71354c7e8050f347177ba3999 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 14 Sep 2022 18:25:11 +0200 Subject: [PATCH] Introducing slowly shells in basis --- src/ao_basis/EZFIO.cfg | 5 +- src/ao_basis/aos.irp.f | 59 ++++++++------- src/basis/EZFIO.cfg | 1 - src/basis/basis.irp.f | 159 +++++++++++++++++++++++++---------------- 4 files changed, 134 insertions(+), 90 deletions(-) diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index 51d726da..2099ad59 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -36,13 +36,13 @@ interface: ezfio, provider type: double precision doc: Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** AOs. size: (ao_basis.ao_num,ao_basis.ao_prim_num_max) -interface: ezfio, provider +interface: ezfio [ao_expo] type: double precision doc: Exponents for each primitive of each |AO| size: (ao_basis.ao_num,ao_basis.ao_prim_num_max) -interface: ezfio, provider +interface: ezfio [ao_md5] type: character*(32) @@ -67,3 +67,4 @@ doc: Use normalized primitive functions interface: ezfio, provider default: true + diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 1cbd3976..9d8cf018 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -1,11 +1,3 @@ -BEGIN_PROVIDER [ integer, ao_prim_num_max ] - implicit none - BEGIN_DOC - ! Max number of primitives. - END_DOC - ao_prim_num_max = maxval(ao_prim_num) -END_PROVIDER - BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ] implicit none BEGIN_DOC @@ -23,6 +15,32 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ] END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_coef , (ao_num,ao_prim_num_max) ] +&BEGIN_PROVIDER [ double precision, ao_expo , (ao_num,ao_prim_num_max) ] + implicit none + BEGIN_DOC +! Primitive coefficients and exponents for each atomic orbital. Copied from shell info. + END_DOC + + integer :: i, l + do i=1,ao_num + l = ao_shell(i) + ao_coef(i,:) = shell_coef(l,:) + ao_expo(i,:) = shell_expo(l,:) + end do + +END_PROVIDER + + +BEGIN_PROVIDER [ integer, ao_prim_num_max ] + implicit none + BEGIN_DOC + ! Max number of primitives. + END_DOC + ao_prim_num_max = shell_prim_num_max +END_PROVIDER + BEGIN_PROVIDER [ integer, ao_first_of_shell, (shell_num) ] implicit none BEGIN_DOC @@ -44,20 +62,20 @@ END_PROVIDER BEGIN_DOC ! Coefficients including the |AO| normalization END_DOC + + do i=1,ao_num + l = ao_shell(i) + ao_coef_normalized(i,:) = shell_coef(l,:) * shell_normalization_factor(l) + end do + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c integer :: l, powA(3), nz integer :: i,j,k nz=100 - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 - ao_coef_normalized = 0.d0 + C_A = 0.d0 do i=1,ao_num -! powA(1) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) -! powA(2) = 0 -! powA(3) = 0 powA(1) = ao_power(i,1) powA(2) = ao_power(i,2) powA(3) = ao_power(i,3) @@ -67,18 +85,9 @@ END_PROVIDER do j=1,ao_prim_num(i) call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j), & powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) - ao_coef_normalized(i,j) = ao_coef(i,j)/dsqrt(norm) - enddo - else - do j=1,ao_prim_num(i) - ao_coef_normalized(i,j) = ao_coef(i,j) + ao_coef_normalized(i,j) = ao_coef_normalized(i,j)/dsqrt(norm) enddo endif - - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) - ! Normalization of the contracted basis functions if (ao_normalized) then norm = 0.d0 diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index a6864418..342fd4cc 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -72,4 +72,3 @@ doc: Exponents in the shell size: (basis.prim_num) interface: ezfio, provider - diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index b750d75a..0fe84506 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -1,67 +1,11 @@ -BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] - implicit none - BEGIN_DOC - ! Number of primitives per |AO| - END_DOC - - logical :: has - PROVIDE ezfio_filename - if (mpi_master) then - if (size(shell_normalization_factor) == 0) return - - call ezfio_has_basis_shell_normalization_factor(has) - if (has) then - write(6,'(A)') '.. >>>>> [ IO READ: shell_normalization_factor ] <<<<< ..' - call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor) - else - - double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c - integer :: l, powA(3), nz - integer :: i,j,k - nz=100 - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 - - do i=1,shell_num - - powA(1) = shell_ang_mom(i) - powA(2) = 0 - powA(3) = 0 - - norm = 0.d0 - do k=1, prim_num - if (shell_index(k) /= i) cycle - do j=1, prim_num - if (shell_index(j) /= i) cycle - call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), & - powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*prim_coef(j)*prim_coef(k) * prim_normalization_factor(j) * prim_normalization_factor(k) - enddo - enddo - shell_normalization_factor(i) = 1.d0/dsqrt(norm) - enddo - - endif - endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( shell_normalization_factor, (shell_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read shell_normalization_factor with MPI' - endif - IRP_ENDIF - - call write_time(6) - +BEGIN_PROVIDER [ integer, shell_prim_num_max ] + implicit none + BEGIN_DOC + ! Max number of primitives. + END_DOC + shell_prim_num_max = maxval(shell_prim_num) END_PROVIDER - BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] implicit none BEGIN_DOC @@ -120,3 +64,94 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] call write_time(6) END_PROVIDER + + BEGIN_PROVIDER [ double precision, shell_coef , (shell_num, shell_prim_num_max) ] +&BEGIN_PROVIDER [ double precision, shell_expo , (shell_num, shell_prim_num_max) ] + implicit none + BEGIN_DOC +! Primitive coefficients and exponents for each shell. + END_DOC + + integer :: i, idx + integer :: count(shell_num) + + count(:) = 0 + do i=1, prim_num + idx = shell_index(i) + count(idx) += 1 + shell_coef(idx, count(idx)) = prim_coef(i) + shell_expo(idx, count(idx)) = prim_expo(i) + end do +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, shell_coef_normalized, (shell_num,shell_prim_num_max) ] +&BEGIN_PROVIDER [ double precision, shell_normalization_factor, (shell_num) ] + implicit none + BEGIN_DOC + ! Coefficients including the |shell| normalization + END_DOC + logical :: has + PROVIDE ezfio_filename + + shell_normalization_factor(:) = 1.d0 + if (mpi_master) then + if (size(shell_normalization_factor) == 0) return + + call ezfio_has_basis_shell_normalization_factor(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: shell_normalization_factor ] <<<<< ..' + call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor) + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( shell_normalization_factor, (shell_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read shell_normalization_factor with MPI' + endif + IRP_ENDIF + + call write_time(6) + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c + integer :: l, powA(3), nz + integer :: i,j,k + nz=100 + C_A = 0.d0 + powA = 0 + shell_coef_normalized = 0.d0 + + do i=1,shell_num + + powA(1) = shell_ang_mom(i) + + ! Normalization of the primitives + if (primitives_normalized) then + do j=1,shell_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,shell_expo(i,j),shell_expo(i,j), & + powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + shell_coef_normalized(i,j) = shell_coef(i,j)/dsqrt(norm) + enddo + else + do j=1,shell_prim_num(i) + shell_coef_normalized(i,j) = shell_coef(i,j) + enddo + endif + + ! Normalization of the contracted basis functions + norm = 0.d0 + do j=1,shell_prim_num(i) + do k=1,shell_prim_num(i) + call overlap_gaussian_xyz(C_A,C_A,shell_expo(i,j),shell_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*shell_coef_normalized(i,j)*shell_coef_normalized(i,k) + enddo + enddo + shell_normalization_factor(i) *= 1.d0/dsqrt(norm) + enddo +END_PROVIDER +