TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity

Fixed sign error in jastrow gradient

Browse Source
This commit is contained in:
Anthony Scemama 2023-05-24 12:41:45 +02:00
parent 04d599649b
commit 5c019b06e3
1 changed files with 31 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
org/qmckl_jastrow_champ.org
View File

@ -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);