From 5c019b06e30dba05622c5431820067b54346eb7d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 24 May 2023 12:41:45 +0200 Subject: [PATCH] Fixed sign error in jastrow gradient --- org/qmckl_jastrow_champ.org | 62 ++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org index 65e87dc..e0e3068 100644 --- a/org/qmckl_jastrow_champ.org +++ b/org/qmckl_jastrow_champ.org @@ -2771,20 +2771,20 @@ integer function qmckl_compute_jastrow_champ_factor_ee_gl_doc_f( & do nw =1, walk_num factor_ee_gl(:,:,nw) = 0.0d0 - do i = 1, elec_num - do j = 1, elec_num - if (j==i) cycle + do j = 1, elec_num + do i = 1, elec_num + if (i == j) cycle - x = ee_distance_rescaled(j,i,nw) + x = ee_distance_rescaled(i,j,nw) denom = 1.0d0 + b_vector(2) * x invdenom = 1.0d0 / denom invdenom2 = invdenom * invdenom - dx(1) = ee_distance_rescaled_gl(1, j, i, nw) - dx(2) = ee_distance_rescaled_gl(2, j, i, nw) - dx(3) = ee_distance_rescaled_gl(3, j, i, nw) - dx(4) = ee_distance_rescaled_gl(4, j, i, nw) + dx(1) = ee_distance_rescaled_gl(1, i, j, nw) + dx(2) = ee_distance_rescaled_gl(2, i, j, nw) + dx(3) = ee_distance_rescaled_gl(3, i, j, nw) + dx(4) = ee_distance_rescaled_gl(4, i, j, nw) grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3) @@ -2855,20 +2855,20 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_ee_gl_hpc( memset(factor_ee_gl, 0, elec_num*4*walk_num*sizeof(double)); for (int nw = 0; nw < walk_num; ++nw) { - for (int i = 0; i < elec_num; ++i) { - const double* dxj = &ee_distance_rescaled_gl[4*elec_num*(i+nw*elec_num)]; - const double* xj = &ee_distance_rescaled [ elec_num*(i+nw*elec_num)]; - - for (int j = 0; j < elec_num; ++j) { + for (int j = 0; j < elec_num; ++j) { + const double* dxj = &ee_distance_rescaled_gl[4*elec_num*(j+nw*elec_num)]; + const double* xj = &ee_distance_rescaled [ elec_num*(j+nw*elec_num)]; + + for (int i = 0; i < elec_num; ++i) { if (j == i) continue; - double x = xj[j]; + double x = xj[i]; const double denom = 1.0 + b_vector[1]*x; const double invdenom = 1.0 / denom; const double invdenom2 = invdenom * invdenom; - const double* restrict dx = dxj + 4*j; + const double* restrict dx = dxj + 4*i; const double grad_c2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]; @@ -3045,16 +3045,16 @@ def make_dist_deriv(elec_coord): for i in range(elec_num): for j in range(elec_num): rij = np.linalg.norm(elec_coord[i] - elec_coord[j]) - rijm = np.linalg.norm(elec_coord[i] - elec_coord[j]+np.array((dx,0.,0.))) - rijp = np.linalg.norm(elec_coord[i] - elec_coord[j]-np.array((dx,0.,0.))) + rijm = np.linalg.norm(elec_coord[i]+np.array((dx,0.,0.)) - elec_coord[j]) + rijp = np.linalg.norm(elec_coord[i]-np.array((dx,0.,0.)) - elec_coord[j]) elec_dist_d[0, i, j] = (rijp-rijm)/(2.*dx) elec_dist_d[3, i, j] = (rijp+rijm-2.*rij)/(dx**2) - rijm = np.linalg.norm(elec_coord[i] - elec_coord[j]+np.array((0.,dx,0.))) - rijp = np.linalg.norm(elec_coord[i] - elec_coord[j]-np.array((0.,dx,0.))) + rijm = np.linalg.norm(elec_coord[i]+np.array((0.,dx,0.)) - elec_coord[j]) + rijp = np.linalg.norm(elec_coord[i]-np.array((0.,dx,0.)) - elec_coord[j]) elec_dist_d[1, i, j] = (rijp-rijm)/(2.*dx) elec_dist_d[3, i, j] += (rijp+rijm-2.*rij)/(dx**2) - rijm = np.linalg.norm(elec_coord[i] - elec_coord[j]+np.array((0.,0.,dx))) - rijp = np.linalg.norm(elec_coord[i] - elec_coord[j]-np.array((0.,0.,dx))) + rijm = np.linalg.norm(elec_coord[i]+np.array((0.,0.,dx)) - elec_coord[j]) + rijp = np.linalg.norm(elec_coord[i]-np.array((0.,0.,dx)) - elec_coord[j]) elec_dist_d[2, i, j] = (rijp-rijm)/(2.*dx) elec_dist_d[3, i, j] += (rijp+rijm-2.*rij)/(dx**2) return elec_dist_d @@ -3068,7 +3068,7 @@ def func(elec_coord): for i in range(elec_num): rij_inv = 1.0 / elec_dist[i, j] for ii in range(3): - elec_dist_gl[ii, i, j] = (elec_coord[j][ii] - elec_coord[i][ii]) * rij_inv + elec_dist_gl[ii, i, j] = (elec_coord[i][ii] - elec_coord[j][ii]) * rij_inv elec_dist_gl[3, i, j] = 2.0 * rij_inv elec_dist_gl[:, j, j] = 6.0 @@ -3134,9 +3134,9 @@ print("factor_ee_gl[3][0]:",factor_ee_gl[3][0]) : asym_one : 0.43340325572525706 : asymp_jasb[0] : 0.31567342786262853 : asymp_jasb[1] : 0.5323750557252571 - : factor_ee_gl[0][0]: 0.16364894652107934 - : factor_ee_gl[1][0]: -0.6927548119830084 - : factor_ee_gl[2][0]: 0.073267755223968 + : factor_ee_gl[0][0]: -0.16364894652107934 + : factor_ee_gl[1][0]: 0.6927548119830084 + : factor_ee_gl[2][0]: -0.073267755223968 : factor_ee_gl[3][0]: 1.5111672803213185 @@ -3149,14 +3149,14 @@ double factor_ee_gl[walk_num][4][elec_num]; rc = qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]),walk_num*4*elec_num); // check factor_ee_gl -printf("%f %f\n", factor_ee_gl[0][0][0], 0.16364894652107934); -assert(fabs(factor_ee_gl[0][0][0]-0.16364894652107934) < 1.e-12); +printf("%f %f\n", factor_ee_gl[0][0][0], -0.16364894652107934); +assert(fabs(factor_ee_gl[0][0][0]+0.16364894652107934) < 1.e-12); -printf("%f %f\n", factor_ee_gl[0][1][0],-0.6927548119830084 ); -assert(fabs(factor_ee_gl[0][1][0]+0.6927548119830084 ) < 1.e-12); +printf("%f %f\n", factor_ee_gl[0][1][0], 0.6927548119830084 ); +assert(fabs(factor_ee_gl[0][1][0]-0.6927548119830084 ) < 1.e-12); -printf("%f %f\n", factor_ee_gl[0][2][0], 0.073267755223968 ); -assert(fabs(factor_ee_gl[0][2][0]-0.073267755223968 ) < 1.e-12); +printf("%f %f\n", factor_ee_gl[0][2][0],-0.073267755223968 ); +assert(fabs(factor_ee_gl[0][2][0]+0.073267755223968 ) < 1.e-12); printf("%f %f\n", factor_ee_gl[0][3][0], 1.5111672803213185 ); assert(fabs(factor_ee_gl[0][3][0]-1.5111672803213185 ) < 1.e-12);