WIP jast deriv

Makefile

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

deriv_num.irp.f

@@ -0,0 +1,47 @@
+program jastrow
+  implicit none
+  print *, 'Nabla J1'
+  integer::k
+  double precision :: j1, j2, j0, deriv, dt, lapl
+  dt = 1.d-4
+
+BEGIN_TEMPLATE 
+lapl = 0.d0
+j0 = $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 
+
+  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
+
+enddo
+  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 ;;
+!elnuc_dist ; (1,1); 1,1 ;;
+!elec_dist  ; (1,2); 1,2 ;;
+
+
+end program

el_nuc_el.irp.f

@@ -29,7 +29,8 @@ BEGIN_PROVIDER [ double precision, factor_een ]
                    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 * (rial + rjal) * riam * rjam_cn
+                                  rijk  * rjam_cn
                 enddo
              enddo
           enddo
@@ -68,7 +69,6 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
              cn = cord_vect_lkp(l, k, p, typenuc_arr(a))
 
              do j = 1, nelec
-                factor_een_deriv_e(:, j) = 0.d0
                 rjal = rescale_een_n(j, a, l)
                 rjam_cn = rescale_een_n(j, a, m) * cn
 
@@ -83,7 +83,7 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
                    rijk = rescale_een_e(i, j, k)
 
                    do ii = 1, 4
-                      drijk(ii) = rescale_een_e_deriv_e(ii, i, j, k)
+                      drijk(ii) = rescale_een_e_deriv_e(ii, j, i, k)
                    enddo
 
                    lap = 0.0d0
@@ -92,14 +92,22 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
                    v1 = rijk * (rial + rjal)
                    v2 = rjam_cn * riam
 
-                   do ii = 1, 4
+                   do ii = 1, 3
+                      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
+!                      factor_een_deriv_e(ii, j) = factor_een_deriv_e(ii, j) + &
+!                         drijk(ii)
+                      lap = lap + d1 * d2
+                   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)
-                      lap = lap + d1 * d2
-                   enddo
 
                 enddo
              enddo
@@ -108,6 +116,4 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
     enddo
  enddo
 
- factor_een_deriv_e = 0.5d0 * factor_een_deriv_e
-
 END_PROVIDER

electrons.irp.f

@@ -48,6 +48,7 @@ BEGIN_PROVIDER [ double precision, elec_dist, (nelec, nelec) ]
        z = elec_coord(i, 3) - elec_coord(j, 3)
        elec_dist(i, j) = dsqrt( x*x + y*y + z*z )
     enddo
+!    elec_dist(j, j) = 1.d-10
  enddo
 END_PROVIDER
 

rescale.irp.f

@@ -188,7 +188,7 @@ BEGIN_PROVIDER [double precision, elec_dist_deriv_e, (4, nelec, nelec)]
  BEGIN_DOC
  ! Derivative of R_{ij} wrt x
  ! Dimensions 1-3 : dx, dy, dz
- ! Dimension 4 : d2x + d2y + d2z
+ ! Dimension 4 : d2x + d2y + d2z = 2/rij
  END_DOC
  implicit none
  integer :: i, ii, j
@@ -196,17 +196,14 @@ BEGIN_PROVIDER [double precision, elec_dist_deriv_e, (4, nelec, nelec)]
 
  do j = 1, nelec
     do i = 1, nelec
-       rij_inv = sign(1.0d0, dble(i - j)) / elec_dist(i, j)
+       rij_inv = 1.0d0 / elec_dist(i, j)
-       lap = 0.0d0
        do ii = 1, 3
           ! \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
-          ! 1 / \sqrt{c+(x-x0)^2} - (x-x0)^2 /\left(c+(x-x0)^2\right)^{3/2}
-          lap = lap + rij_inv - elec_dist_deriv_e(ii, i, j) * elec_dist_deriv_e(ii, i, j) * rij_inv
        end do
-       elec_dist_deriv_e(4, i, j) = lap
+       elec_dist_deriv_e(4, i, j) = 2.d0 * rij_inv
-       elec_dist_deriv_e(:, i, i) = 0.0d0
     end do
+    elec_dist_deriv_e(:, j, j) = 0.0d0
  end do
 
 END_PROVIDER