Jen gradient. Not laplacian

deriv_num.irp.f

@@ -1,9 +1,9 @@
 program jastrow
   implicit none
-  print *, 'Nabla Jeen'
+  print *, 'Derivatives test'
   integer :: k
   double precision :: j1, j2, j0, deriv, dt, lapl
-  dt = 1.0d-5
+  dt = 1.0d-4
 
 BEGIN_TEMPLATE 
 lapl = 0.0d0
@@ -35,11 +35,14 @@ print *, $X_deriv_e(4, $Z)
 print *, ''
 
 SUBST [X, Y, Z]
-factor_een ; ; 1;;
+factor_en ; ; 1;;
 END_TEMPLATE
+!factor_een ; ; 1;;
 !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 ;;
+!factor_en ; ; 1;;
+!rescale_en ; (1, 2) ; 1, 2 ;;
 !elnuc_dist ; (1,1); 1,1 ;;
 !elec_dist  ; (1,2); 1,2 ;;
 

el_nuc_el.irp.f

@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
  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
- double precision :: cn, v1, v2, d1, d2, lap
+ double precision :: cn, v1, v2, d1, d2, lap1, lap2
 
  factor_een_deriv_e = 0.0d0
 
@@ -84,22 +84,24 @@ BEGIN_PROVIDER [ double precision, factor_een_deriv_e, (4, nelec) ]
                       drijk(ii) = rescale_een_e_deriv_e(ii, j, i, k)
                    enddo
 
-                   v1 = rijk * (rial + rjal)
+                   v1 = rijk * (rial + rjal) ! v(x)
-                   v2 = rjam_cn * riam
+                   v2 = rjam_cn * riam ! u(x)
 
-                   lap = 0.0d0
+                   lap1 = 0.0d0
+                   lap2 = 0.0d0
                    do ii = 1, 3
                       d1 = drijk(ii) * (rial + rjal) + rijk * drjal(ii)
                       d2 = drjam_cn(ii) * riam
-                      lap = lap + d1 * d2
+                      lap1 = lap1 + d1 * d2
+                      lap2 = lap2 + drijk(ii) * drjal(ii)
                       factor_een_deriv_e(ii, j) += v1 * d2 + d1 * v2
                    enddo
 
                    ! v(x) u''(x) + 2 * u'(x) v'(x) + u(x) v''(x)
                    ii = 4
-                   d1 = drijk(ii) * (rial + rjal) + rijk * drjal(ii)
+                   d1 = drijk(ii) * (rial + rjal) + rijk * drjal(ii) + 2.0d0 * lap2
                    d2 = drjam_cn(ii) * riam
-                   factor_een_deriv_e(ii, j) += v1 * d2 + d1 * v2 + 2.0d0 * lap
+                   factor_een_deriv_e(ii, j) += v1 * d2 + d1 * v2 + 2.0d0 * lap1
 
                 enddo
              enddo

nuclei.irp.f

@@ -58,23 +58,76 @@ BEGIN_PROVIDER [double precision, factor_en]
  BEGIN_DOC
  ! Electron-nuclei contribution to Jastrow factor
  END_DOC
- integer :: i, j, p, q
+ integer :: i, a, p, q
  double precision :: pow_ser, x
 
  factor_en = 0.0d0
 
- do j = 1 , nnuc
+ do a = 1 , nnuc
     do i = 1, nelec
-       x = rescale_en(i, j) 
+       x = rescale_en(i, a) 
        pow_ser = 0.0d0
 
        do p = 2, naord
-          x = x * rescale_en(i, j) 
+          x = x * rescale_en(i, a) 
-          pow_ser = pow_ser + aord_vect(p + 1, typenuc_arr(j)) * x
+          pow_ser = pow_ser + aord_vect(p + 1, typenuc_arr(a)) * x
        end do
 
-       factor_en = factor_en + aord_vect(1, typenuc_arr(j)) * rescale_en(i, j) &
+       factor_en = factor_en + aord_vect(1, typenuc_arr(a)) * rescale_en(i, a) &
-            / (1 + aord_vect(2, typenuc_arr(j)) * rescale_en(i, j)) + pow_ser
+            / (1.0d0 + aord_vect(2, typenuc_arr(a)) * rescale_en(i, a)) + pow_ser
+
+    end do
+ end do
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, factor_en_deriv_e, (4, nelec) ]
+ implicit none
+ BEGIN_DOC
+ ! Dimensions 1-3 : dx, dy, dz
+ ! Dimension 4 : d2x + d2y + d2z
+ END_DOC
+ integer :: i, ii, a, p, q
+ double precision :: x, x_inv, y, den, invden, lap1, lap2
+ double precision, dimension(4) :: dx, pow_ser_g
+
+ factor_en_deriv_e = 0.0d0
+
+ do a = 1 , nnuc
+    do i = 1, nelec
+       pow_ser_g = 0.0d0
+       den = 1.0d0 + aord_vect(2, typenuc_arr(a)) * rescale_en(i, a)
+       invden = 1.0d0 / den
+
+       do ii = 1, 4
+          dx(ii) = rescale_en_deriv_e(ii, i, a)
+       enddo
+
+       lap1 = 0.0d0
+       lap2 = 0.0d0
+       do ii = 1, 3
+          x = rescale_en(i, a)
+          x_inv = 1.0d0 / x
+          do p = 2, naord
+             pow_ser_g(ii) += p * aord_vect(p + 1, typenuc_arr(a)) * x * dx(ii)
+             pow_ser_g(4) += p * (p - 1) * aord_vect(p + 1, typenuc_arr(a)) * x * x_inv * dx(ii) * dx(ii)
+             lap2 += p * aord_vect(p + 1, typenuc_arr(a)) * x
+             x = x * rescale_en(i, a) 
+          end do
+
+          ! (a1 (-2 a2 r'[i,a]^2+(1+a2 r[i,a]) r''[i,a]))/(1+a2 r[i,a])^3
+          lap1 += -2.0d0 * aord_vect(2, typenuc_arr(a)) * dx(ii) * dx(ii)
+
+          ! \frac{\text{a1} r'(i,a)}{(\text{a2} r(i,a)+1)^2}
+          factor_en_deriv_e(ii, i) += aord_vect(1, typenuc_arr(a)) &
+               * dx(ii) * invden * invden + pow_ser_g(ii)
+       enddo
+
+       ii = 4
+       lap1 += den * dx(ii)
+       lap1 = lap1 * aord_vect(1, typenuc_arr(a)) * invden * invden * invden
+       pow_ser_g(ii) += lap1 + lap2 * dx(ii)
+       factor_en_deriv_e(ii, i) += pow_ser_g(ii)
 
     end do
  end do

rescale.irp.f

@@ -33,15 +33,42 @@ BEGIN_PROVIDER [ double precision, rescale_en, (nelec, nnuc) ]
  BEGIN_DOC
  ! R = (1 - exp(-kappa r))/kappa for electron-nucleus for $J_{en}$
  END_DOC
- integer :: i, j
+ integer :: i, a
 
- do j = 1, nnuc
+ do a = 1, nnuc
     do i = 1, nelec
-       rescale_en(i, j) = (1.d0 - dexp(-kappa * elnuc_dist(i, j))) * kappa_inv
+       rescale_en(i, a) = (1.0d0 - dexp(-kappa * elnuc_dist(i, a))) * kappa_inv
     enddo
  enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, rescale_en_deriv_e, (4, nelec, nnuc) ]
+ implicit none
+ BEGIN_DOC
+ ! R = (1 - exp(-kappa r))/kappa derived wrt x
+ ! Dimensions 1-3 : dx, dy, dz
+ ! Dimension 4 : d2x + d2y + d2z
+ END_DOC
+ integer :: i, ii, a
+
+ do a = 1, nnuc
+    do i = 1, nelec
+       do ii = 1, 4
+          rescale_en_deriv_e(ii, i, a) = elnuc_dist_deriv_e(ii, i, a)
+       end do
+       rescale_en_deriv_e(4, i, a) = rescale_en_deriv_e(4, i, a) + &
+       (-kappa * rescale_en_deriv_e(1, i, a) * rescale_en_deriv_e(1, i, a)) + &
+       (-kappa * rescale_en_deriv_e(2, i, a) * rescale_en_deriv_e(2, i, a)) + &
+       (-kappa * rescale_en_deriv_e(3, i, a) * rescale_en_deriv_e(3, i, a))
+       do ii = 1, 4
+          rescale_en_deriv_e(ii, i, a) = rescale_en_deriv_e(ii, i, a) &
+               * dexp(-kappa * elnuc_dist(i, a))
+       enddo
+    enddo
+ enddo
+
+END_PROVIDER
+
 BEGIN_PROVIDER [double precision, rescale_een_e, (nelec, nelec, 0:ncord)]
  implicit none
  BEGIN_DOC