Integral over contracted basis functions OK.

Makefile

@@ -1,10 +1,10 @@
-IRPF90 = irpf90  #-a -d
+IRPF90 = irpf90  -a #-d
 FC     = ifort 
 FCFLAGS= -O3 -xP
 
 SRC=
 OBJ=
-LIB=
+LIB=-lqcio
 
 include irpf90.make
 

ao.irp.f

@@ -0,0 +1,190 @@
+BEGIN_PROVIDER [ integer, ao_num ]
+ implicit none
+
+ BEGIN_DOC
+! Number of atomic orbitals
+ END_DOC
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_num_contr(ao_num)
+!$OMP END CRITICAL (qcio_critical)
+ assert (ao_num > 0)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_prim_num, (ao_num) ]
+ implicit none
+
+ BEGIN_DOC
+! Number of primitives per atomic orbital
+ END_DOC
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_num_prim(ao_prim_num)
+!$OMP END CRITICAL (qcio_critical)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_nucl, (ao_num) ]
+ implicit none
+
+ BEGIN_DOC
+! Nucleus on which the atomic orbital is centered
+ END_DOC
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_atom(ao_nucl)
+!$OMP END CRITICAL (qcio_critical)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_power, (ao_num,3) ]
+ implicit none
+
+ BEGIN_DOC
+! x,y,z powers of the atomic orbital
+ END_DOC
+ integer :: buffer(3,ao_num)
+ integer :: i,j
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_power(buffer)
+!$OMP END CRITICAL (qcio_critical)
+
+ do i=1,3
+  do j=1,ao_num
+   ao_power(j,i) = buffer(i,j)
+  enddo
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer , ao_power_max ]
+ BEGIN_DOC
+! Maximum power among x, y and z
+ END_DOC
+ ao_power_max = maxval(ao_power_max_nucl)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer , ao_power_max_nucl, (nucl_num,3) ]
+ implicit none
+ BEGIN_DOC
+! Maximum powers of x, y and z per nucleus
+ END_DOC
+ integer :: i, j
+ do j=1,3
+  do i=1,nucl_num
+   ao_power_max_nucl(i,j) = 0
+  enddo
+ enddo
+
+ integer :: inucl
+ do j=1,3
+  do i=1,ao_num
+   inucl = ao_nucl(i)
+   ao_power_max_nucl(inucl,j) = max(ao_power(i,j),ao_power_max_nucl(inucl,j))
+  enddo
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, ao_prim_num_max ]
+ implicit none
+
+ BEGIN_DOC
+! Max Number of primitives per atomic orbital
+ END_DOC
+
+ ao_prim_num_max = maxval(ao_prim_num)
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ real, ao_expo, (ao_prim_num_max,ao_num) ]
+&BEGIN_PROVIDER [ real, ao_coef, (ao_prim_num_max,ao_num) ]
+ implicit none
+ BEGIN_DOC
+! Exponents and coefficients of the atomic orbitals
+ END_DOC
+
+ double precision :: buffer(ao_prim_num_max,ao_num)
+ integer :: i,j
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_exponent(buffer)
+!$OMP END CRITICAL (qcio_critical)
+ do i=1,ao_num
+  do j=1,ao_prim_num(i)
+   ao_expo(j,i) = buffer(j,i)
+  enddo
+ enddo
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_basis_coefficient(buffer)
+!$OMP END CRITICAL (qcio_critical)
+
+ double precision :: norm, norm2
+ double precision :: overlap
+ do i=1,ao_num
+  do j=1,ao_prim_num(i)
+   norm = overlap(ao_expo(j,i),ao_expo(j,i),ao_power(i,:))
+   norm = sqrt(norm)
+   ao_coef(j,i) = buffer(j,i)/norm
+  enddo
+ enddo
+
+END_PROVIDER
+
+double precision function ddfact2(n)
+  implicit none
+  integer :: n
+  
+  ASSERT (mod(n,2) /= 0)
+
+  integer :: i
+  ddfact2 = 1.
+  do i=1,n,2
+    ddfact2 = ddfact2 * float(i)
+  enddo
+  
+end function
+
+double precision function rintgauss(n)
+  implicit none
+
+  integer :: n
+  double precision :: pi
+  pi = acos(-1.)
+
+  rintgauss = sqrt(pi)
+  if ( n == 0 ) then
+    return
+  else if ( n == 1 ) then
+    rintgauss = 0.
+  else if ( mod(n,2) == 1) then
+    rintgauss = 0.
+  else
+    double precision :: ddfact2
+    rintgauss = rintgauss/(2.**(n/2))
+    rintgauss = rintgauss * ddfact2(n-1)
+  endif
+end function
+
+double precision function overlap(gamA,gamB,nA)
+  implicit none
+
+  real    :: gamA, gamB
+  integer :: nA(3)
+
+  double precision :: gamtot
+  gamtot = gamA+gamB
+
+  overlap=1.0
+
+  integer :: l
+  double precision :: rintgauss
+  do l=1,3
+    overlap = overlap * rintgauss(nA(l)+nA(l))/ (gamtot**(0.5+float(nA(l))))
+  enddo
+
+end function
+

debug.irp.f

@@ -1,6 +1,16 @@
 program debug
- double precision a
- double precision primitive_overlap
- a = primitive_overlap(3.d0,-1.d0,3,2.d0,1.d0,4)
- print *,  a
+ double precision :: ao_overlap
+ print *,  ao_overlap(1,1)
+ integer :: i
+ do i=1,ao_prim_num(1)
+  print *,  ao_expo(i,1), ao_coef(i,1) 
+ enddo
+ print *,  ''
+ print *,  ao_overlap(1,5)
+ print *,  ao_power(5,1:3)
+ do i=1,ao_prim_num(5)
+  print *,  ao_expo(i,5), ao_coef(i,5)
+ enddo
+ print *,  ''
+ print *,  ao_overlap(5,5)
 end

integral.irp.f

@@ -1,73 +0,0 @@
-subroutine gaussian_product(a,xa,b,xb,k,p,xp)
- implicit none
-! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
-
- double precision, intent(in) :: a,b   ! Exponents
- double precision, intent(in) :: xa,xb ! Centers
- double precision, intent(out) :: p    ! New exponent
- double precision, intent(out) :: xp   ! New center
- double precision, intent(out) :: k  ! Constant
-
- double precision:: p_inv
-
- p = a+b
- xp = (a*xa+b*xb)
-
- p_inv = 1./p
-
- xp = xp*p_inv
- k = dexp(-a*b*p_inv*(xa-xb)**2)
-
-end subroutine
-
-double precision function primitive_overlap(a,xa,i,b,xb,j)
- implicit none
- include 'constants.F'
- double precision, intent(in)    :: a,b   ! Exponents
- double precision, intent(in)    :: xa,xb ! Centers
- integer, intent(in) :: i,j   ! Powers of xa and xb
-
- integer :: ii, jj, kk, ll
- double precision:: xp
- double precision:: p
- double precision :: S(0:i,0:j)
- double precision :: inv_2p, di(i), dj(j)
-
- call gaussian_product(a,xa,b,xb,S(0,0),p,xp)
- S(0,0) = S(0,0) * dsqrt(pi/p)
-
- if (i>0) then
-   S(1,0) = (xp-xa) * S(0,0) ! TODO verifier signe
- endif
- if (j>0) then
-   S(0,1) = (xp-xb) * S(0,0) ! TODO verifier signe
- endif
-
- inv_2p = 1./(2.*p)
-
- do ii=1,max(i,j)
-  di(ii) = inv_2p * dble(ii)
- enddo
-
- if (i>1) then
-  do ii=1,i-1
-   S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
-  enddo
- endif
-
- if (j>1) then
-  do jj=1,j-1
-   S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
-  enddo
- endif
-
- do jj=1,j
-  do ii=1,i
-   S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
-  enddo
- enddo
-
- primitive_overlap = S(i,j)
-
-end function
-

nuclei.irp.f

@@ -0,0 +1,95 @@
+BEGIN_PROVIDER [ integer, nucl_num ]
+ implicit none
+
+ BEGIN_DOC
+! Number of nuclei
+ END_DOC
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_geometry_num_atom(nucl_num)
+!$OMP END CRITICAL (qcio_critical)
+ assert (nucl_num > 0)
+
+END_PROVIDER
+
+BEGIN_PROVIDER  [ real, nucl_charge, (nucl_num) ]
+ implicit none
+
+ BEGIN_DOC
+! Nuclear charge
+ END_DOC
+
+ double precision :: buffer(nucl_num)
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_geometry_charge(buffer)
+!$OMP END CRITICAL (qcio_critical)
+
+ integer :: i
+ do i=1,nucl_num
+   nucl_charge(i) = buffer(i)
+ enddo
+END_PROVIDER 
+
+BEGIN_PROVIDER [ real, nucl_coord,  (nucl_num,3) ]
+ implicit none
+
+ BEGIN_DOC
+! Nuclear coordinates
+ END_DOC
+ double precision :: buffer(3,nucl_num)
+
+!$OMP CRITICAL (qcio_critical)
+ call qcio_get_geometry_coord(buffer)
+!$OMP END CRITICAL (qcio_critical)
+
+ integer :: i,j
+ do i=1,3
+  do j=1,nucl_num
+   nucl_coord(j,i) = buffer(i,j)
+  enddo
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ real, nucl_dist_2, (nucl_num,nucl_num) ]
+&BEGIN_PROVIDER [ real, nucl_dist_vec, (nucl_num,nucl_num,3) ]
+ implicit none
+ BEGIN_DOC
+! nucl_dist_2   : Nucleus-nucleus distances squared
+! nucl_dist_vec : Nucleus-nucleus distances vectors
+ END_DOC
+
+ integer :: ie1, ie2, l
+
+ do ie2 = 1,nucl_num
+  do ie1 = 1,nucl_num
+   nucl_dist_2(ie1,ie2) = 0.d0
+  enddo
+ enddo
+
+ do l=1,3
+  do ie2 = 1,nucl_num
+   do ie1 = 1,nucl_num
+    nucl_dist_vec(ie1,ie2,l) = nucl_coord(ie1,l) - nucl_coord(ie2,l)
+   enddo
+   do ie1 = 1,nucl_num
+    nucl_dist_2(ie1,ie2) = nucl_dist_2(ie1,ie2) &
+      + nucl_dist_vec(ie1,ie2,l)*nucl_dist_vec(ie1,ie2,l) 
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ real, nucl_dist, (nucl_num,nucl_num) ]
+ implicit none
+ BEGIN_DOC
+! Nucleus-nucleus distances
+ END_DOC
+ integer :: ie1, ie2
+ do ie2 = 1,nucl_num
+  do ie1 = 1,nucl_num
+   nucl_dist(ie1,ie2) = sqrt(nucl_dist_2(ie1,ie2))
+  enddo
+ enddo
+END_PROVIDER
+

overlap.irp.f

@@ -0,0 +1,141 @@
+subroutine gaussian_product(a,xa,b,xb,k,p,xp)
+ implicit none
+! e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2}
+
+ real, intent(in) :: a,b   ! Exponents
+ real, intent(in) :: xa,xb ! Centers
+ real, intent(out) :: p    ! New exponent
+ real, intent(out) :: xp   ! New center
+ double precision, intent(out) :: k  ! Constant
+
+ double precision:: p_inv
+
+ ASSERT (a>0.)
+ ASSERT (b>0.)
+
+ p = a+b
+ xp = (a*xa+b*xb)
+
+ p_inv = 1.d0/p
+
+ xp = xp*p_inv
+ k = dexp(-a*b*p_inv*(xa-xb)**2)
+
+end subroutine
+
+double precision function primitive_overlap_oneD(a,xa,i,b,xb,j)
+ implicit none
+ include 'constants.F'
+
+ real, intent(in)    :: a,b   ! Exponents
+ real, intent(in)    :: xa,xb ! Centers
+ integer, intent(in) :: i,j   ! Powers of xa and xb
+ integer :: ii, jj, kk, ll
+ real :: xp
+ real :: p
+ double precision :: S(0:i,0:j)
+ double precision :: inv_p, di(i), dj(j)
+
+ ASSERT (a>0.)
+ ASSERT (b>0.)
+ ASSERT (i>=0)
+ ASSERT (j>=0)
+
+ ! Gaussian product
+ p = a+b
+ xp = (a*xa+b*xb)
+ inv_p = 1.d0/p
+ xp = xp*inv_p
+ S(0,0) = dexp(-a*b*inv_p*(xa-xb)**2)* dsqrt(pi*inv_p)
+
+ ! Obara-Saika recursion
+
+ if (i>0) then
+   S(1,0) = (xp-xa) * S(0,0)
+ endif
+ if (j>0) then
+   S(0,1) = (xp-xb) * S(0,0)
+ endif
+
+ do ii=1,max(i,j)
+  di(ii) = 0.5d0*inv_p*dble(ii)
+ enddo
+
+ if (i>1) then
+  do ii=1,i-1
+   S(ii+1,0) = (xp-xa) * S(ii,0) + di(ii)*S(ii-1,0)
+  enddo
+ endif
+
+ if (j>1) then
+  do jj=1,j-1
+   S(0,jj+1) = (xp-xb) * S(0,jj) + di(jj)*S(0,jj-1)
+  enddo
+ endif
+
+ do jj=1,j
+  do ii=1,i
+   S(ii,jj) = (xp-xa) * S(ii-1,jj) + di(ii-1) * S(ii-2,jj) + di(jj) * S(ii-1,jj-1)
+  enddo
+ enddo
+
+ primitive_overlap_oneD = S(i,j)
+
+end function
+
+double precision function ao_overlap(i,j)
+ implicit none
+ integer, intent(in) :: i, j
+ integer :: p,q,k
+ double precision :: integral(ao_prim_num_max,ao_prim_num_max)
+
+ double precision :: primitive_overlap_oneD
+
+ ASSERT(i>0)
+ ASSERT(j>0)
+ ASSERT(i<=ao_num)
+ ASSERT(j<=ao_num)
+
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   integral(p,q) =                         &
+   primitive_overlap_oneD (                &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),1),              &
+    ao_power(i,1),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),1),              &
+    ao_power(j,1) ) *                      &
+   primitive_overlap_oneD (                &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),2),              &
+    ao_power(i,2),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),2),              &
+    ao_power(j,2) ) *                      &
+   primitive_overlap_oneD (                &
+    ao_expo(p,i),                          &
+    nucl_coord(ao_nucl(i),3),              &
+    ao_power(i,3),                         &
+    ao_expo(q,j),                          &
+    nucl_coord(ao_nucl(j),3),              &
+    ao_power(j,3) ) 
+  enddo
+ enddo
+
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   integral(p,q) = integral(p,q)*ao_coef(p,i)*ao_coef(q,j)
+  enddo
+ enddo
+
+ ao_overlap = 0.
+ do q=1,ao_prim_num(j)
+  do p=1,ao_prim_num(i)
+   ao_overlap = ao_overlap + integral(p,q)
+  enddo
+ enddo
+ 
+end function
+
+

test/QCIO_File/basis/is_cartesian

@@ -0,0 +1 @@
+T

test/QCIO_File/basis/num_contr

@@ -0,0 +1 @@
+                  35

test/QCIO_File/geometry/nuc_energy

@@ -0,0 +1 @@
+   2.032531471843000E+01

test/QCIO_File/geometry/num_atom

@@ -0,0 +1 @@
+                   3

test/QCIO_File/mo/fitcusp

@@ -0,0 +1 @@
+T

test/QCIO_File/random/number

@@ -0,0 +1 @@
+   0.980820489500000E+10

test/QCIO_File/simulation/method

@@ -0,0 +1 @@
+VMC

test/QCIO_File/simulation/num_step

@@ -0,0 +1 @@
+                 500

test/QCIO_File/simulation/sampling

@@ -0,0 +1 @@
+Langevin 

test/QCIO_File/simulation/time_step

@@ -0,0 +1 @@
+   2.000000000000000E-01

test/QCIO_File/symmetry/num_sym

@@ -0,0 +1 @@
+                   4

test/QCIO_File/system/energy

@@ -0,0 +1 @@
+  -7.614196870000001E+01

test/QCIO_File/system/num_alpha

@@ -0,0 +1 @@
+                   7

test/QCIO_File/system/num_beta

@@ -0,0 +1 @@
+                   6

test/QCIO_File/system/title

@@ -0,0 +1 @@
+titre

test/QCIO_File/walker/num_walk

@@ -0,0 +1 @@
+               1000