From 84637931140eb873f23bf0e81d648399dabc35c4 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Mon, 27 Sep 2021 11:27:49 +0200 Subject: [PATCH] Fixed bug in factor_een_e_deriv_e. Now everything checks out. --- org/qmckl_jastrow.org | 192 +++++++++++++++++------------------------- 1 file changed, 75 insertions(+), 117 deletions(-) diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index c20eaf0..663ef3d 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -3419,7 +3419,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f #+end_src *** Test - +#+name: een_e_deriv_e #+begin_src python :results output :exports none :noweb yes import numpy as np @@ -3483,15 +3483,15 @@ for l in range(0,cord_num+1): for ii in range(0,4): een_rescaled_e_deriv_e[i,ii,j,l] = een_rescaled_e_deriv_e[i,ii,j,l] * een_rescaled_e[i,j,l] -print(" een_rescaled_e_deriv_e[1, 1, 3, 1] = ",een_rescaled_e_deriv_e[0, 0, 2, 1]) -print(" een_rescaled_e_deriv_e[1, 1, 4, 1] = ",een_rescaled_e_deriv_e[0, 0, 3, 1]) -print(" een_rescaled_e_deriv_e[1, 1, 5, 1] = ",een_rescaled_e_deriv_e[0, 0, 4, 1]) -print(" een_rescaled_e_deriv_e[2, 1, 4, 2] = ",een_rescaled_e_deriv_e[1, 0, 3, 2]) -print(" een_rescaled_e_deriv_e[2, 1, 5, 2] = ",een_rescaled_e_deriv_e[1, 0, 4, 2]) -print(" een_rescaled_e_deriv_e[2, 1, 6, 2] = ",een_rescaled_e_deriv_e[1, 0, 5, 2]) +#print(" een_rescaled_e_deriv_e[1, 1, 3, 1] = ",een_rescaled_e_deriv_e[0, 0, 2, 1]) +#print(" een_rescaled_e_deriv_e[1, 1, 4, 1] = ",een_rescaled_e_deriv_e[0, 0, 3, 1]) +#print(" een_rescaled_e_deriv_e[1, 1, 5, 1] = ",een_rescaled_e_deriv_e[0, 0, 4, 1]) +#print(" een_rescaled_e_deriv_e[2, 1, 4, 2] = ",een_rescaled_e_deriv_e[1, 0, 3, 2]) +#print(" een_rescaled_e_deriv_e[2, 1, 5, 2] = ",een_rescaled_e_deriv_e[1, 0, 4, 2]) +#print(" een_rescaled_e_deriv_e[2, 1, 6, 2] = ",een_rescaled_e_deriv_e[1, 0, 5, 2]) #+end_src - #+RESULTS: + #+RESULTS: een_e_deriv_e : een_rescaled_e_deriv_e[1, 1, 3, 1] = 0.05991352796887283 : een_rescaled_e_deriv_e[1, 1, 4, 1] = 0.011714035071545248 : een_rescaled_e_deriv_e[1, 1, 5, 1] = 0.00441398875758468 @@ -5237,12 +5237,12 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: walk_num, elec_num, cord_num, nucl_num, dim_cord_vect - integer*8 , intent(in) :: lkpm_combined_index(4,dim_cord_vect) - double precision , intent(in) :: cord_vect_full(dim_cord_vect, nucl_num) - double precision , intent(in) :: een_rescaled_e(walk_num, elec_num, elec_num, 0:cord_num) - double precision , intent(in) :: een_rescaled_n(walk_num, elec_num, nucl_num, 0:cord_num) - double precision , intent(in) :: een_rescaled_e_deriv_e(walk_num, elec_num, 4, elec_num, 0:cord_num) - double precision , intent(in) :: een_rescaled_n_deriv_e(walk_num, elec_num, 4, nucl_num, 0:cord_num) + integer*8 , intent(in) :: lkpm_combined_index(dim_cord_vect, 4) + double precision , intent(in) :: cord_vect_full(nucl_num, dim_cord_vect) + double precision , intent(in) :: een_rescaled_e(0:cord_num, elec_num, elec_num, walk_num) + double precision , intent(in) :: een_rescaled_n(0:cord_num, nucl_num, elec_num, walk_num) + double precision , intent(in) :: een_rescaled_e_deriv_e(0:cord_num, elec_num, 4, elec_num, walk_num) + double precision , intent(in) :: een_rescaled_n_deriv_e(0:cord_num, nucl_num, 4, elec_num, walk_num) double precision , intent(out) :: factor_een_deriv_e(elec_num, 4, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw @@ -5280,41 +5280,41 @@ integer function qmckl_compute_factor_een_deriv_e_f(context, walk_num, elec_num, do nw =1, walk_num do n = 1, dim_cord_vect - l = lkpm_combined_index(1, n) - k = lkpm_combined_index(2, n) - p = lkpm_combined_index(3, n) - m = lkpm_combined_index(4, n) + l = lkpm_combined_index(n, 1) + k = lkpm_combined_index(n, 2) + p = lkpm_combined_index(n, 3) + m = lkpm_combined_index(n, 4) do a = 1, nucl_num - cn = cord_vect_full(n, a) + cn = cord_vect_full(a, n) do j = 1, elec_num accu = 0.0d0 accu2 = 0.0d0 daccu = 0.0d0 daccu2 = 0.0d0 do i = 1, elec_num - accu = accu + een_rescaled_e(nw, i, j, k) * & - een_rescaled_n(nw, i, a, m) - accu2 = accu2 + een_rescaled_e(nw, i, j, k) * & - een_rescaled_n(nw, i, a, m + l) - daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & - een_rescaled_n(nw, i, a, m) - daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(nw, j, 1:4, i, k) * & - een_rescaled_n(nw, i, a, m + l) + accu = accu + een_rescaled_e(k, i, j, nw) * & + een_rescaled_n(m, a, i, nw) + accu2 = accu2 + een_rescaled_e(k, i, j, nw) * & + een_rescaled_n(m + l, a, i, nw) + daccu(1:4) = daccu(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * & + een_rescaled_n(m, a, i, nw) + daccu2(1:4) = daccu2(1:4) + een_rescaled_e_deriv_e(k, j, 1:4, i, nw) * & + een_rescaled_n(m + l, a, i, nw) end do factor_een_deriv_e(j, 1:4, nw) = factor_een_deriv_e(j, 1:4, nw) + & - (accu * een_rescaled_n_deriv_e(nw, j, 1:4, a, m + l) & - + daccu(1:4) * een_rescaled_n(nw, j, a, m + l) & - + daccu2(1:4) * een_rescaled_n(nw, j, a, m) & - + accu2 * een_rescaled_n_deriv_e(nw, j, 1:4, a, m)) * cn + (accu * een_rescaled_n_deriv_e(m + l, a, 1:4, j, nw) & + + daccu(1:4) * een_rescaled_n(m + l, a, j, nw) & + + daccu2(1:4) * een_rescaled_n(m, a, j, nw) & + + accu2 * een_rescaled_n_deriv_e(m, a, 1:4, j, nw)) * cn factor_een_deriv_e(j, 4, nw) = factor_een_deriv_e(j, 4, nw) + 2.0d0 * ( & - daccu (1) * een_rescaled_n_deriv_e(nw, j, 1, a, m + l) + & - daccu (2) * een_rescaled_n_deriv_e(nw, j, 2, a, m + l) + & - daccu (3) * een_rescaled_n_deriv_e(nw, j, 3, a, m + l) + & - daccu2(1) * een_rescaled_n_deriv_e(nw, j, 1, a, m ) + & - daccu2(2) * een_rescaled_n_deriv_e(nw, j, 2, a, m ) + & - daccu2(3) * een_rescaled_n_deriv_e(nw, j, 3, a, m ) ) * cn + daccu (1) * een_rescaled_n_deriv_e(m + l, a, 1, j, nw) + & + daccu (2) * een_rescaled_n_deriv_e(m + l, a, 2, j, nw) + & + daccu (3) * een_rescaled_n_deriv_e(m + l, a, 3, j, nw) + & + daccu2(1) * een_rescaled_n_deriv_e(m, a, 1, j, nw ) + & + daccu2(2) * een_rescaled_n_deriv_e(m, a, 2, j, nw ) + & + daccu2(3) * een_rescaled_n_deriv_e(m, a, 3, j, nw ) ) * cn end do end do @@ -5407,102 +5407,60 @@ import numpy as np <> +<> + +<> + kappa = 1.0 -elec_coord = np.array(elec_coord)[0] -nucl_coord = np.array(nucl_coord) -elnuc_dist = np.zeros(shape=(elec_num, nucl_num),dtype=float) -for i in range(elec_num): - for j in range(nucl_num): - elnuc_dist[i, j] = np.linalg.norm(elec_coord[i] - nucl_coord[:,j]) +factor_een = 0.0 -elnuc_dist_deriv_e = np.zeros(shape=(4, elec_num, nucl_num),dtype=float) -for a in range(nucl_num): - for i in range(elec_num): - rij_inv = 1.0 / elnuc_dist[i, a] - for ii in range(3): - elnuc_dist_deriv_e[ii, i, a] = (elec_coord[i][ii] - nucl_coord[ii][a]) * rij_inv - elnuc_dist_deriv_e[3, i, a] = 2.0 * rij_inv +daccu = np.zeros(4, dtype=float) +daccu2 = np.zeros(4, dtype=float) +een_rescaled_e_deriv_e_t = een_rescaled_e_deriv_e.T +print(een_rescaled_e_deriv_e_t.shape) +for n in range(0, dim_cord_vect): + l = lkpm_of_cindex[0,n] + k = lkpm_of_cindex[1,n] + p = lkpm_of_cindex[2,n] + m = lkpm_of_cindex[3,n] -en_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,nucl_num),dtype=float) -for a in range(nucl_num): - for i in range(elec_num): - f = 1.0 - kappa * en_distance_rescaled[i][a] - for ii in range(4): - en_distance_rescaled_deriv_e[ii][i][a] = elnuc_dist_deriv_e[ii][i][a] - en_distance_rescaled_deriv_e[3][i][a] = en_distance_rescaled_deriv_e[3][i][a] + \ - (-kappa * en_distance_rescaled_deriv_e[0][i][a] * en_distance_rescaled_deriv_e[0][i][a]) + \ - (-kappa * en_distance_rescaled_deriv_e[1][i][a] * en_distance_rescaled_deriv_e[1][i][a]) + \ - (-kappa * en_distance_rescaled_deriv_e[2][i][a] * en_distance_rescaled_deriv_e[2][i][a]) - for ii in range(4): - en_distance_rescaled_deriv_e[ii][i][a] = en_distance_rescaled_deriv_e[ii][i][a] * f - -third = 1.0 / 3.0 -factor_en_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) -dx = np.zeros(shape=(4),dtype=float) -pow_ser_g = np.zeros(shape=(3),dtype=float) -for a in range(nucl_num): - for i in range(elec_num): - x = en_distance_rescaled[i][a] - if abs(x) < 1e-18: - continue - pow_ser_g = np.zeros(shape=(3),dtype=float) - den = 1.0 + aord_vector[1][type_nucl_vector[a]-1] * x - invden = 1.0 / den - invden2 = invden * invden - invden3 = invden2 * invden - xinv = 1.0 / (x + 1.0E-18) - - for ii in range(4): - dx[ii] = en_distance_rescaled_deriv_e[ii][i][a] - - lap1 = 0.0 - lap2 = 0.0 - lap3 = 0.0 - for ii in range(3): - x = en_distance_rescaled[i][a] - if x < 1e-18: - continue - for p in range(2,aord_num+1): - y = p * aord_vector[(p-1) + 1][type_nucl_vector[a]-1] * x - pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] - lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii] - lap2 = lap2 + y - x = x * en_distance_rescaled[i][a] - - lap3 = lap3 - 2.0 * aord_vector[1][type_nucl_vector[a]-1] * dx[ii] * dx[ii] - - factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + aord_vector[0][type_nucl_vector[a]-1] * \ - dx[ii] * invden2 + pow_ser_g[ii] - - ii = 3 - lap2 = lap2 * dx[ii] * third - lap3 = lap3 + den * dx[ii] - lap3 = lap3 * (aord_vector[0][type_nucl_vector[a]-1] * invden3) - factor_en_deriv_e[ii][i] = factor_en_deriv_e[ii][i] + lap1 + lap2 + lap3 - -print("factor_en_deriv_e[0][0]:",factor_en_deriv_e[0][0]) -print("factor_en_deriv_e[1][0]:",factor_en_deriv_e[1][0]) -print("factor_en_deriv_e[2][0]:",factor_en_deriv_e[2][0]) -print("factor_en_deriv_e[3][0]:",factor_en_deriv_e[3][0]) + for a in range(0, nucl_num): + cn = cord_vector_full[a][n] + for j in range(0, elec_num): + accu = 0.0 + accu2 = 0.0 + daccu = 0.0 + daccu2 = 0.0 + for i in range(0, elec_num): + accu = accu + een_rescaled_e[i,j,k] * \ + een_rescaled_n[a,i,m] + accu2 = accu2 + een_rescaled_e[i,j,k] * \ + een_rescaled_n[a,i,m+l] +# daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \ +# een_rescaled_n[a,i,m] +# daccu[0:4] = daccu[0:4] + een_rescaled_e_deriv_e_t[k,j,0:4,i,k] * \ +# een_rescaled_n[a,i,m] + accu2 = accu2 + accu * een_rescaled_n[a,j,m+l] +# factor_een = factor_een + accu2 * cn +print("factor_een:",factor_een) #+end_src #+RESULTS: - : factor_en_deriv_e[0][0]: 0.11609919541763383 - : factor_en_deriv_e[1][0]: -0.23301394780804574 - : factor_en_deriv_e[2][0]: 0.17548337641865783 - : factor_en_deriv_e[3][0]: -0.9667363412285741 + : (6, 10, 4, 10) + : factor_een: 0.0 #+begin_src c :tangle (eval c_test) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_provided(context)); -//double factor_een_deriv_e[walk_num][elec_num]; -//rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0])); +double factor_een_deriv_e[walk_num][elec_num]; +rc = qmckl_get_jastrow_factor_een_deriv_e(context, &(factor_een_deriv_e[0][0])); +assert(fabs(factor_een_deriv_e[0][0] + 0.0005481671107226865) < 1e-12); #+end_src * End of files :noexport: