Bug in finalization of MO basis

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                |
@@ -1220,7 +1220,7 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context)
 
     #+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       |
@@ -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)

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;
+  qmckl_exit_code rc = QMCKL_SUCCESS;
-  rc = qmckl_provide_mo_vgl(context);
-  if (rc != QMCKL_SUCCESS) {
-    return qmckl_failwith( context,
-                           QMCKL_NOT_PROVIDED,
-                           "qmckl_finalize_mo_basis",
-                            NULL);
-  }
   return rc;
 }
    #+end_src