Introducing slowly shells in basis

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
 
+

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 = 0.d0
-  C_A(2) = 0.d0
-  C_A(3) = 0.d0
-  ao_coef_normalized = 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)
+        ao_coef_normalized(i,j) = ao_coef_normalized(i,j)/dsqrt(norm)
-      enddo
-    else
-      do j=1,ao_prim_num(i)
-        ao_coef_normalized(i,j) = ao_coef(i,j)
       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

src/basis/EZFIO.cfg

@@ -72,4 +72,3 @@ doc: Exponents in the shell
 size: (basis.prim_num)
 interface: ezfio, provider
 
-

src/basis/basis.irp.f

@@ -1,67 +1,11 @@
-BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
+BEGIN_PROVIDER [ integer, shell_prim_num_max ]
-  implicit none
+ implicit none
-  BEGIN_DOC
+ BEGIN_DOC
-  ! Number of primitives per |AO|
+ ! Max number of primitives.
-  END_DOC
+ END_DOC
-
+ shell_prim_num_max = maxval(shell_prim_num)
-  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)
-
 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
+