From ba0c93d0154978e98ca85c47c4249b9cf7d18d04 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 17 Jan 2022 15:54:21 +0100 Subject: [PATCH] Bug in finalization of MO basis --- org/qmckl_electron.org | 149 ++++++++++++++++++++++++++++++++++++++--- org/qmckl_mo.org | 10 +-- 2 files changed, 140 insertions(+), 19 deletions(-) diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 0f814ef..1f18dcc 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -394,7 +394,7 @@ qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const *** Electron coordinates Returns the current electron coordinates. The pointer is assumed - to point on a memory block of size ~3 * elec_num * walk_num~. + to point on a memory block of size ~size_max~ \ge ~3 * elec_num * walk_num~. The order of the indices is: | | Normal | Transposed | @@ -1219,14 +1219,14 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) :END: #+NAME: qmckl_ee_distance_rescaled_args - | Variable | Type | In/Out | Description | - |----------------------------------------+---------------------------+--------+--------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | + | Variable | Type | In/Out | Description | + |---------------------------+----------------------------------------+--------+--------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~coord~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, & @@ -2077,7 +2077,7 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); ** Electron-nucleus rescaled distances ~en_distance_rescaled~ stores the matrix of the rescaled distances between - electrons and nucleii. + electrons and nuclei. \[ C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa @@ -2746,7 +2746,7 @@ qmckl_exit_code qmckl_provide_en_potential(qmckl_context context) |---------------+----------------------------------------+--------+--------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nucleii | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~charge~ | ~double[nucl_num]~ | in | charge of nucleus | | ~en_distance~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-electron rescaled distances | @@ -2847,6 +2847,133 @@ rc = qmckl_get_electron_en_potential(context, &(en_pot[0])); assert (rc == QMCKL_SUCCESS); #+end_src +** Generate initial coordinates + +*** Compute :noexport: + + # begin_src f90 :comments org :tangle (eval f) :noweb yes +subroutine draw_init_points + implicit none + BEGIN_DOC +! Place randomly electrons around nuclei + END_DOC + integer :: iwalk + logical, allocatable :: do_elec(:) + integer :: acc_num + + real, allocatable :: xmin(:,:) + + integer :: i, j, k, l, kk + + real :: norm + allocate (do_elec(elec_num), xmin(3,elec_num)) + xmin = -huge(1.) + norm = 0. + do i=1,elec_alpha_num + do j=1,ao_num + norm += mo_coef_transp(i,j)*mo_coef_transp(i,j) + enddo + enddo + norm = sqrt(norm/float(elec_alpha_num)) + call rinfo( irp_here, 'Norm : ', norm ) + call rinfo( irp_here, 'mo_scale: ' , mo_scale ) + mo_coef_transp = mo_coef_transp/norm + + double precision :: qmc_ranf + real :: mo_max + do i=1,elec_alpha_num + l=1 + xmin(1,i) = mo_coef_transp(i,1)*mo_coef_transp(i,1) - 0.001*qmc_ranf() + do j=2,ao_num + xmin(2,i) = mo_coef_transp(i,j)*mo_coef_transp(i,j) - 0.001*qmc_ranf() + if (xmin(2,i) > xmin(1,i) ) then + xmin(1,i) = xmin(2,i) + l = ao_nucl(j) + endif + enddo + xmin(1,i) = nucl_coord(l,1) + xmin(2,i) = nucl_coord(l,2) + xmin(3,i) = nucl_coord(l,3) + enddo + + call iinfo(irp_here, 'Det num = ', det_num ) + do k=1,elec_beta_num + i = k+elec_alpha_num + l=1 + xmin(1,i) = mo_coef_transp(k,1)*mo_coef_transp(k,1) - 0.001*qmc_ranf() + do j=2,ao_num + xmin(2,i) = mo_coef_transp(k,j)*mo_coef_transp(k,j) - 0.001*qmc_ranf() + if (xmin(2,i) > xmin(1,i) ) then + xmin(1,i) = xmin(2,i) + l = ao_nucl(j) + endif + enddo + xmin(1,i) = nucl_coord(l,1) + xmin(2,i) = nucl_coord(l,2) + xmin(3,i) = nucl_coord(l,3) + enddo + + call rinfo( irp_here, 'time step =', time_step ) + do iwalk=1,walk_num + print *, 'Generating initial positions for walker', iwalk + acc_num = 0 + do_elec = .True. + integer :: iter + do iter = 1,10000 + if (acc_num >= elec_num) then + exit + endif + double precision :: gauss + real :: re_compute + re_compute = 0. + do while (re_compute < 1.e-6) + do i=1,elec_num + if (do_elec(i)) then + do l=1,3 + elec_coord(i,l) = xmin(l,i) + 1.5*(0.5-qmc_ranf()) + enddo + endif + enddo + TOUCH elec_coord + re_compute = minval(nucl_elec_dist(1:nucl_num,1:elec_num)) + enddo + + do i=1,elec_alpha_num + if (do_elec(i)) then + if ( mo_value_transp(i,i)**2 >= qmc_ranf()) then + acc_num += 1 + do_elec(i) = .False. + endif + endif + enddo + + do i=1,elec_beta_num + if (do_elec(i+elec_alpha_num)) then + if ( mo_value_transp(i,i+elec_alpha_num)**2 >= qmc_ranf()) then + acc_num += 1 + do_elec(i+elec_alpha_num) = .False. + endif + endif + enddo + + enddo + + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,iwalk) = elec_coord(i,l) + enddo + enddo + enddo + if (.not.is_worker) then + call ezfio_set_electrons_elec_coord_pool_size(walk_num) + call ezfio_set_electrons_elec_coord_pool(elec_coord_full) + endif + SOFT_TOUCH elec_coord elec_coord_full + deallocate (do_elec, xmin) + +end + # end_src + * End of files :noexport: #+begin_src c :tangle (eval h_private_type) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 219d0b2..1da14b1 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -344,6 +344,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context); #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, @@ -354,14 +355,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - qmckl_exit_code rc; - rc = qmckl_provide_mo_vgl(context); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_finalize_mo_basis", - NULL); - } + qmckl_exit_code rc = QMCKL_SUCCESS; return rc; } #+end_src