Jee derivatives fully functional

Makefile

@@ -1,11 +1,8 @@
 IRPF90 = irpf90  #-a -d
-FC     = gfortran -g
+FC     = gfortran
 FCFLAGS= -O2 -ffree-line-length-none -I .
-FC     = ifort -C -traceback -g
-#FCFLAGS= -O2 -ffree-line-length-none -I .
-#NINJA  = ninja
+NINJA  = ninja
 AR     = ar
-ARCHIVE= ar crs
 RANLIB = ranlib
 
 SRC=

deriv_num.irp.f

@@ -1,45 +1,45 @@
 program jastrow
   implicit none
-  print *, 'Nabla J1'
-  integer::k
+  print *, 'Nabla Jeen'
+  integer :: k
   double precision :: j1, j2, j0, deriv, dt, lapl
-  dt = 1.d-4
+  dt = 1.0d-5
 
 BEGIN_TEMPLATE 
-lapl = 0.d0
-j0 = $X $Y 
-do k=1,3
+lapl = 0.0d0
+j0 = $X $Y
 
-  elec_coord(1,k) -= dt
-  TOUCH elec_coord
-  j1 = $X $Y 
+do k = 1, 3
+   elec_coord(1, k) -= dt
+   TOUCH elec_coord
+   j1 = $X $Y 
 
-  elec_coord(1,k) += 2.d0*dt
-  TOUCH elec_coord
-  j2 = $X $Y 
+   elec_coord(1, k) += 2.0d0*dt
+   TOUCH elec_coord
+   j2 = $X $Y 
 
-  deriv = (j2 - j1)/(2.d0*dt)
-  lapl += (j2 - 2.d0*j0 + j1)/(dt*dt)
-  print *, 'deriv $X '
-  print *, deriv
-  print *, $X_deriv_e(k,$Z)
-  print *, ''
-
-  elec_coord(1,k) -= dt
-  TOUCH elec_coord
+   deriv = (j2 - j1) / (2.0d0 * dt)
+   lapl += (j2 - 2.0d0*j0 + j1) / (dt*dt)
+   print *, 'Deriv $X '
+   print *, deriv
+   print *, $X_deriv_e(k, $Z)
+   print *, ''
 
+   elec_coord(1, k) -= dt
+   TOUCH elec_coord
 enddo
-  print *, 'lapl $X '
-  print *, lapl
-  print *, $X_deriv_e(4 ,$Z)
-  print *, ''
 
-SUBST [X,Y,Z]
-factor_een ;     ;  1   ;;
+print *, 'Lapl $X '
+print *, lapl
+print *, $X_deriv_e(4, $Z)
+print *, ''
+
+SUBST [X, Y, Z]
+factor_een ; ; 1;;
 END_TEMPLATE
 !rescale_een_e ; (1,3,1) ; 1,3,1 ;;
 !rescale_een_n ; (1,1,2) ; 1,1,2 ;;
-!rescale_een_e ; (1,2,2) ; 1,2,2 ;;
+!rescale_een_e ; (1, 2, 2) ; 1, 2, 2 ;;
 !elnuc_dist ; (1,1); 1,1 ;;
 !elec_dist  ; (1,2); 1,2 ;;
 

el_nuc_el.irp.f

@@ -28,9 +28,7 @@ BEGIN_PROVIDER [ double precision, factor_een ]
                    rial = rescale_een_n(i, a, l)
                    riam = rescale_een_n(i, a, m)
                    rijk = rescale_een_e(i, j, k)
-                   factor_een = factor_een + &
-!                                  rijk * (rial + rjal) * riam * rjam_cn
-                                  rijk  * rjam_cn
+                   factor_een = factor_een + rijk * (rial + rjal) * riam * rjam_cn
                 enddo
              enddo
           enddo
@@ -48,7 +46,7 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
  END_DOC
  integer :: i, ii, j, a, p, k, l, lmax, m
  double precision :: riam, rjam_cn, rial, rjal, rijk
- double precision, dimension(4) :: driam, drjam_cn, drial, drjal, drijk, x
+ double precision, dimension(4) :: driam, drjam_cn, drial, drjal, drijk
  double precision :: cn, v1, v2, d1, d2, lap
 
  factor_een_deriv_e = 0.0d0
@@ -86,28 +84,22 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
                       drijk(ii) = rescale_een_e_deriv_e(ii, j, i, k)
                    enddo
 
-                   lap = 0.0d0
-                   x(1:3) = 0.0d0
-                   x(4) = 2.0d0
                    v1 = rijk * (rial + rjal)
                    v2 = rjam_cn * riam
 
+                   lap = 0.0d0
                    do ii = 1, 3
-                      d1 = drijk(ii) * (rial + rjal) + rijk * (rial + drjal(ii))
+                      d1 = drijk(ii) * (rial + rjal) + rijk * drjal(ii)
                       d2 = drjam_cn(ii) * riam
-                      factor_een_deriv_e(ii, j) = factor_een_deriv_e(ii, j) + &
-                         v1 * d2 + d1 * v2
-!                      factor_een_deriv_e(ii, j) = factor_een_deriv_e(ii, j) + &
-!                         drijk(ii)
                       lap = lap + d1 * d2
+                      factor_een_deriv_e(ii, j) += v1 * d2 + d1 * v2
                    enddo
 
-                   ii = 4
-                   d1 = drijk(ii) * (rial + rjal) + rijk * (rial + drjal(ii))
-                   d2 = drjam_cn(ii) * riam
-                   factor_een_deriv_e(ii, j) = factor_een_deriv_e(ii, j) + &
-                         v1 * d2 + d1 * v2 + x(ii) * lap
                    ! v(x) u''(x) + 2 * u'(x) v'(x) + u(x) v''(x)
+                   ii = 4
+                   d1 = drijk(ii) * (rial + rjal) + rijk * drjal(ii)
+                   d2 = drjam_cn(ii) * riam
+                   factor_een_deriv_e(ii, j) += v1 * d2 + d1 * v2 + 2.0d0 * lap
 
                 enddo
              enddo

rescale.irp.f

@@ -132,19 +132,16 @@ BEGIN_PROVIDER [double precision, elnuc_dist_deriv_e, (4, nelec, nnuc)]
  END_DOC
  implicit none
  integer :: i, ii, a
- double precision :: ria_inv, lap
+ double precision :: ria_inv
 
  do a = 1, nnuc 
     do i = 1, nelec
        ria_inv = 1.0d0 / elnuc_dist(i, a)
-       lap = 0.0d0
        do ii = 1, 3
           ! \frac{x-x0}{\sqrt{c+(x-x0)^2}}
           elnuc_dist_deriv_e(ii, i, a) = (elec_coord(i, ii) - nuc_coord(a, ii)) * ria_inv
-          ! 1 / \sqrt{c+(x-x0)^2} - (x-x0)^2 /\left(c+(x-x0)^2\right)^{3/2}
-          lap = lap + ria_inv - elnuc_dist_deriv_e(ii, i, a) * elnuc_dist_deriv_e(ii, i, a) * ria_inv
        end do
-       elnuc_dist_deriv_e(4, i, a) = lap
+       elnuc_dist_deriv_e(4, i, a) = 2.0d0 * ria_inv
     end do
  end do
 END_PROVIDER
@@ -192,7 +189,7 @@ BEGIN_PROVIDER [double precision, elec_dist_deriv_e, (4, nelec, nelec)]
  END_DOC
  implicit none
  integer :: i, ii, j
- double precision :: rij_inv, lap
+ double precision :: rij_inv
 
  do j = 1, nelec
     do i = 1, nelec
@@ -201,7 +198,7 @@ BEGIN_PROVIDER [double precision, elec_dist_deriv_e, (4, nelec, nelec)]
           ! \frac{x-x0}{\sqrt{c+(x-x0)^2}}
           elec_dist_deriv_e(ii, i, j) = (elec_coord(i, ii) - elec_coord(j, ii)) * rij_inv
        end do
-       elec_dist_deriv_e(4, i, j) = 2.d0 * rij_inv
+       elec_dist_deriv_e(4, i, j) = 2.0d0 * rij_inv
     end do
     elec_dist_deriv_e(:, j, j) = 0.0d0
  end do