Reading input file. BROKENgit add .!

Makefile

@@ -1,10 +1,8 @@
-IRPF90 = ~/irpf90/bin/irpf90 --codelet factor_een:200
+IRPF90 = irpf90  #-a -d
-#FC = gfortran
+FC     = gfortran
-#FCFLAGS= -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+FCFLAGS= -O2 -ffree-line-length-none -I .
-FC     = ifort -g
-FCFLAGS= -O2 -xHost -I .
 NINJA  = ninja
-AR     = ar crs
+AR     = ar
 RANLIB = ranlib
 
 SRC=

OUT_GF0

@@ -0,0 +1,5 @@
+  0.17703281415941130     
+  0.17332587478375222     
+  0.51043161486617861     
+ The total Jastrow factor
+   2.3650290467901534     

OUT_REF

@@ -0,0 +1,5 @@
+  0.278125485678554     
+  0.308815718675770     
+  0.533040790952379     
+ The total Jastrow factor
+   3.06479902202983     

el_nuc_el.irp.f

@@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, factor_een]
  BEGIN_DOC
  ! Electron-electron nucleus contribution to Jastrow factor
  END_DOC
- integer :: i, j, alpha, p, k, l, lmax = 0
+ integer :: i, j, alpha, p, k, l, lmax
  double precision :: x, y, z, t, c_inv, u, a, b, a2, b2, c, t0
  PROVIDE cord_vect
  factor_een = 0.0d0
@@ -12,14 +12,14 @@ BEGIN_PROVIDER [double precision, factor_een]
     do j = 1, nelec
        b = rescale_een_n(j, alpha)
        do i = 1, nelec
-        u = rescale_een_e(i,j)
+        u = rescale_een_e(i, j)
         a = rescale_een_n(i, alpha)
-        a2 = a*a
+        a2 = a * a
-        b2 = b*b
+        b2 = b * b
         c = rescale_een_n(i, alpha) * rescale_een_n(j, alpha)
-        c_inv = 1.d0/c
+        c_inv = 1.0d0 / c
         do p = 2, ncord
-           x = 1.d0
+           x = 1.0d0
            do k = 0, p - 1
               if ( k /= 0 ) then
                  lmax = p - k
@@ -27,23 +27,24 @@ BEGIN_PROVIDER [double precision, factor_een]
                  lmax = p - k - 2
               end if
               t = x
-              do l=1,rshift(p - k,1)
+              do l = 1, rshift(p - k, 1)
-                t = t*c
+                t = t * c
               end do
               ! We have suppressed this from the following loop:
               ! if ( iand(p - k - l, 1) == 0 ) then
               !
               ! Start from l=0 when p-k is even
               ! Start from l=1 when p-k is odd
-              if (iand(p-k,1) == 0) then
+              if (iand(p - k, 1) == 0) then
-                y = 1.d0
+                y = 1.0d0
-                z = 1.d0
+                z = 1.0d0
               else
                 y = a
                 z = b
               endif
-              do l = iand(p-k,1), lmax, 2
+              do l = iand(p - k, 1), lmax, 2
-                factor_een = factor_een + cord_vect(l, k, p, alpha) * (y+z) * t
+                ! factor_een = factor_een + cord_vect(l, k, p, alpha) * (y + z) * t
+                factor_een = factor_een + cord_vect(l + k + p + alpha) * (y + z) * t
                 t = t * c_inv
                 y = y * a2
                 z = z * b2

electrons.irp.f

@@ -46,17 +46,21 @@ BEGIN_PROVIDER [double precision, factor_ee]
  ! Electron-electron contribution to Jastrow factor
  END_DOC
  integer :: i, j, p
- double precision :: pow_ser = 0.0d0
+ double precision :: pow_ser, x
  factor_ee = 0.0d0
+ pow_ser = 0.0d0
 
  do j = 1 , nelec
     do i = 1, nelec
+       x = rescale_ee(i, j) 
        do p = 2, nbord
-          pow_ser = pow_ser + bord_vect(p) * rescale_ee(i, j) ** p
+          x = x * rescale_ee(i, j)
+          pow_ser = pow_ser + bord_vect(p) * x
        end do
        factor_ee = factor_ee + bord_vect(1) * rescale_ee(i, j) &
             / (1 + bord_vect(2) * rescale_ee(i, j)) + pow_ser
     end do
  end do
+
  factor_ee = 0.5d0 * factor_ee
 END_PROVIDER

jastrow_provider.irp.f

@@ -3,6 +3,7 @@ BEGIN_PROVIDER [ double precision, jastrow_full ]
  BEGIN_DOC
  ! Complete jastrow factor 
  END_DOC
+
  print *, factor_ee
  print *, factor_en
  print *, factor_een

nuclei.irp.f

@@ -46,13 +46,17 @@ BEGIN_PROVIDER [double precision, factor_en]
  ! Electron-nuclei contribution to Jastrow factor
  END_DOC
  integer :: i, j, p
- double precision :: pow_ser = 0.0d0
+ double precision :: pow_ser, x
+
  factor_en = 0.0d0
+ pow_ser = 0.0d0
 
  do j = 1 , nnuc
     do i = 1, nnuc
+       x = rescale_en(i, j)
        do p = 2, naord
-          pow_ser = pow_ser + aord_vect(p) * rescale_en(i, j) ** p
+          x = x * rescale_en(i, j) 
+          pow_ser = pow_ser + aord_vect(p) * x
        end do
        factor_en = factor_en + aord_vect(1) * rescale_en(i, j) &
             / (1 + aord_vect(2) * rescale_en(i, j)) + pow_ser

orders.irp.f

@@ -23,36 +23,64 @@ BEGIN_PROVIDER [integer, ncord]
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, aord_vect, (naord)]
- implicit none
+&BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
- BEGIN_DOC
+&BEGIN_PROVIDER [double precision, cord_vect, (ncord * ncord * ncord * nnuc)]
- ! Vector of the `a' coefficients
+implicit none
- END_DOC
+PROVIDE naord
- integer :: i
+PROVIDE nbord
- PROVIDE seed
+PROVIDE ncord
- call random_number(aord_vect)
+BEGIN_DOC
- aord_vect = aord_vect*.1d-2
+! Read Jastow coefficients from file (NEEDS OPTIMIZATION!)
- FREE seed
+END_DOC
+character(len=*), parameter :: FILE_NAME = "orders_inp"
+integer :: i, fu, rc
+double precision, dimension(naord + nbord + ncord * ncord * ncord * nnuc) :: allord_vect
+
+open(action='read', file=FILE_NAME, iostat=rc, newunit=fu)
+
+do i = 1, naord + nbord + ncord * ncord * ncord * nnuc
+   read(fu, *) allord_vect(i)
+end do
+
+aord_vect = allord_vect(1:naord)
+bord_vect = allord_vect(naord + 1: naord + nbord)
+cord_vect = allord_vect(naord + nbord + 1:)
+
+close(fu)
+
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
+! BEGIN_PROVIDER [double precision, aord_vect, (naord)]
- implicit none
+!  implicit none
- BEGIN_DOC
+!  BEGIN_DOC
- ! Vector of the `b' coefficients
+!  ! Vector of the `a' coefficients
- END_DOC
+!  END_DOC
- integer :: i
+!  integer :: i
- PROVIDE seed
+!  PROVIDE seed
- call random_number(bord_vect)
+!  call random_number(aord_vect)
- bord_vect = bord_vect*.1d-6
+!  aord_vect = aord_vect*.1d-2
- FREE seed
+!  FREE seed
-END_PROVIDER
+! END_PROVIDER
-
+! 
-BEGIN_PROVIDER [double precision, cord_vect, (0:ncord,0:ncord,ncord,nnuc)]
+! BEGIN_PROVIDER [double precision, bord_vect, (nbord)]
- implicit none
+!  implicit none
- BEGIN_DOC
+!  BEGIN_DOC
- ! Vector of the `c' coefficients
+!  ! Vector of the `b' coefficients
- END_DOC
+!  END_DOC
- PROVIDE seed
+!  integer :: i
- call random_number(cord_vect)
+!  PROVIDE seed
- cord_vect = cord_vect*.1d-4
+!  call random_number(bord_vect)
- FREE seed
+!  bord_vect = bord_vect*.1d-6
-END_PROVIDER
+!  FREE seed
+! END_PROVIDER
+! 
+! BEGIN_PROVIDER [double precision, cord_vect, (0:ncord,0:ncord,ncord,nnuc)]
+!  implicit none
+!  BEGIN_DOC
+!  ! Vector of the `c' coefficients
+!  END_DOC
+!  PROVIDE seed
+!  call random_number(cord_vect)
+!  cord_vect = cord_vect*.1d-4
+!  FREE seed
+! END_PROVIDER