Merge gitlab.com:abdammar/qmcchem into abdallah

bin/jast_opt.py

@@ -8,10 +8,13 @@ import os
 import time
 import subprocess
 from math import sqrt
+#
+#from bayes_opt import BayesianOptimization
 
 QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
-
+# this should give: /home/ammar/qmcchem
 sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# this should add: /home/ammar/qmcchem/EZFIO/Python/
 
 from ezfio import ezfio
 
@@ -22,7 +25,7 @@ block_time = 20
 def main():
     if len(sys.argv) != 2:
       print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
-      sys.exti(1)
+      sys.exit(1)
 
     filename = sys.argv[1]
     ezfio.set_file(filename)
@@ -60,6 +63,21 @@ def main():
         else:
             return None, None
 
+#####################################################################################################
+    def get_energy_deriv_nucPar():
+    	# !!!
+        buffer = subprocess.check_output(['qmcchem', 'result', filename], encoding='UTF-8')
+        # !!!
+        if buffer.strip()!= "":
+        	e_der1   = [ buffer.splitlines()[13].split()[2], buffer.splitlines()[14].split()[2], buffer.splitlines()[15].split()[2] ]
+        	err_der1 = [ buffer.splitlines()[13].split()[4], buffer.splitlines()[14].split()[4], buffer.splitlines()[15].split()[4] ]
+        	e_der2   = [ buffer.splitlines()[2].split()[2], buffer.splitlines()[3].split()[2],  buffer.splitlines()[4].split()[2] ]
+        	err_der2 = [ buffer.splitlines()[2].split()[4], buffer.splitlines()[3].split()[4],  buffer.splitlines()[4].split()[4] ]
+        	return e_der1, e_der2
+        else:
+            return None
+#####################################################################################################
+
     def get_variance():
         buffer = subprocess.check_output(['qmcchem', 'result', '-e',
                                           'e_loc_qmcvar', filename],
@@ -99,28 +117,53 @@ def main():
 #                                 '-w', '10',
 #                                 filename])
         subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
-                                 '-l', str(block_time),
+                                 '-l', str(block_time), filename])
-                                 filename])
+        # !!!
-
+    # !!!
-
+    #####################################################################################################
+    # Only for CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr.
+    # !!!
+    def fprime(x):
+    	# !!!
+    	print ("derivative on: x = %s"%str(x))
+    	# !!!
+    	e_der1, e_der2 = get_energy_deriv_nucPar()
+    	e, _ = get_energy()
+    	energy_deriv_0 = 2. * ( float(e_der1[0]) - e * float(e_der2[0]) )
+    	energy_deriv_1 = 2. * ( float(e_der1[1]) - e * float(e_der2[1]) )
+    	energy_deriv = np.array([  energy_deriv_0 , energy_deriv_1 ])
+    	# !!!
+    	print(energy_deriv)
+    	# !!!
+    	return energy_deriv
+    # !!!
+    #####################################################################################################
+    # !!!
     memo_energy = {'fmin': 100000000.}
     def f(x):
     	print ("x = %s"%str(x))
+    	# !!!
     	h = str(x)
     	if h in memo_energy:
     		return memo_energy[h]
+    	# !!!
     	set_params_pen(x)
     	set_vmc_params()
+    	# !!!
     	pid = os.fork()
     	if pid == 0:
+    		# In child process
     		run_qmc()
+    		# Exit with status os.EX_OK using os._exit() method. The value of os.EX_OK is 0
     		os._exit(os.EX_OK)
     	else:
+    		# In parent process
     		import atexit
     		atexit.register(stop_qmc)
-
+    		# !!!
-            err = thresh+1.
     		time.sleep(3.*block_time/4.)
+    		# !!!
+    		err = thresh+1.
     		local_thresh = thresh
     		while err > local_thresh:
     			time.sleep(block_time)
@@ -128,14 +171,20 @@ def main():
     			variance, v_err = get_variance()
     			if e is None or variance is None:
     				continue
-                energy = e + variance
+    			# !!!
-                err = sqrt(e_err*e_err+v_err*v_err)
+    			#energy = e + variance
-                print(" %f  %f    %f  %f  %f  %f"%(e, e_err, variance, v_err, energy, err))
+    			#err = sqrt(e_err*e_err+v_err*v_err)
+    			energy = e
+    			err = e_err
+    			# !!!
+    			#print(" %f  %f    %f  %f  %f  %f"%(e, e_err, variance, v_err, energy, err))
+    			# !!!
+    			print(" %f  %f"%(energy, err))
     			if (energy-2.*err) > memo_energy['fmin']+thresh:
     				local_thresh = 10.*thresh
     			elif (energy+2.*err) < memo_energy['fmin']-thresh:
     				local_thresh = 10.*thresh
-
+    			# !!!
     			# Check if PID is still running
     			try:
     				os.kill(pid,0)
@@ -144,22 +193,34 @@ def main():
     				break
     		stop_qmc()
     	os.wait()
-        memo_energy[h] = energy + err
+    	# !!!
+    	#memo_energy[h] = energy + err
+    	memo_energy[h] = energy
+    	# !!!
     	memo_energy['fmin'] = min(energy, memo_energy['fmin'])
     	return energy
-
+    # !!!
-
     def run():
     	x = get_params_pen()
+    	# !!!
     	if sum(x) == 0.:
     		jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
     		x += jast_a_up_dn
-        opt = sp.optimize.minimize(f,x,method="Powell",
+    	# !!!
-                             options= {'disp':True, 'ftol':thresh,'xtol':0.02})
+    	opt = sp.optimize.minimize(f,x,method="Powell", options= {'disp':True, 'ftol':thresh,'xtol':0.02})
+    	# !!!
+    	bnds = ((0.001, 100), (0.001, 100)) 
+    	#opt = sp.optimize.minimize(f, x, method='CG', jac=fprime, options={'gtol': thresh, 'disp': True})
+    	#opt = sp.optimize.minimize(f, x, method='TNC', jac=fprime, options={'gtol': thresh, 'disp': True}, bounds=bnds)
+    	#opt = sp.optimize.minimize(f, x, method='Newton-CG', jac=fprime) 
+    	#opt = sp.optimize.minimize(f, x, method='trust-constr', jac=fprime)
     	print("x = "+str(opt))
     	set_params_pen(opt['x'])
-
+    	# !!!
+    # !!!
     run()
+    # !!!
+# !!!
 
 if __name__ == '__main__':
     main()

ezfio_config/qmc.config

@@ -43,7 +43,20 @@ spindeterminants
   psi_coef_matrix_rows    integer (spindeterminants_n_det)
   psi_coef_matrix_columns integer (spindeterminants_n_det)
   psi_coef_matrix_values  double precision (spindeterminants_n_det,spindeterminants_n_states)
-
+  n_svd_coefs_unique         integer
+  n_svd_coefs                integer
+  n_svd_selected             integer
+  n_svd_toselect             integer
+  psi_svd_alpha_unique       double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs_unique,spindeterminants_n_states)
+  psi_svd_beta_unique        double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs_unique,spindeterminants_n_states)
+  psi_svd_coefs_unique       double precision (spindeterminants_n_svd_coefs_unique,spindeterminants_n_states)
+  psi_svd_alpha              double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs,spindeterminants_n_states)
+  psi_svd_beta               double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs,spindeterminants_n_states)
+  psi_svd_coefs              double precision (spindeterminants_n_svd_coefs,spindeterminants_n_states)
+  psi_svd_alpha_numselected  integer          (spindeterminants_n_svd_selected,spindeterminants_n_states)
+  psi_svd_beta_numselected   integer          (spindeterminants_n_svd_selected,spindeterminants_n_states)
+  psi_svd_alpha_numtoselect  integer          (spindeterminants_n_svd_toselect,spindeterminants_n_states)
+  psi_svd_beta_numtoselect   integer          (spindeterminants_n_svd_toselect,spindeterminants_n_states)
 
 simulation
   do_run              integer

src/JASTROW/jastrow_simple.irp.f

@@ -4,7 +4,7 @@
 BEGIN_PROVIDER [ double precision , jast_elec_Simple_value, (elec_num_8)  ]
 implicit none
  BEGIN_DOC  
-! J(i) = \sum_j a.rij/(1+b^2.rij) - \sum_A (a.riA/(1+a.riA))^2
+! J(i) = \sum_j a.rij/(1+b.rij) - \sum_A (a.riA/(1+a.riA))^2
  END_DOC
  integer :: i,j
  double precision :: a, b, rij, tmp
@@ -20,10 +20,11 @@ implicit none
   enddo
  enddo
 
- a = 0.5*jast_a_up_up
+ a = 0.5d0*jast_a_up_up
  b = jast_b_up_up
 
- do j=1,elec_alpha_num
+ ! parallele spin up-up
+ do j=1,elec_alpha_num-1
   !DIR$ LOOP COUNT (50)
   do i=j+1,elec_alpha_num
     rij = elec_dist(i,j)
@@ -33,7 +34,8 @@ implicit none
   enddo
  enddo
 
- do j=elec_alpha_num+1,elec_num
+ ! parallele spin dn-dn
+ do j=elec_alpha_num+1,elec_num-1
   !DIR$ LOOP COUNT (50)
   do i=j+1,elec_num
     rij = elec_dist(i,j)
@@ -43,9 +45,11 @@ implicit none
   enddo
  enddo
 
- a = 0.5*jast_a_up_dn
+
+ a = 0.5d0*jast_a_up_dn
  b = jast_b_up_dn
 
+ ! anti-parallele spin
  do j=1,elec_alpha_num
   !DIR$ LOOP COUNT (50)
   do i=elec_alpha_num+1,elec_num
@@ -190,5 +194,31 @@ BEGIN_PROVIDER [ double precision , jast_elec_Simple_lapl, (elec_num_8) ]
   enddo
  enddo
 
+END_PROVIDER
+
+
+
+
+BEGIN_PROVIDER[ double precision, jast_elec_Simple_deriv_nucPar, (nucl_num) ]
+  implicit none
+  BEGIN_DOC  
+  ! Variation of the Jastrow factor with respect to nuclear parameters
+  END_DOC
+
+  integer          :: i, j
+  double precision :: a, rij, tmp1, tmp2
+
+  do j = 1, nucl_num
+    a    = jast_pen(j)
+    tmp2 = 0.d0
+    !DIR$ LOOP COUNT (100)
+    do i = 1, elec_num
+      rij  = nucl_elec_dist(j,i)
+      tmp1 = (1.d0+a*rij)*(1.d0+a*rij)*(1.d0+a*rij)
+      tmp2 += rij*rij/tmp1
+    end do
+    jast_elec_Simple_deriv_nucPar(j) = -2.d0 * a * tmp2
+  end do
 
 END_PROVIDER
+

src/MAIN/qmcchem_info.irp.f

@@ -13,6 +13,9 @@ program qmcchem_info
   endif
   print *,  'Number of determinants                   : ', det_num
   print *,  'Number of unique alpha/beta determinants : ', det_alpha_num, det_beta_num
+  if (utilise_svd) then
+    print *,  'SVD rank                                 : ', n_svd_coefs
+  endif
   print *,  'Closed-shell MOs                         : ', mo_closed_num
   print *,  'Number of MOs in determinants            : ', num_present_mos
 !  print *,  'Det alpha norm:'

src/MAIN/test_debug.irp.f

@@ -0,0 +1,51 @@
+program qmcchem_info
+  implicit none
+  PROVIDE ezfio_filename
+  double precision               :: cpu0, cpu1
+  character*(8)                  :: str_n
+  integer                        :: iargc
+  integer                        :: imax
+  if (command_argument_count() > 1) then
+    call get_command_argument(2,str_n)
+    read(str_n,*) imax
+  else
+    imax = 100
+  endif
+  print *,  'Number of determinants                   : ', det_num
+  print *,  'Number of unique alpha/beta determinants : ', det_alpha_num, det_beta_num
+  print *,  'Closed-shell MOs                         : ', mo_closed_num
+  print *,  'Number of MOs in determinants            : ', num_present_mos
+!  print *,  'Det alpha norm:'
+!  print *,  det_alpha_norm
+!  print *,  'Det beta norm:'
+!  print *,  det_beta_norm
+  call step1
+  call cpu_time (cpu0)
+  call step2(imax)
+  call cpu_time (cpu1)
+  print *,  'Time for the calculation of E_loc (ms)   : ', 1000.*(cpu1-cpu0)/float(imax)
+end
+
+subroutine step1
+  implicit none
+  print *,  'E_loc                                    : ', E_loc
+  print *,  'E_loc_svd                                : ', E_loc_svd
+  integer :: i
+  do i=1,elec_num
+    print *, 'x', psi_grad_psi_inv_x(i), psi_grad_psi_inv_x_SVD(i)
+    print *, 'y', psi_grad_psi_inv_y(i), psi_grad_psi_inv_y_SVD(i)
+    print *, 'z', psi_grad_psi_inv_z(i), psi_grad_psi_inv_z_SVD(i)
+    print *, 'l', psi_lapl_psi_inv(i), psi_lapl_psi_inv_SVD(i)
+    print *, ''
+  enddo
+end
+
+subroutine step2(imax)
+  implicit none
+  integer, intent(in)            :: imax
+  integer                        :: i
+  do i=1,imax
+    PROVIDE E_loc
+    TOUCH elec_coord
+  enddo
+end

src/PROPERTIES/properties_buildpsi.irp.f

@@ -0,0 +1,312 @@
+
+BEGIN_PROVIDER [ double precision, hij_fm, (size_hij_fm) ]
+
+  implicit none
+
+  BEGIN_DOC
+  ! !!!
+  ! hij  =  < psi_svd J | H | J l l' > / < psi_svd J | psi_svd J > 
+  !      =  < H (J l l')/(psi_svd J) >    ( first method: fm )
+  !      =  < E_loc (l l') / psi_svd >    ( secnd method: sm )
+  ! Dimensions : n_svd_toselect
+  END_DOC
+
+  integer          :: i, l, lp, e
+  double precision :: f, g, h, T, V
+
+  do i = 1, n_svd_toselect
+
+    l  = psi_svd_alpha_numtoselect(i,1)
+    lp = psi_svd_beta_numtoselect (i,1)
+
+    ! Lapl D 
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(4,e,l) * det_beta_value_SVD_unique(lp)
+    enddo
+    do e = elec_alpha_num+1, elec_num
+      g += det_alpha_value_SVD_unique(l) * det_beta_grad_lapl_SVD_unique(4,e,lp)
+    enddo
+    T = g
+
+    ! D (Lapl J)/J
+    g = 0.d0
+    do e = 1, elec_num
+      g += jast_lapl_jast_inv(e)
+    enddo
+    T += det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp) * g
+
+    ! 2 (grad D).(Grad J)/J
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(1,e,l) * jast_grad_jast_inv_x(e) + &
+           det_alpha_grad_lapl_SVD_unique(2,e,l) * jast_grad_jast_inv_y(e) + &
+           det_alpha_grad_lapl_SVD_unique(3,e,l) * jast_grad_jast_inv_z(e)
+    enddo
+    h = 0.d0
+    do e = elec_alpha_num+1, elec_num
+      h += det_beta_grad_lapl_SVD_unique(1,e,lp) * jast_grad_jast_inv_x(e) + &
+           det_beta_grad_lapl_SVD_unique(2,e,lp) * jast_grad_jast_inv_y(e) + &
+           det_beta_grad_lapl_SVD_unique(3,e,lp) * jast_grad_jast_inv_z(e)
+    enddo
+    T += 2.d0 * ( g * det_beta_value_SVD_unique(lp) + h * det_alpha_value_SVD_unique(l) )
+
+    g = det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp)
+    V = E_pot * g
+    ! TODO
+    !do e = 1, elec_alpha_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_pseudo_SVD_unique(e,l) * det_beta_value_SVD_unique(lp)
+    !enddo
+    !do e = elec_alpha_num+1, elec_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_value_SVD_unique(l) * det_beta_pseudo_SVD_unique(e,lp)
+    !enddo
+    f = -0.5d0*T + V
+    f *= psidet_inv_SVD
+
+    hij_fm(i) = f
+
+  enddo
+
+  hij_fm_min = min( hij_fm_min, minval(hij_fm) )
+  hij_fm_max = max( hij_fm_max, maxval(hij_fm) )
+
+  SOFT_TOUCH hij_fm_min hij_fm_max
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, hij_sm, (size_hij_sm) ]
+
+  implicit none
+
+  BEGIN_DOC
+  ! !!!
+  ! hij  =  < psi_svd J | H | J l l' > / < psi_svd J | psi_svd J > 
+  !      =  < H (J l l')/(psi_svd J) >    ( first method: fm)
+  !      =  < E_loc (l l') / psi_svd >    ( secnd method: sm)
+  ! Dimensions : n_svd_toselect
+  END_DOC
+
+  integer :: i, l, lp
+
+  do i = 1, n_svd_toselect
+
+    l  = psi_svd_alpha_numtoselect(i,1)
+    lp = psi_svd_beta_numtoselect (i,1)
+
+    hij_sm(i) = E_loc * det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp) * psidet_inv_SVD
+
+  enddo
+
+  hij_sm_min = min( hij_sm_min, minval(hij_sm) )
+  hij_sm_max = max( hij_sm_max, maxval(hij_sm) )
+
+  SOFT_TOUCH hij_sm_min hij_sm_max
+
+END_PROVIDER
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, xij_diag, (size_xij_diag) ]
+
+  implicit none
+
+  BEGIN_DOC
+  ! !!!
+  ! < l l' | H | l l' >
+  ! Dimensions : n_svd_toselect
+  END_DOC
+
+  integer          :: i, l, lp, e
+  double precision :: f, g, h, T, V
+
+  do i = 1, n_svd_toselect
+  
+    l  = psi_svd_alpha_numtoselect(i,1)
+    lp = psi_svd_beta_numtoselect (i,1)
+  
+    ! Lapl D 
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(4,e,l) * det_beta_value_SVD_unique(lp)
+    enddo
+    do e = elec_alpha_num+1, elec_num
+      g += det_alpha_value_SVD_unique(l) * det_beta_grad_lapl_SVD_unique(4,e,lp)
+    enddo
+    T = g
+
+    ! D (Lapl J)/J
+    g = 0.d0
+    do e = 1, elec_num
+      g += jast_lapl_jast_inv(e)
+    enddo
+    T += det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp) * g
+
+    ! 2 (grad D).(Grad J)/J
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(1,e,l) * jast_grad_jast_inv_x(e) + &
+           det_alpha_grad_lapl_SVD_unique(2,e,l) * jast_grad_jast_inv_y(e) + &
+           det_alpha_grad_lapl_SVD_unique(3,e,l) * jast_grad_jast_inv_z(e)
+    enddo
+    h = 0.d0
+    do e = elec_alpha_num+1, elec_num
+      h += det_beta_grad_lapl_SVD_unique(1,e,lp) * jast_grad_jast_inv_x(e) + &
+           det_beta_grad_lapl_SVD_unique(2,e,lp) * jast_grad_jast_inv_y(e) + &
+           det_beta_grad_lapl_SVD_unique(3,e,lp) * jast_grad_jast_inv_z(e)
+    enddo
+    T += 2.d0 * ( g * det_beta_value_SVD_unique(lp) + h * det_alpha_value_SVD_unique(l) )
+
+    g = det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp)
+    V = E_pot * g
+    ! TODO
+    !do e = 1, elec_alpha_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_pseudo_SVD_unique(e,l) * det_beta_value_SVD_unique(lp)
+    !enddo
+    !do e = elec_alpha_num+1, elec_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_value_SVD_unique(l) * det_beta_pseudo_SVD_unique(e,lp)
+    !enddo
+    f = -0.5d0*T + V
+    f *= psidet_inv_SVD * psidet_inv_SVD
+
+    xij_diag(i) = f * det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp)
+
+  enddo
+
+  xij_diag_min = min( xij_diag_min, minval(xij_diag) )
+  xij_diag_max = max( xij_diag_max, maxval(xij_diag) )
+
+  SOFT_TOUCH xij_diag_min xij_diag_max
+
+END_PROVIDER
+
+
+
+
+BEGIN_PROVIDER [ double precision, overlop_selected_matrix, (size_overlop_selected_matrix) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < k_selected k'_selected | l_selected l'_selected >
+  ! Dimensions : n_svd_selected * n_svd_selected
+  END_DOC
+
+  integer          :: k, kp, l, lp
+  integer          :: i, j, ii0, ii
+  double precision :: f
+
+  do i = 1, n_svd_selected
+    ii0 = (i-1)*n_svd_selected
+
+    l  = psi_svd_alpha_numselected(i,1)
+    lp = psi_svd_beta_numselected (i,1)
+
+    f = det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp) * psidet_inv_SVD * psidet_inv_SVD
+
+    do j = 1, n_svd_selected
+      ii = ii0 + j
+
+      k  = psi_svd_alpha_numselected(j,1)
+      kp = psi_svd_beta_numselected (j,1)
+
+      overlop_selected_matrix(ii) = det_alpha_value_SVD_unique(k) * det_beta_value_SVD_unique(kp) * f
+
+    enddo
+  enddo
+
+  overlop_selected_matrix_min = min(overlop_selected_matrix_min,minval(overlop_selected_matrix))
+  overlop_selected_matrix_max = max(overlop_selected_matrix_max,maxval(overlop_selected_matrix))
+  SOFT_TOUCH overlop_selected_matrix_min overlop_selected_matrix_max
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, h_selected_matrix, (size_h_selected_matrix) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < k_selected k'_selected | H | l_selected l'_selected >
+  ! Dimensions : n_svd_selected * n_svd_selected
+  END_DOC
+
+  integer          :: k, kp, l, lp
+  integer          :: i, j, ii0, ii
+  integer          :: e
+  double precision :: f, g, h, T, V
+
+  do i = 1, n_svd_selected
+    ii0 = (i-1)*n_svd_selected
+
+    l  = psi_svd_alpha_numselected(i,1)
+    lp = psi_svd_beta_numselected (i,1)
+
+    ! Lapl D 
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(4,e,l) * det_beta_value_SVD_unique(lp)
+    enddo
+    do e = elec_alpha_num+1, elec_num
+      g += det_alpha_value_SVD_unique(l) * det_beta_grad_lapl_SVD_unique(4,e,lp)
+    enddo
+    T = g
+
+    ! D (Lapl J)/J
+    g = 0.d0
+    do e = 1, elec_num
+      g += jast_lapl_jast_inv(e)
+    enddo
+    T += det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp) * g
+
+    ! 2 (grad D).(Grad J)/J
+    g = 0.d0
+    do e = 1, elec_alpha_num
+      g += det_alpha_grad_lapl_SVD_unique(1,e,l) * jast_grad_jast_inv_x(e) + &
+           det_alpha_grad_lapl_SVD_unique(2,e,l) * jast_grad_jast_inv_y(e) + &
+           det_alpha_grad_lapl_SVD_unique(3,e,l) * jast_grad_jast_inv_z(e)
+    enddo
+    h = 0.d0
+    do e = elec_alpha_num+1, elec_num
+      h += det_beta_grad_lapl_SVD_unique(1,e,lp) * jast_grad_jast_inv_x(e) + &
+           det_beta_grad_lapl_SVD_unique(2,e,lp) * jast_grad_jast_inv_y(e) + &
+           det_beta_grad_lapl_SVD_unique(3,e,lp) * jast_grad_jast_inv_z(e)
+    enddo
+    T += 2.d0 * ( g * det_beta_value_SVD_unique(lp) + h * det_alpha_value_SVD_unique(l) )
+
+    g = det_alpha_value_SVD_unique(l) * det_beta_value_SVD_unique(lp)
+    V = E_pot * g
+    ! TODO
+    ! ajouter le terme pseudo
+    !do e = 1, elec_alpha_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_pseudo_SVD_unique(e,l) * det_beta_value_SVD_unique(lp)
+    !enddo
+    !do e = elec_alpha_num+1, elec_num
+    !  V -= pseudo_non_local_SVD(e) * g
+    !  V += det_alpha_value_SVD_unique(l) * det_beta_pseudo_SVD_unique(e,lp)
+    !enddo
+    f = -0.5d0*T + V
+    f *= psidet_inv_SVD * psidet_inv_SVD
+
+    do j = 1, n_svd_selected
+      ii = ii0 + j
+
+      k  = psi_svd_alpha_numselected(j,1)
+      kp = psi_svd_beta_numselected (j,1)
+
+      h_selected_matrix(ii) = f * det_alpha_value_SVD_unique(k) * det_beta_value_SVD_unique(kp)
+
+    enddo
+  enddo
+
+  h_selected_matrix_min = min(h_selected_matrix_min,minval(h_selected_matrix))
+  h_selected_matrix_max = max(h_selected_matrix_max,maxval(h_selected_matrix))
+  SOFT_TOUCH h_selected_matrix_min h_selected_matrix_max
+END_PROVIDER
+

src/PROPERTIES/properties_ci.irp.f

@@ -23,7 +23,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ]
   implicit none
   BEGIN_DOC
-  ! < Phi_0 | det(j) >
+  ! < Phi_0 | H | det(j) >
   !
   ! Dimensions : det_num
   END_DOC
@@ -54,7 +54,7 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, ci_overlap_matrix, (size_ci_overlap_matrix) ]
   implicit none
   BEGIN_DOC
-  ! < det(i) |H| det(j) >
+  ! < det(i) | det(j) >
   !
   ! Dimensions : det_num*det_num
   END_DOC

src/PROPERTIES/properties_ci_SVD.irp.f

@@ -0,0 +1,221 @@
+
+
+
+BEGIN_PROVIDER [ double precision, ci_overlap_psidet_SVD, (size_ci_overlap_psidet_SVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < psi_0 | det(j) >
+  ! Dimensions : n_svd_coefs
+  END_DOC
+
+  integer :: k
+  do k = 1, n_svd_coefs
+    ci_overlap_psidet_SVD(k) = det_alpha_value_SVD(k) * det_beta_value_SVD(k) * psidet_inv_SVD
+  enddo
+
+  ci_overlap_psidet_SVD_min = min(ci_overlap_psidet_SVD_min,minval(ci_overlap_psidet_SVD))
+  ci_overlap_psidet_SVD_max = max(ci_overlap_psidet_SVD_max,maxval(ci_overlap_psidet_SVD))
+  SOFT_TOUCH ci_overlap_psidet_SVD_min ci_overlap_psidet_SVD_max
+END_PROVIDER
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_psidet_SVD, (size_ci_h_psidet_SVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < psi_0 |H| det(j) >
+  ! Dimensions : n_svd_coefs
+  END_DOC
+
+  integer          :: k, e
+  double precision :: T
+
+  do k = 1, n_svd_coefs
+    T = 0.d0
+    do e = 1, elec_alpha_num
+      T += det_alpha_grad_lapl_SVD(4,e,k) * det_beta_value_SVD(k)
+    enddo
+    do e = elec_beta_num+1, elec_num
+      T += det_beta_grad_lapl_SVD(4,e,k) * det_alpha_value_SVD(k)
+    enddo
+    ci_h_psidet_SVD(k) = -0.5d0*T + E_pot * det_alpha_value_SVD(k) * det_beta_value_SVD(k)
+    ci_h_psidet_SVD(k) *= psidet_inv_SVD
+  enddo
+
+  ci_h_psidet_SVD_min = min(ci_h_psidet_SVD_min,minval(ci_h_psidet_SVD))
+  ci_h_psidet_SVD_max = max(ci_h_psidet_SVD_max,maxval(ci_h_psidet_SVD))
+  SOFT_TOUCH ci_h_psidet_SVD_min ci_h_psidet_SVD_max
+END_PROVIDER
+
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_overlap_matrix_SVD, (size_ci_overlap_matrix_SVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < det(i) | det(j) >
+  ! Dimensions : n_svd_coefs * n_svd_coefs
+  END_DOC
+
+  integer          :: k, l
+  double precision :: f
+
+  do k = 1, n_svd_coefs
+    f = det_alpha_value_SVD(k) * det_beta_value_SVD(k) * psidet_inv_SVD * psidet_inv_SVD
+    do l = 1, n_svd_coefs
+      ci_overlap_matrix_SVD( n_svd_coefs*(k-1) + l) = det_alpha_value_SVD(l) * det_beta_value_SVD(l) * f
+    enddo
+  enddo
+
+  ci_overlap_matrix_SVD_min = min(ci_overlap_matrix_SVD_min,minval(ci_overlap_matrix_SVD))
+  ci_overlap_matrix_SVD_max = max(ci_overlap_matrix_SVD_max,maxval(ci_overlap_matrix_SVD))
+  SOFT_TOUCH ci_overlap_matrix_SVD_min ci_overlap_matrix_SVD_max
+END_PROVIDER
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix_SVD, (size_ci_h_matrix_SVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < det(i) | H | det(j) >
+  ! Dimensions : n_svd_coefs * n_svd_coefs
+  END_DOC
+
+  integer          :: k, l, e
+  double precision :: f, g, h, T, V
+
+  do l = 1, n_svd_coefs
+      ! Lapl D 
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(l)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,l)
+      enddo
+      T = g
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e = 1, elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value_SVD(l) * det_beta_value_SVD(l) * g
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += &
+          det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e = elec_alpha_num+1, elec_num
+        h += &
+          det_beta_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0 * ( g * det_beta_value_SVD(l) + h * det_alpha_value_SVD(l) )
+      g = det_alpha_value_SVD(l) * det_beta_value_SVD(l)
+      V = E_pot * g
+      do e = 1, elec_alpha_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(l)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,l)
+      enddo
+      f = -0.5d0*T + V
+      f *= psidet_inv_SVD * psidet_inv_SVD
+      do k = 1, n_svd_coefs
+        ci_h_matrix_SVD( n_svd_coefs*(l-1) + k) = f * det_alpha_value_SVD(k) * det_beta_value_SVD(k)
+      enddo
+  enddo
+
+  ci_h_matrix_SVD_min = min(ci_h_matrix_SVD_min,minval(ci_h_matrix_SVD))
+  ci_h_matrix_SVD_max = max(ci_h_matrix_SVD_max,maxval(ci_h_matrix_SVD))
+  SOFT_TOUCH ci_h_matrix_SVD_min ci_h_matrix_SVD_max
+END_PROVIDER
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix_diag_SVD, (size_ci_h_matrix_diag_SVD) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) |H| det(j) >
+  !
+  ! Dimensions : n_svd_coefs
+  END_DOC
+
+  integer          :: l, e
+  double precision :: f, g, h, T, V
+
+  do l = 1, n_svd_coefs
+      ! Lapl D 
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(l)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,l)
+      enddo
+      T = g
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e = 1, elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value(l) * det_beta_value_SVD(l) * g
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e = 1 , elec_alpha_num
+        g += &
+          det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e = elec_alpha_num+1, elec_num
+        h += &
+          det_beta_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0 * ( g * det_beta_value_SVD(l) + h * det_alpha_value_SVD(l) )
+      g = det_alpha_value_SVD(l) * det_beta_value_SVD(l)
+      V = E_pot * g
+      do e = 1, elec_alpha_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(l)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,l)
+      enddo
+      f = -0.5d0*T + V
+      f *= psidet_inv_SVD * psidet_inv_SVD
+      ci_h_matrix_diag_SVD(l) = f * det_alpha_value_SVD(l) * det_beta_value_SVD(l)
+  enddo
+
+  ci_h_matrix_diag_SVD_min = min(ci_h_matrix_diag_SVD_min,minval(ci_h_matrix_diag_SVD))
+  ci_h_matrix_diag_SVD_max = max(ci_h_matrix_diag_SVD_max,maxval(ci_h_matrix_diag_SVD))
+  SOFT_TOUCH ci_h_matrix_diag_SVD_min ci_h_matrix_diag_SVD_max
+END_PROVIDER
+
+

src/PROPERTIES/properties_ci_postSVD.irp.f

@@ -0,0 +1,258 @@
+
+
+
+BEGIN_PROVIDER [ double precision, ci_overlap_psidet_postSVD, (size_ci_overlap_psidet_postSVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < psi_0 | det(j) >
+  ! Dimensions : n_svd_coefs2
+  END_DOC
+
+  integer :: k, kp
+  do k = 1, n_svd_coefs
+    do kp = 1, n_svd_coefs
+      ci_overlap_psidet_postSVD(kp+(k-1)*n_svd_coefs) = det_alpha_value_SVD(k) * det_beta_value_SVD(kp) * psidet_inv_SVD
+    enddo
+  enddo
+
+  ci_overlap_psidet_postSVD_min = min(ci_overlap_psidet_postSVD_min,minval(ci_overlap_psidet_postSVD))
+  ci_overlap_psidet_postSVD_max = max(ci_overlap_psidet_postSVD_max,maxval(ci_overlap_psidet_postSVD))
+  SOFT_TOUCH ci_overlap_psidet_postSVD_min ci_overlap_psidet_postSVD_max
+END_PROVIDER
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_psidet_postSVD, (size_ci_h_psidet_postSVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < psi_0 |H| det(j) >
+  ! Dimensions : n_svd_coefs2
+  END_DOC
+
+  integer          :: k, kp, e
+  double precision :: T
+
+  do k = 1, n_svd_coefs
+    do kp = 1, n_svd_coefs
+      T = 0.d0
+      do e = 1, elec_alpha_num
+        T += det_alpha_grad_lapl_SVD(4,e,k) * det_beta_value_SVD(kp)
+      enddo
+      do e = elec_beta_num+1, elec_num
+        T += det_alpha_value_SVD(k) * det_beta_grad_lapl_SVD(4,e,kp)
+      enddo
+      ci_h_psidet_postSVD(kp+(k-1)*n_svd_coefs) = -0.5d0*T + E_pot * det_alpha_value_SVD(k) * det_beta_value_SVD(kp)
+      ci_h_psidet_postSVD(kp+(k-1)*n_svd_coefs) *= psidet_inv_SVD
+    enddo
+  enddo
+
+  ci_h_psidet_postSVD_min = min(ci_h_psidet_postSVD_min,minval(ci_h_psidet_postSVD))
+  ci_h_psidet_postSVD_max = max(ci_h_psidet_postSVD_max,maxval(ci_h_psidet_postSVD))
+  SOFT_TOUCH ci_h_psidet_postSVD_min ci_h_psidet_postSVD_max
+END_PROVIDER
+
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_overlap_matrix_postSVD, (size_ci_overlap_matrix_postSVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < det(i) | det(j) >
+  ! Dimensions : n_svd_coefs2 * n_svd_coefs2
+  END_DOC
+
+  integer          :: k, kp, l, lp
+  integer          :: ii0, ii1, ii2, ii
+  double precision :: f
+
+  do k = 1, n_svd_coefs
+    ii0 = (k-1)*n_svd_coefs3
+    do kp = 1, n_svd_coefs
+      ii1 = ii0 + (kp-1)*n_svd_coefs2
+
+      f = det_alpha_value_SVD(k) * det_beta_value_SVD(kp) * psidet_inv_SVD * psidet_inv_SVD
+
+      do l = 1, n_svd_coefs
+        ii2 = ii1 + (l-1)*n_svd_coefs
+        do lp = 1, n_svd_coefs
+          ii = ii2 + lp
+          ci_overlap_matrix_postSVD(ii) = det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * f
+        enddo
+      enddo
+    enddo
+  enddo
+
+  ci_overlap_matrix_postSVD_min = min(ci_overlap_matrix_postSVD_min,minval(ci_overlap_matrix_postSVD))
+  ci_overlap_matrix_postSVD_max = max(ci_overlap_matrix_postSVD_max,maxval(ci_overlap_matrix_postSVD))
+  SOFT_TOUCH ci_overlap_matrix_postSVD_min ci_overlap_matrix_postSVD_max
+END_PROVIDER
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix_postSVD, (size_ci_h_matrix_postSVD) ]
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! < det(i) | H | det(j) >
+  ! Dimensions : n_svd_coefs2 * n_svd_coefs2
+  END_DOC
+
+  integer          :: k, kp, l, lp, e
+  integer          :: ii0, ii1, ii2, ii
+  double precision :: f, g, h, T, V
+
+  do l = 1, n_svd_coefs
+    ii0 = (l-1)*n_svd_coefs3
+    do lp = 1, n_svd_coefs
+      ii1 = ii0 + (lp-1)*n_svd_coefs2
+
+      ! Lapl D 
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(lp)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,lp)
+      enddo
+      T = g
+
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e = 1, elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * g
+
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += &
+          det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e = elec_alpha_num+1, elec_num
+        h += &
+          det_beta_grad_lapl_SVD(1,e,lp) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl_SVD(2,e,lp) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl_SVD(3,e,lp) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0 * ( g * det_beta_value_SVD(lp) + h * det_alpha_value_SVD(l) )
+
+      g = det_alpha_value_SVD(l) * det_beta_value_SVD(lp)
+      V = E_pot * g
+      do e = 1, elec_alpha_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(lp)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,lp)
+      enddo
+
+      f = -0.5d0*T + V
+      f *= psidet_inv_SVD * psidet_inv_SVD
+
+      do k = 1, n_svd_coefs
+        ii2 = ii1 + (k-1)*n_svd_coefs
+        do kp = 1, n_svd_coefs
+          ii = ii2 + kp
+          ci_h_matrix_postSVD(ii) = f * det_alpha_value_SVD(k) * det_beta_value_SVD(kp)
+        enddo
+      enddo
+
+    enddo
+  enddo
+
+  ci_h_matrix_postSVD_min = min(ci_h_matrix_postSVD_min,minval(ci_h_matrix_postSVD))
+  ci_h_matrix_postSVD_max = max(ci_h_matrix_postSVD_max,maxval(ci_h_matrix_postSVD))
+  SOFT_TOUCH ci_h_matrix_postSVD_min ci_h_matrix_postSVD_max
+END_PROVIDER
+
+
+
+
+
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix_diag_postSVD, (size_ci_h_matrix_diag_postSVD) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) |H| det(j) >
+  !
+  ! Dimensions : n_svd_coefs2
+  END_DOC
+
+  integer          :: l, lp, e
+  integer          :: ii0, ii
+  double precision :: f, g, h, T, V
+
+  do l = 1, n_svd_coefs
+    ii0 = (l-1)*n_svd_coefs
+    do lp = 1, n_svd_coefs
+      ii = ii0 + lp
+
+      ! Lapl D 
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += det_alpha_grad_lapl_SVD(4,e,l) * det_beta_value_SVD(lp)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        g += det_alpha_value_SVD(l) * det_beta_grad_lapl_SVD(4,e,lp)
+      enddo
+      T = g
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e = 1, elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value_SVD(l) * det_beta_value_SVD(lp) * g
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e = 1, elec_alpha_num
+        g += &
+          det_alpha_grad_lapl_SVD(1,e,l) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl_SVD(2,e,l) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl_SVD(3,e,l) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e = elec_alpha_num+1, elec_num
+        h += &
+          det_beta_grad_lapl_SVD(1,e,lp) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl_SVD(2,e,lp) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl_SVD(3,e,lp) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0 * ( g * det_beta_value_SVD(lp) + h * det_alpha_value_SVD(l) )
+      g = det_alpha_value_SVD(l) * det_beta_value_SVD(lp)
+      V = E_pot * g
+      do e = 1, elec_alpha_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_pseudo_SVD(e,l) * det_beta_value_SVD(lp)
+      enddo
+      do e = elec_alpha_num+1, elec_num
+        V -= pseudo_non_local_SVD(e) * g
+        V += det_alpha_value_SVD(l) * det_beta_pseudo_SVD(e,lp)
+      enddo
+      f = -0.5d0*T + V
+      f *= psidet_inv_SVD * psidet_inv_SVD
+      ci_h_matrix_diag_postSVD(ii) = f * det_alpha_value_SVD(l) * det_beta_value_SVD(lp)
+    enddo
+  enddo
+
+  ci_h_matrix_diag_postSVD_min = min(ci_h_matrix_diag_postSVD_min,minval(ci_h_matrix_diag_postSVD))
+  ci_h_matrix_diag_postSVD_max = max(ci_h_matrix_diag_postSVD_max,maxval(ci_h_matrix_diag_postSVD))
+  SOFT_TOUCH ci_h_matrix_diag_postSVD_min ci_h_matrix_diag_postSVD_max
+END_PROVIDER
+
+

src/PROPERTIES/properties_deriv_jast.irp.f

@@ -0,0 +1,235 @@
+
+BEGIN_PROVIDER  [ double precision, E_deriv_nucPar_loc1, (size_E_deriv_nucPar_loc1) ]
+  implicit none
+  BEGIN_DOC
+  ! Local energy variation  with respect to nuclear parameters
+  ! 
+  ! Dimensions : nucl_num
+  END_DOC
+
+  integer :: i
+
+  do i = 1, nucl_num
+      E_deriv_nucPar_loc1(i) = E_loc*jast_elec_Simple_deriv_nucPar(i)
+  enddo
+
+  E_deriv_nucPar_loc1_min = min(E_deriv_nucPar_loc1_min,minval(E_deriv_nucPar_loc1))
+  E_deriv_nucPar_loc1_max = max(E_deriv_nucPar_loc1_max,maxval(E_deriv_nucPar_loc1))
+  SOFT_TOUCH E_deriv_nucPar_loc1_min E_deriv_nucPar_loc1_max
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER  [ double precision, E_deriv_nucPar_loc2, (size_E_deriv_nucPar_loc2) ]
+  implicit none
+  BEGIN_DOC
+  ! Local energy variation  with respect to nuclear parameters
+  ! 
+  ! Dimensions : nucl_num
+  END_DOC
+
+  integer :: i
+
+  do i=1,nucl_num
+      E_deriv_nucPar_loc2(i) = jast_elec_Simple_deriv_nucPar(i)
+  enddo
+
+  E_deriv_nucPar_loc2_min = min(E_deriv_nucPar_loc2_min,minval(E_deriv_nucPar_loc2))
+  E_deriv_nucPar_loc2_max = max(E_deriv_nucPar_loc2_max,maxval(E_deriv_nucPar_loc2))
+  SOFT_TOUCH E_deriv_nucPar_loc2_min E_deriv_nucPar_loc2_max
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER  [ double precision, E_deriv_bPar_loc1 ]
+  implicit none
+  BEGIN_DOC
+  ! Local energy variation  with respect to parameter b
+  END_DOC
+
+  E_deriv_bPar_loc1 = E_loc*J_deriv_bPar_ex
+
+  E_deriv_bPar_loc1_min = min(E_deriv_bPar_loc1_min, E_deriv_bPar_loc1)
+  E_deriv_bPar_loc1_max = max(E_deriv_bPar_loc1_max, E_deriv_bPar_loc1)
+  SOFT_TOUCH E_deriv_bPar_loc1_min E_deriv_bPar_loc1_max
+END_PROVIDER
+
+BEGIN_PROVIDER  [ double precision, E_deriv_bPar_loc2 ]
+  implicit none
+  BEGIN_DOC
+  ! Local energy variation  with respect to parameter b
+  END_DOC
+
+  E_deriv_bPar_loc2 = J_deriv_bPar_ex
+
+  E_deriv_bPar_loc2_min = min(E_deriv_bPar_loc2_min,  E_deriv_bPar_loc2)
+  E_deriv_bPar_loc2_max = max(E_deriv_bPar_loc2_max,  E_deriv_bPar_loc2)
+  SOFT_TOUCH E_deriv_bPar_loc2_min E_deriv_bPar_loc2_max
+END_PROVIDER
+
+
+
+
+
+BEGIN_PROVIDER  [ double precision, J_deriv_nucPar_verif, (size_J_deriv_nucPar_verif) ]
+  implicit none
+  BEGIN_DOC
+  ! Jastrow variation  with respect to nuclear parameters
+  ! 
+  ! Dimensions : nucl_num
+  END_DOC
+
+  integer          :: i, j
+  double precision :: eps = 1d-7, der1, der2, pos0
+
+  do j = 1, nucl_num
+    !!!
+    pos0 = eps * jast_pen(j)
+    !!!
+    jast_pen(j) = jast_pen(j) + pos0
+    TOUCH jast_pen
+    der1 = 0.d0
+    !DIR$ LOOP COUNT (100)
+    do i = 1, elec_num
+      der1 += jast_elec_Simple_value(i)
+    end do
+    !!!
+    jast_pen(j) = jast_pen(j) - 2.d0 * pos0
+    TOUCH jast_pen
+    der2 = 0.d0
+    !DIR$ LOOP COUNT (100)
+    do i = 1, elec_num
+      der2 += jast_elec_Simple_value(i)
+    end do
+    !!!
+    J_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0
+    !!!
+  end do
+
+  J_deriv_nucPar_verif_min = min(J_deriv_nucPar_verif_min,minval(J_deriv_nucPar_verif))
+  J_deriv_nucPar_verif_max = max(J_deriv_nucPar_verif_max,maxval(J_deriv_nucPar_verif))
+  SOFT_TOUCH J_deriv_nucPar_verif_min J_deriv_nucPar_verif_max
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER  [ double precision, J_deriv_nucPar_ex, (size_J_deriv_nucPar_ex) ]
+  implicit none
+  BEGIN_DOC
+  ! Jastrow variation  with respect to nuclear parameters
+  ! 
+  ! Dimensions : nucl_num
+  END_DOC
+
+  integer          :: j
+
+  do j = 1, nucl_num
+    J_deriv_nucPar_ex(j) = jast_elec_Simple_deriv_nucPar(j)
+  end do
+
+  J_deriv_nucPar_ex_min = min(J_deriv_nucPar_ex_min,minval(J_deriv_nucPar_ex))
+  J_deriv_nucPar_ex_max = max(J_deriv_nucPar_ex_max,maxval(J_deriv_nucPar_ex))
+  SOFT_TOUCH J_deriv_nucPar_ex_min J_deriv_nucPar_ex_max
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER  [ double precision, J_deriv_bPar_ex ]
+  implicit none
+  BEGIN_DOC
+  ! Jastrow variation with respect to parameter b
+  ! we suppose that jast_b_up_up = jast_b_up_dn
+  END_DOC
+
+  double precision :: tmp, Jtmp, a, b, rij
+  integer          :: i, j
+
+  b = jast_b_up_up ! and also = jast_b_up_dn
+
+  ! parallele spin up-up
+  Jtmp = 0.d0
+  a    = jast_a_up_up
+  do j = 1, elec_alpha_num-1
+    !DIR$ LOOP COUNT (50)
+    do i = j+1, elec_alpha_num
+      rij  = elec_dist(i,j)
+      tmp  = rij/(1.d0+b*rij) 
+      Jtmp = Jtmp + tmp * tmp
+    enddo
+  enddo
+  do j = elec_alpha_num+1, elec_num-1
+    !DIR$ LOOP COUNT (50)
+    do i = j+1, elec_num
+      rij  = elec_dist(i,j)
+      tmp  = rij/(1.d0+b*rij) 
+      Jtmp = Jtmp + tmp * tmp
+    enddo
+  enddo
+  J_deriv_bPar_ex = -1.d0 * a * Jtmp
+
+
+
+  ! anti-parallele spin
+  Jtmp = 0.d0
+  a    = jast_a_up_dn
+  do j = 1, elec_alpha_num
+  !DIR$ LOOP COUNT (50)
+    do i = elec_alpha_num+1, elec_num
+      rij  = elec_dist(i,j)
+      tmp  = rij/(1.d0+b*rij) 
+      Jtmp = Jtmp + tmp * tmp
+    enddo
+  enddo
+  J_deriv_bPar_ex = J_deriv_bPar_ex - a * Jtmp
+
+
+  J_deriv_bPar_ex_min = min(J_deriv_bPar_ex_min, J_deriv_bPar_ex)
+  J_deriv_bPar_ex_max = max(J_deriv_bPar_ex_max, J_deriv_bPar_ex)
+  SOFT_TOUCH J_deriv_bPar_ex_min J_deriv_bPar_ex_max
+
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER  [ double precision, J_deriv_bPar_verif ]
+  implicit none
+  BEGIN_DOC
+  ! verification: Jastrow variation with respect to parameter b
+  ! we suppose that jast_b_up_up = jast_b_up_dn
+  END_DOC
+
+  integer          :: i
+  double precision :: eps = 1d-7, der1, der2, pos0
+
+  pos0 = eps * jast_b_up_up
+  ! !!!
+  jast_b_up_up = jast_b_up_up + pos0
+  TOUCH jast_b_up_up
+  jast_b_up_dn = jast_b_up_up
+  TOUCH jast_b_up_dn
+  der1 = 0.d0
+  !DIR$ LOOP COUNT (100)
+  do i = 1, elec_num
+    der1 += jast_elec_Simple_value(i)
+  end do
+  ! !!!
+  jast_b_up_up = jast_b_up_up - 2.d0 * pos0
+  TOUCH jast_b_up_up
+  jast_b_up_dn = jast_b_up_up
+  TOUCH jast_b_up_dn
+  der2 = 0.d0
+  !DIR$ LOOP COUNT (100)
+  do i = 1, elec_num
+    der2 += jast_elec_Simple_value(i)
+  end do
+  ! !!!
+  J_deriv_bPar_verif = 0.5d0 * (der1 - der2) / pos0
+
+
+  J_deriv_bPar_verif_min = min(J_deriv_bPar_verif_min, J_deriv_bPar_verif)
+  J_deriv_bPar_verif_max = max(J_deriv_bPar_verif_max, J_deriv_bPar_verif)
+  SOFT_TOUCH J_deriv_bPar_verif_min J_deriv_bPar_verif_max
+
+END_PROVIDER

src/PROPERTIES/properties_energy.irp.f

@@ -277,9 +277,3 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
 
 END_PROVIDER
 
-
-
-
-
-
-

src/PROPERTIES/properties_energy_SVD.irp.f

@@ -0,0 +1,49 @@
+!==========================================================================!
+!                             DIMENSIONS                                   !
+!==========================================================================!
+
+BEGIN_PROVIDER [ double precision, E_kin_elec_SVD, (elec_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Electronic Kinetic energy : -1/2 (Lapl.Psi_SVD)/Psi_SVD
+  END_DOC
+  integer :: i
+  do i = 1, elec_num
+    E_kin_elec_SVD(i) = -0.5d0 * psi_lapl_psi_inv_SVD(i)
+  enddo
+END_PROVIDER
+
+
+
+!==========================================================================!
+!                             PROPERTIES                                   !
+!==========================================================================!
+
+
+BEGIN_PROVIDER  [ double precision, E_loc_SVD ]
+  implicit none
+  include '../types.F'
+  BEGIN_DOC
+  ! Local energy : E_kin + E_pot + E_nucl
+  END_DOC
+
+  integer :: i
+
+  E_loc_SVD = E_nucl
+  !DIR$ VECTOR ALIGNED
+  !DIR$ LOOP COUNT(200)
+  do i = 1, elec_num
+    E_loc_SVD += E_kin_elec_SVD(i) + E_pot_elec(i)
+  enddo
+
+  ! Avoid divergence of E_loc_SVD and population explosion
+  if (do_pseudo) then
+    double precision :: delta_e
+    E_loc_SVD = max(2.d0*E_ref, E_loc_SVD)
+  endif
+
+  E_loc_SVD_min = min(E_loc_SVD,E_loc_SVD_min)
+  E_loc_SVD_max = max(E_loc_SVD,E_loc_SVD_max)
+  SOFT_TOUCH E_loc_SVD_min E_loc_SVD_max
+
+END_PROVIDER

src/QMC_SVD/v0/test_SVD.py

@@ -0,0 +1,262 @@
+#!/usr/bin/env python3
+# !!!
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from math import sqrt
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+from scipy.linalg import eig, eigh
+
+# PARAMETERS
+block_time = 20
+eps = 1.
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_energy():
+    buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],  encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+      return None, None
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def set_vmc_params():
+    #subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+    #                                 '-m', 'VMC',
+    #                                 '-l', str(20),
+    #                                 '--time-step=0.3',
+    #                                 '--stop-time=36000',
+    #                                 '--norm=1.e-5',
+    #                                 '-w', '10',
+    #                                 filename])
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'None', '-l', str(block_time), filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_svd(n_svd):
+    # !!!
+    Ci_h_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_h_matrix_svd = results.find('Ci_h_matrix_svd : [ ') + len( 'Ci_h_matrix_svd : [ ' ) 
+    end_Ci_h_matrix_svd = len(results)#results.find('E_loc :') 
+    Ci_h_matrix_svd_buf = results[beg_Ci_h_matrix_svd:end_Ci_h_matrix_svd]
+    Ci_h_matrix_svd_buf = Ci_h_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_h_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_svd')
+            stop
+        else:
+          #Ci_h_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+        # !!!
+    # !!!
+    # Ci_h_matrix_svd = np.reshape(Ci_h_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_h_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_svd(n_svd):
+    # !!!
+    Ci_overlap_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_svd = results.find('Ci_overlap_matrix_svd : [ ') + len( 'Ci_overlap_matrix_svd : [ ' ) 
+    end_Ci_overlap_matrix_svd = len(results)#results.find('Ci_h_matrix_svd : [') 
+    Ci_overlap_matrix_svd_buf = results[beg_Ci_overlap_matrix_svd:end_Ci_overlap_matrix_svd]
+    Ci_overlap_matrix_svd_buf = Ci_overlap_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_overlap_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        # !!!
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_svd')
+            stop
+            # !!!
+        # !!!
+        else:
+          #Ci_overlap_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+          # !!!
+      # !!!
+    # !!!
+    #Ci_overlap_matrix_svd = np.reshape(Ci_overlap_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_postsvd(n_svd):
+    # !!!
+    Ci_h_matrix_postsvd = np.zeros( (n_svd*n_svd*n_svd*n_svd) )
+    # !!!
+    beg_Ci_h_matrix_postsvd = results.find('Ci_h_matrix_postsvd : [ ') + len( 'Ci_h_matrix_postsvd : [ ' ) 
+    end_Ci_h_matrix_postsvd = len(results)
+    Ci_h_matrix_postsvd_buf = results[beg_Ci_h_matrix_postsvd:end_Ci_h_matrix_postsvd]
+    Ci_h_matrix_postsvd_buf = Ci_h_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_h_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_postsvd')
+            stop
+        else:
+            Ci_h_matrix_postsvd[indc-1] = float( line[2] )
+        # !!!
+    Ci_h_matrix_postsvd = np.reshape(Ci_h_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    return(Ci_h_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_postsvd(n_svd):
+    # !!!
+    Ci_overlap_matrix_postsvd = np.zeros( (n_svd*n_svd*n_svd*n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_postsvd = results.find('Ci_overlap_matrix_postsvd : [ ') + len( 'Ci_overlap_matrix_postsvd : [ ' ) 
+    end_Ci_overlap_matrix_postsvd = len(results)
+    Ci_overlap_matrix_postsvd_buf = results[beg_Ci_overlap_matrix_postsvd:end_Ci_overlap_matrix_postsvd]
+    Ci_overlap_matrix_postsvd_buf = Ci_overlap_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_overlap_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_postsvd')
+            stop
+        else:
+            Ci_overlap_matrix_postsvd[indc-1] = float( line[2] )
+        # !!!
+    Ci_overlap_matrix_postsvd = np.reshape(Ci_overlap_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    t0 = time.time()
+    # !!!
+    filename = "/home/ammar/qp2/src/svdwf/h2o.ezfio"
+    ezfio.set_file(filename)
+    # !!!
+    print("Today's date:", datetime.now() )
+    print("filename = {}".format(filename))
+    print("start QMC:")
+    #set_vmc_params()
+    #run_qmc() 
+    #
+    ##stop_qmc()
+    results =  subprocess.check_output(['qmcchem', 'result', filename],  encoding='UTF-8')
+    # !!!
+    #file = open('results.txt','a') 
+    #file.write('\n \n \n')
+    #file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #file.write("Today's date: {}\n".format(datetime.now()) )
+    #file.write("filename = {}\n".format(filename))
+    #file.write('\n')
+    #file.write( results )
+    #file.write('\n')
+    #print("end QMC after {} minutes\n".format((time.time()-t0)/60.) )
+    #file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
+    #file.close()
+    # !!!
+    E_loc, ErrEloc = get_energy()
+    print('Eloc = {} +/- {}'.format(E_loc, ErrEloc))
+    # !!!
+    n_svd                 = 10
+    Ci_h_matrix_svd       = get_Ci_h_matrix_svd(n_svd)
+    Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd(n_svd)
+    # !!!
+    aa = Ci_h_matrix_svd
+    bb = Ci_overlap_matrix_svd
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, check_finite=True, homogeneous_eigvals=False)
+    print( eigvals_svd+9.194966082434476 )
+    # !!!
+    Ci_h_matrix_postsvd       = get_Ci_h_matrix_postsvd(n_svd)
+    Ci_overlap_matrix_postsvd = get_Ci_overlap_matrix_postsvd(n_svd)
+    # !!!
+    aa = Ci_h_matrix_postsvd
+    bb = Ci_overlap_matrix_postsvd
+    eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, check_finite=True, homogeneous_eigvals=False)
+    print( eigvals_postsvd+9.194966082434476 )
+    # !!!
+    print("end code after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!! 
+# !!!
+
+

src/QMC_SVD/v1/QR.py

@@ -0,0 +1,10 @@
+# !!!
+import numpy as np
+# !!!
+def QR_fact(X):
+    Q, R     = np.linalg.qr(X, mode="reduced")
+    D        = np.diag( np.sign( np.diag(R) ) )
+    Qunique  = np.dot(Q,D)
+    #Runique = np.dot(D,R)
+    return(Qunique)
+# !!!

src/QMC_SVD/v1/RSVD.py

@@ -0,0 +1,20 @@
+# !!!
+import numpy as np
+from QR import QR_fact
+# !!!
+def powit_RSVD(X, new_r, nb_powit, nb_oversamp):
+    # !!!
+    G  = np.random.randn(X.shape[1], new_r+nb_oversamp)
+    Q  = QR_fact( np.dot(X,G) )
+    # !!!
+    for _ in range(nb_powit):
+        Q = QR_fact( np.dot(X.T,Q) )
+        Q = QR_fact( np.dot(X,Q) )
+        # !!!
+    Y = np.dot(Q.T,X)
+    # !!!
+    U, S, VT = np.linalg.svd(Y, full_matrices=0)
+    U = np.dot(Q,U)
+    return U[:,:(new_r)], S[:(new_r)], VT[:(new_r),:]
+    # !!!
+# !!!

src/QMC_SVD/v1/diag_after_QMCCHEM.py

@@ -0,0 +1,513 @@
+#!/usr/bin/env python3
+# !!!
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from math import sqrt
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+from scipy.linalg import eig, eigh
+from RSVD import powit_RSVD  
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_energy():
+    buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', EZFIO_file],  encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+      return None, None
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', EZFIO_file])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', EZFIO_file])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def set_vmc_params():
+    #subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+    #                                 '-m', 'VMC',
+    #                                 '-l', str(20),
+    #                                 '--time-step=0.3',
+    #                                 '--stop-time=36000',
+    #                                 '--norm=1.e-5',
+    #                                 '-w', '10',
+    #                                 EZFIO_file])
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'None', '-l', str(block_time), EZFIO_file])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_svd():
+    # !!!
+    Ci_h_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_h_matrix_svd = results.find('Ci_h_matrix_svd : [ ') + len( 'Ci_h_matrix_svd : [ ' ) 
+    end_Ci_h_matrix_svd = len(results)
+    Ci_h_matrix_svd_buf = results[beg_Ci_h_matrix_svd:end_Ci_h_matrix_svd]
+    Ci_h_matrix_svd_buf = Ci_h_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_h_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_svd')
+            stop
+        else:
+          #Ci_h_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+        # !!!
+    # !!!
+    # Ci_h_matrix_svd = np.reshape(Ci_h_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_h_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_svd():
+    # !!!
+    Ci_overlap_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_svd = results.find('Ci_overlap_matrix_svd : [ ') + len( 'Ci_overlap_matrix_svd : [ ' ) 
+    end_Ci_overlap_matrix_svd = len(results) 
+    Ci_overlap_matrix_svd_buf = results[beg_Ci_overlap_matrix_svd:end_Ci_overlap_matrix_svd]
+    Ci_overlap_matrix_svd_buf = Ci_overlap_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_overlap_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        # !!!
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_svd')
+            stop
+            # !!!
+        # !!!
+        else:
+          #Ci_overlap_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+          # !!!
+      # !!!
+    # !!!
+    #Ci_overlap_matrix_svd = np.reshape(Ci_overlap_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_postsvd():
+    #file = open('verif_order.txt','w')
+    # !!!
+    Ci_h_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_h_matrix_postsvd = results.find('Ci_h_matrix_postsvd : [ ') + len( 'Ci_h_matrix_postsvd : [ ' ) 
+    end_Ci_h_matrix_postsvd = len(results)
+    Ci_h_matrix_postsvd_buf = results[beg_Ci_h_matrix_postsvd:end_Ci_h_matrix_postsvd]
+    Ci_h_matrix_postsvd_buf = Ci_h_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_h_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_postsvd')
+            stop
+        else:
+            # !!!
+            kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            #indcrep = kp + (k-1)*n_svd + (lp-1)*n_svd**2 + (l-1)*n_svd**3
+            #file.write( '{:5} {:5} {:5} {:5} {:5} {:5} \n'.format(indc, indc-indcrep, kp, k, lp, l ) )
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_h_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_h_matrix_postsvd[indc-1] = float( line[2] )
+        # !!!
+    #Ci_h_matrix_postsvd = np.reshape(Ci_h_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    #file.close()
+    return(Ci_h_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_postsvd():
+    # !!!
+    Ci_overlap_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_postsvd = results.find('Ci_overlap_matrix_postsvd : [ ') + len( 'Ci_overlap_matrix_postsvd : [ ' ) 
+    end_Ci_overlap_matrix_postsvd = len(results)
+    Ci_overlap_matrix_postsvd_buf = results[beg_Ci_overlap_matrix_postsvd:end_Ci_overlap_matrix_postsvd]
+    Ci_overlap_matrix_postsvd_buf = Ci_overlap_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_overlap_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_postsvd')
+            stop
+        else:
+            # !!!
+            kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_overlap_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_overlap_matrix_postsvd[indc-1] = float( line[2] )
+            # !!!
+        # !!!
+    #Ci_overlap_matrix_postsvd = np.reshape(Ci_overlap_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def check_symmetric(a, tol=1e-3):
+    return np.all(np.abs(a-a.T) < tol)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def save_results_to_resultsQMC():
+    file = open('resultsQMC.txt','a') 
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write("Today's date: {}\n".format(datetime.now()))
+    file.write("EZFIO file = {}\n".format(EZFIO_file))
+    file.write('\n')
+    file.write( results )
+    file.write('\n')
+    file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
+    file.close()
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Esvd():
+    # !!!
+    # read CI_SVD matrices
+    Ci_h_matrix_svd       = get_Ci_h_matrix_svd()
+    Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd()
+    #print( 'Ci_h_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_svd)) )
+    #print( 'Ci_overlap_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_svd)) )
+    # !!!
+    # symmetrise and diagonalise
+    aa = Ci_h_matrix_svd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_svd
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False)
+    #print( eigvals_svd + E_toadd )
+    recouvre_svd = np.abs(psi_svd_coeff @ vr)
+    ind_gssvd = np.argmax(recouvre_svd)
+    # !!!
+    E_svd = eigvals_svd[ind_gssvd] + E_toadd
+    return( E_svd, vr[:,ind_gssvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Epostsvd():
+    # !!!
+    # read CI_postSVD matrices
+    Ci_h_matrix_postsvd       = get_Ci_h_matrix_postsvd()
+    Ci_overlap_matrix_postsvd = get_Ci_overlap_matrix_postsvd()
+    #print( 'Ci_h_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_postsvd)) )
+    #print( 'Ci_overlap_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_postsvd)) )
+    # !!!
+    # symmetrise and diagonalise
+    aa = Ci_h_matrix_postsvd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_postsvd
+    eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False)
+    #print( eigvals_postsvd + E_toadd )
+    d_postsvd = np.diagflat(psi_svd_coeff)
+    d_postsvd = d_postsvd.reshape( (1,n_svd*n_svd) )
+    recouvre_postsvd = np.abs(d_postsvd @ vr)
+    ind_gspostsvd = np.argmax(recouvre_postsvd)
+    #print(recouvre_postsvd, ind_gspostsvd)
+    # !!!
+    E_postsvd = eigvals_postsvd[ind_gspostsvd] + E_toadd
+    # !!!
+    return( E_postsvd,  vr[:,ind_gspostsvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def SVD_postsvd(sigma_postsvd):
+    # !!!
+    print( 'performing new SVD for the post SVD eigenvector:' )
+    # !!!
+    #sigma_postsvd = sigma_postsvd.reshape( (n_svd,n_svd) )
+    sigma_postsvd = sigma_postsvd.reshape( n_svd, n_svd, order='F' )
+    #print( 'sigma_postsvd is symmetric ? {}' .format(check_symmetric(sigma_postsvd)) )
+    # !!!
+    # construct the new matrix Y
+    Y     = U_svd @ sigma_postsvd @ V_svd.T
+    normY = np.linalg.norm(Y, ord='fro')
+    # !!!
+    # parameters of RSVD
+    rank        = n_svd
+    npow        = 10
+    nb_oversamp = 10
+    # !!!
+    # call RSV
+    U_postSVD, sigma_postsvd_diag, VT_postsvd = powit_RSVD(Y, rank, npow, nb_oversamp)
+    # !!!
+    # check precision
+    Y_SVD = np.dot( U_postSVD , np.dot( np.diag(sigma_postsvd_diag) , VT_postsvd ) )
+    energy = np.sum( np.square(sigma_postsvd_diag) ) / normY**2
+    err_SVD = 100. * np.linalg.norm( Y - Y_SVD, ord="fro") / normY
+    print('energy = {}, error = {}\n'.format(energy, err_SVD))
+    # !!!
+    return(U_postSVD, sigma_postsvd_diag, VT_postsvd)
+    # !!!
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_recouv_svd():
+    recouv_svd = np.abs( sigma0 @ sigma_svd )
+    return( recouv_svd )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_recouv_postsvd():
+    d_postsvd        = np.diagflat( sigma0 )
+    d_postsvd        = d_postsvd.reshape( (1, n_svd*n_svd) )
+    recouv_postsvd = np.abs( d_postsvd @ sigma_postsvd )
+    return( recouv_postsvd )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_Hsvd_QP(Hsvd_qp_txt):
+    Hsvd_qp      = np.zeros( (n_svd,n_svd) )
+    Hsvd_qp_file = open(Hsvd_qp_txt, 'r')
+    for line in Hsvd_qp_file:
+        line         = line.split()
+        i            = int(line[0]) - 1
+        j            = int(line[1]) - 1
+        Hsvd_qp[i,j] = float(line[2])
+    return(Hsvd_qp)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_Esvd_QP(Hsvd_qp):
+    # !!!
+    # symmetrise and diagonalise
+    aa = Hsvd_qp
+    aa = 0.5*( aa + aa.T )
+    bb = np.identity(n_svd)
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False) 
+    recouvre_svd = np.abs(psi_svd_coeff @ vr)
+    ind_gssvd = np.argmax(recouvre_svd)
+    E_svd = eigvals_svd[ind_gssvd] + E_toadd
+    return( E_svd, vr[:,ind_gssvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+if __name__ == '__main__':
+    t0 = time.time()
+    # !!!
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    EZFIO_file = "/home/ammar/qp2/src/svdwf/h2o_QPsvd.ezfio"
+    E_toadd    =  9.194966082434476 #6.983610961797779    
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    # !!!
+    ezfio.set_file(EZFIO_file)
+    n_svd         = ezfio.get_spindeterminants_n_svd_coefs()
+    psi_svd_coeff = np.array(ezfio.get_spindeterminants_psi_svd_coefs())
+    U_svd         = np.array(ezfio.get_spindeterminants_psi_svd_alpha())
+    V_svd         = np.array(ezfio.get_spindeterminants_psi_svd_beta())
+    # !!!
+    U_svd = U_svd[0,:,:].T
+    V_svd = V_svd[0,:,:].T
+    # !!!
+    print("Today's date:", datetime.now() )
+    print("EZFIO file = {}".format(EZFIO_file))
+    print("nuclear energy = {}".format(E_toadd) )
+    print("n_svd = {} \n".format(n_svd) )
+    # !!!
+    #set_vmc_params()
+    #run_qmc()
+    #print("start QMC:")
+    #stop_qmc()
+    # !!!
+    print( 'getting QMCCHEM results from {}'.format(EZFIO_file) )
+    results =  subprocess.check_output(['qmcchem', 'result', EZFIO_file],  encoding='UTF-8')
+    # !!!
+    E_loc, ErrEloc = get_energy()
+    print('Eloc = {} +/- {}\n'.format(E_loc, ErrEloc))
+    # !!!
+    save_resultsQMC = input( 'save QMC results from {}? (y/n) '.format(EZFIO_file) )
+    if( save_resultsQMC == 'y' ):
+        print('saving in resultsQMC.txt')
+        save_results_to_resultsQMC()
+    print('\n')
+    # !!!
+    sigma0 = psi_svd_coeff
+    # !!!
+    read_QPsvd = input( 'read QP Hsvd matrix ? (y/n) ')
+    if( read_QPsvd == 'y' ):
+        Hsvd_qp_txt  = input('name of file with QM H_svd matrix:')
+        Hsvd_qp      = get_Hsvd_QP(Hsvd_qp_txt)
+        E_svd_QP, _ = get_Esvd_QP(Hsvd_qp)
+        print('QP SVD enegry = {} \n'.format(E_svd_QP) )
+    # !!!
+    E_svd, sigma_svd = get_Esvd()
+    recouv_svd       = get_recouv_svd()
+    print('QMC=CHEM SVD enegry = {} '.format(E_svd) )
+    print('QMC=CHEM recouvrement SVD ={} \n'.format(recouv_svd) )
+    # !!!
+    E_postsvd, sigma_postsvd = get_Epostsvd()
+    recouv_postsvd           = get_recouv_postsvd()
+    print('QMC=CHEM post SVD energy = {} '.format(E_postsvd) )
+    print('QMC=CHEM recouvrement post SVD ={} \n'.format(recouv_postsvd) )
+    # !!!
+    save_to_EZFIO = input( "modify EZFIO (with svd or postsvd)? "  )
+    # !!!
+    if( save_to_EZFIO == 'svd' ):
+        sigma_SVD_toEZFIO      = np.zeros( ( n_svd, 1) )
+        sigma_SVD_toEZFIO[:,0] = sigma_svd
+        ezfio.set_spindeterminants_psi_svd_coefs( sigma_SVD_toEZFIO )
+    # !!!
+    elif( save_to_EZFIO == 'postsvd' ):
+        # SVD first
+        U_postSVD, sigma_postsvd_diag, V_postSVD = SVD_postsvd(sigma_postsvd)
+        V_postSVD = V_postSVD.T
+        # save next in ezfio file   
+        sigma_postSVD_toEZFIO      = np.zeros( ( n_svd, 1) )
+        U_postSVD_toEZFIO          = np.zeros( ( U_postSVD.shape[0], U_postSVD.shape[1], 1) )
+        V_postSVD_toEZFIO          = np.zeros( ( V_postSVD.shape[0], V_postSVD.shape[1], 1) )
+        sigma_postSVD_toEZFIO[:,0] = sigma_postsvd_diag
+        U_postSVD_toEZFIO[:,:,0]   = U_postSVD  
+        V_postSVD_toEZFIO[:,:,0]   = V_postSVD  
+        #
+        ezfio.set_spindeterminants_psi_svd_alpha( U_postSVD_toEZFIO )        
+        ezfio.set_spindeterminants_psi_svd_coefs( sigma_postSVD_toEZFIO )
+        ezfio.set_spindeterminants_psi_svd_beta( V_postSVD_toEZFIO )
+    # !!!
+    else:
+        print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+        exit()
+    # !!!
+    print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!!
+# !!!
+
+

src/QMC_SVD/v1/save_svd.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+# !!!
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from datetime import datetime
+import time
+
+
+# ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+file_to_save = 'postsvd_1.txt'
+# ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+
+
+if __name__ == '__main__':
+    # !!!
+    if input("are you sure you want before saving in {}? (y/n)".format(file_to_save)) != "y":
+        exit()
+    # !!!
+    t0 = time.time()
+    # !!!
+    filename = "/home/ammar/qp2/src/svdwf/h2o_pseudo_1_postsvd"
+    ezfio.set_file(filename)
+    # !!!
+    n_svd = ezfio.get_spindeterminants_n_svd_coefs()
+    d_svd = ezfio.get_spindeterminants_psi_svd_coefs()
+    u_svd = ezfio.get_spindeterminants_psi_svd_alpha()
+    v_svd = ezfio.get_spindeterminants_psi_svd_beta()
+    # !!!
+    print("Today's date:", datetime.now() )
+    print("EZFIO filename: {}".format(filename))
+    # !!!   
+    file = open(file_to_save, 'a') 
+    #
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write("Today's date: {}\n".format(datetime.now()) )
+    file.write("EZFIO filename = {}\n".format(filename))
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write( 'n_svd :')
+    file.write( str(n_svd) )
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write( 'd_svd :')
+    file.write( str(d_svd) )
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write( 'u_svd :')
+    file.write( str(u_svd) )
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write( 'v_svd :')
+    file.write( str(v_svd) )
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #
+    file.write('\n')
+    file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
+    #   
+    file.close()
+    # !!!
+    print("SVD save in {} after {:.3f} minutes".format( file_to_save, (time.time()-t0)/60.) )
+    # !!! 
+# !!!
+
+

src/det.irp.f

@@ -1287,6 +1287,7 @@ END_PROVIDER
 
     ddet = 0.d0
 
+    ! n_to_do < 2*elec_alpha_num
     if (n_to_do < shiftl(elec_alpha_num,1)) then
 
       do while ( n_to_do > 0 )
@@ -1350,6 +1351,7 @@ END_PROVIDER
   ASSERT (ddet /= 0.d0)
 
   det_alpha_value_curr = ddet
+
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
@@ -1677,8 +1679,8 @@ END_PROVIDER
   DaC = 0.d0
   CDb = 0.d0
 
+  !if( det_num .ne. 1 )  then
   if( det_num < shiftr(det_alpha_num*det_beta_num,2))  then
-
     det_num4 = iand(det_num,not(3))
     !DIR$ VECTOR ALIGNED
     do k=1,det_num4,4
@@ -1696,7 +1698,6 @@ END_PROVIDER
         CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
         CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
         CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*f
-
         if ( ((i2-i1) == 1).and.((i3-i1) == 2).and.((i4-i1) == 3) ) then
           DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
           + det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
@@ -1719,40 +1720,61 @@ END_PROVIDER
         CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
       endif
     enddo
-
     do k=det_num4+1,det_num
       i = det_coef_matrix_rows(k)
       j = det_coef_matrix_columns(k)
       DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
       CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
     enddo
-
   else
 
+    ! DaC(det_beta_num) = Trans( det_coef_matrix_dense(det_alpha_num,det_beta_num) )
+    !                   @ det_alpha_value(det_alpha_num)
     call dgemv('T',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
       size(det_coef_matrix_dense,1), det_alpha_value, 1, 0.d0, DaC, 1)
+    ! CDb(det_alpha_num) = det_coef_matrix_dense(det_alpha_num,det_beta_num)
+    !                    @ det_beta_value(det_beta_num)
     call dgemv('N',det_alpha_num,det_beta_num,1.d0,det_coef_matrix_dense, &
       size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1)
-
   endif
+  !endif
+
 
   ! Value
   ! -----
-
   psidet_value = 0.d0
+  if(det_num .eq. 1) then
+    psidet_value = det_alpha_value_curr * det_beta_value_curr
+  else
     do j=1,det_beta_num
       psidet_value = psidet_value + det_beta_value(j) * DaC(j)
     enddo
-
+  endif
 
   if (psidet_value == 0.d0) then
     call abrt(irp_here,'Determinantal component of the wave function is zero.')
   endif
   psidet_inv = 1.d0/psidet_value
 
+
   ! Gradients
   ! ---------
-
+  if(det_num .eq. 1) then
+    do i = 1, elec_alpha_num
+      psidet_grad_lapl(1,i) = det_alpha_grad_lapl_curr(1,i) * det_beta_value_curr
+      psidet_grad_lapl(2,i) = det_alpha_grad_lapl_curr(2,i) * det_beta_value_curr
+      psidet_grad_lapl(3,i) = det_alpha_grad_lapl_curr(3,i) * det_beta_value_curr
+      psidet_grad_lapl(4,i) = det_alpha_grad_lapl_curr(4,i) * det_beta_value_curr
+    enddo
+    do i = elec_alpha_num+1, elec_num
+      psidet_grad_lapl(1,i) = det_beta_grad_lapl_curr(1,i) * det_alpha_value_curr
+      psidet_grad_lapl(2,i) = det_beta_grad_lapl_curr(2,i) * det_alpha_value_curr
+      psidet_grad_lapl(3,i) = det_beta_grad_lapl_curr(3,i) * det_alpha_value_curr
+      psidet_grad_lapl(4,i) = det_beta_grad_lapl_curr(4,i) * det_alpha_value_curr
+    enddo
+  else
+    ! psidet_grad_lapl(4,elec_alpha_num) =
+    ! det_alpha_grad_lapl(4,elec_alpha_num,det_alpha_num) @ CDb(det_alpha_num)
     call dgemv('N',elec_alpha_num*4,det_alpha_num,1.d0,                &
         det_alpha_grad_lapl,                                           &
         size(det_alpha_grad_lapl,1)*size(det_alpha_grad_lapl,2),       &
@@ -1763,6 +1785,8 @@ END_PROVIDER
           size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2),         &
           DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1)
     endif
+  endif
+
 
   if (do_pseudo) then
     call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv,          &
@@ -1777,6 +1801,467 @@ END_PROVIDER
 
 END_PROVIDER
 
+
+
+
+
+
+
+
+
+
+
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~
+
+ BEGIN_PROVIDER [ integer, n_svd_coefs_unique ]
+   implicit none
+   BEGIN_DOC
+   ! !!!
+   ! rank of Full SVD
+   END_DOC
+   call get_spindeterminants_n_svd_coefs_unique(n_svd_coefs_unique)
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, psi_svd_coefs_unique, (               n_svd_coefs_unique, n_states) ]
+&BEGIN_PROVIDER [ double precision, psi_svd_alpha_unique, (det_alpha_num, n_svd_coefs_unique, n_states) ]
+&BEGIN_PROVIDER [ double precision, psi_svd_beta_unique , (det_beta_num , n_svd_coefs_unique, n_states) ]
+   implicit none
+   BEGIN_DOC
+   ! Full SVD:
+   ! SVD coeff unique
+   ! SVD U unique
+   ! SVD Vt unique
+   END_DOC
+   call get_spindeterminants_psi_svd_coefs_unique(psi_svd_coefs_unique)
+   call get_spindeterminants_psi_svd_alpha_unique(psi_svd_alpha_unique)
+   call get_spindeterminants_psi_svd_beta_unique (psi_svd_beta_unique )
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, det_alpha_value_SVD_unique, (                                  n_svd_coefs_unique) ]
+&BEGIN_PROVIDER [ double precision, det_beta_value_SVD_unique , (                                  n_svd_coefs_unique) ]
+&BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_SVD_unique, (4, elec_alpha_num           , n_svd_coefs_unique) ]
+&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_SVD_unique , (4, elec_alpha_num+1:elec_num, n_svd_coefs_unique) ]
+
+   implicit none
+   BEGIN_DOC
+   ! !!!
+   END_DOC
+
+   integer       :: mm
+   integer, save :: ifirst = 0
+   if (ifirst == 0) then
+     ifirst = 1
+     det_alpha_value_SVD_unique     = 0.d0
+     det_beta_value_SVD_unique      = 0.d0
+     det_alpha_grad_lapl_SVD_unique = 0.d0
+     det_beta_grad_lapl_SVD_unique  = 0.d0
+   endif
+
+   ! det_alpha_value_SVD_unique = psi_svd_alpha_unique.T @ det_alpha_value
+   call dgemv('T', det_alpha_num, n_svd_coefs_unique &
+     , 1.d0, psi_svd_alpha_unique(:,:,1), size(psi_svd_alpha_unique,1), det_alpha_value, 1 &
+     , 0.d0, det_alpha_value_SVD_unique, 1)
+
+   ! det_beta_value_SVD_unique = psi_svd_beta_unique.T @ det_beta_value
+   call dgemv('T', det_beta_num, n_svd_coefs_unique &
+     , 1.d0, psi_svd_beta_unique(:,:,1), size(psi_svd_beta_unique,1), det_beta_value, 1 &
+     , 0.d0, det_beta_value_SVD_unique, 1)
+
+   call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs_unique, det_alpha_num, 1.d0 &
+     , det_alpha_grad_lapl, 4*size(det_alpha_grad_lapl,2)                         &
+     , psi_svd_alpha_unique, size(psi_svd_alpha_unique,1)                         &
+     , 0.d0, det_alpha_grad_lapl_SVD_unique, 4*size(det_alpha_grad_lapl_SVD_unique,2) )
+   if (elec_beta_num /= 0) then
+     call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs_unique, det_beta_num, 1.d0   &
+       , det_beta_grad_lapl, 4*size(det_beta_grad_lapl,2)                   &
+       , psi_svd_beta_unique, size(psi_svd_beta_unique,1)                   &
+       , 0.d0, det_beta_grad_lapl_SVD_unique, 4*size(det_beta_grad_lapl_SVD_unique,2) )
+   endif
+
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ logical, utilise_SVD ]
+&BEGIN_PROVIDER [ integer, n_svd_coefs ]
+   implicit none
+   BEGIN_DOC
+   ! truncated SVD rank
+   END_DOC
+   n_svd_coefs = -1
+   call get_spindeterminants_n_svd_coefs(n_svd_coefs)
+   utilise_SVD = n_svd_coefs > 0
+   if (.not.utilise_SVD) then
+     n_svd_coefs = 1
+   endif
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, n_svd_coefs2 ]
+&BEGIN_PROVIDER [ integer, n_svd_coefs3 ]
+   implicit none
+   BEGIN_DOC
+   ! square and cube of n_svd_coefs
+   END_DOC
+   n_svd_coefs2 = n_svd_coefs * n_svd_coefs
+   n_svd_coefs3 = n_svd_coefs * n_svd_coefs * n_svd_coefs
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, psi_svd_coefs, (               n_svd_coefs, n_states) ]
+&BEGIN_PROVIDER [ double precision, psi_svd_alpha, (det_alpha_num, n_svd_coefs, n_states) ]
+&BEGIN_PROVIDER [ double precision, psi_svd_beta , (det_beta_num , n_svd_coefs, n_states) ]
+   implicit none
+   BEGIN_DOC
+   ! !!!
+   ! truncated SVD
+   END_DOC
+   integer :: l 
+   !do l = 1, n_svd_coefs
+   !  psi_svd_coefs(l,1)   = psi_svd_coefs_unique(l,1)
+   !  psi_svd_alpha(:,l,1) = psi_svd_alpha_unique(:,l,1)
+   !  psi_svd_beta (:,l,1) = psi_svd_beta_unique (:,l,1)
+   !enddo 
+   call get_spindeterminants_psi_svd_coefs(psi_svd_coefs)
+   call get_spindeterminants_psi_svd_alpha(psi_svd_alpha)
+   call get_spindeterminants_psi_svd_beta (psi_svd_beta )
+ END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, det_alpha_value_SVD, (n_svd_coefs) ]
+&BEGIN_PROVIDER [ double precision, det_beta_value_SVD , (n_svd_coefs) ]
+
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! D_alpha_SVD_{k} = sum_i U_{i,k} D_alpha_{i}
+  ! D_beta_SVD_{k } = sum_j V_{j,k} D_beta_{j }
+  ! !!!
+  END_DOC
+
+  integer, save :: ifirst = 0
+  if (ifirst == 0) then
+    ifirst = 1
+    det_alpha_value_SVD = 0.d0
+    det_beta_value_SVD  = 0.d0
+  endif
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  ! det_alpha_value_SVD = psi_svd_alpha.T @ det_alpha_value
+  call dgemv('T', det_alpha_num, n_svd_coefs &
+    , 1.d0, psi_svd_alpha(:,:,1), size(psi_svd_alpha,1), det_alpha_value, 1 &
+    , 0.d0, det_alpha_value_SVD, 1)
+  ! !!!
+  ! det_beta_value_SVD = psi_svd_beta.T @ det_beta_value
+  call dgemv('T', det_beta_num, n_svd_coefs &
+    , 1.d0, psi_svd_beta(:,:,1), size(psi_svd_beta,1), det_beta_value, 1 &
+    , 0.d0, det_beta_value_SVD, 1)
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+  !do l = 1, n_svd_coefs
+  !  tmp = 0.d0
+  !  do ii = 1, det_alpha_num
+  !    tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_value(ii)
+  !  enddo
+  !  det_alpha_value_SVD(l) = tmp
+  !  tmp = 0.d0
+  !  do jj = 1, det_beta_num
+  !    tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_value(jj)
+  !  enddo
+  !  det_beta_value_SVD(l) = tmp
+  !enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, psidet_value_SVD ]
+&BEGIN_PROVIDER [ double precision, psidet_inv_SVD ]
+&BEGIN_PROVIDER [ double precision, det_alpha_pseudo_SVD, (elec_alpha_num, n_svd_coefs) ]
+&BEGIN_PROVIDER [ double precision, det_beta_pseudo_SVD, (elec_alpha_num+1:elec_num, n_svd_coefs) ]
+
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  END_DOC
+
+  integer       :: l
+  integer, save :: ifirst = 0
+  if (ifirst == 0) then
+    ifirst = 1
+    det_alpha_pseudo_SVD = 0.d0
+    det_beta_pseudo_SVD  = 0.d0
+  endif
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  psidet_value_SVD = 0.d0
+  do l = 1, n_svd_coefs
+    psidet_value_SVD = psidet_value_SVD + det_beta_value_SVD(l) &
+                     * psi_svd_coefs(l,1) * det_alpha_value_SVD(l)
+  enddo
+  if (psidet_value_SVD == 0.d0) then
+    call abrt(irp_here,'Determinantal component of the SVD wave function is zero.')
+  endif
+  psidet_inv_SVD = 1.d0 / psidet_value_SVD
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  if (do_pseudo) then
+    ! det_alpha_pseudo_SVD = det_alpha_pseudo @ psi_svd_alpha * psidet_inv_SVD
+    call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs , det_alpha_num, psidet_inv_SVD &
+      , det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(1,1,1) &
+      , size(psi_svd_alpha,1), 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) )
+    ! !!!
+      if (elec_beta_num /= 0) then
+        ! det_beta_pseudo_SVD = det_beta_pseudo @ psi_svd_beta * psidet_inv_SVD
+        call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs , det_beta_num, psidet_inv_SVD &
+          , det_beta_pseudo(elec_alpha_num+1,1), size(det_beta_pseudo,1) &
+          , psi_svd_beta(1,1,1), size(psi_svd_beta,1) &
+          , 0.d0, det_beta_pseudo_SVD(elec_alpha_num+1,1), size(det_beta_pseudo_SVD,1) )
+      endif
+  endif
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+    !do l = 1, n_svd_coefs
+    !  do ee = 1, elec_alpha_num
+    !    tmp = 0.d0
+    !    do ii = 1, det_alpha_num
+    !      tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_pseudo(ee,ii)
+    !    enddo
+    !    det_alpha_pseudo_SVD(ee,l) = tmp
+    !  enddo
+    !  do ee = elec_alpha_num+1, elec_num
+    !    tmp = 0.d0
+    !    do jj = 1, det_beta_num
+    !      tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_pseudo(ee,jj)
+    !    enddo
+    !    det_beta_pseudo_SVD(ee,l) = tmp
+    !  enddo
+    !enddo
+
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_SVD, (4, elec_alpha_num           , n_svd_coefs) ]
+&BEGIN_PROVIDER [ double precision, det_beta_grad_lapl_SVD , (4, elec_alpha_num+1:elec_num, n_svd_coefs) ]
+
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! det_alpha_grad_lapl_SVD_{k} = sum_i U_{i,k} det_alpha_grad_lapl_{i}
+  ! det_beta_grad_lapl_SVD_{k } = sum_j V_{j,k} det_beta_grad_lapl_{j }
+  ! !!!
+  END_DOC
+
+  integer       :: mm
+  integer, save :: ifirst = 0
+  if (ifirst == 0) then
+    ifirst = 1
+    det_alpha_grad_lapl_SVD = 0.d0
+    det_beta_grad_lapl_SVD  = 0.d0
+  endif
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  call dgemm('N', 'N', 4*elec_alpha_num, n_svd_coefs, det_alpha_num, 1.d0 &
+    , det_alpha_grad_lapl, 4*size(det_alpha_grad_lapl,2)          &
+    , psi_svd_alpha, size(psi_svd_alpha,1)                       &
+    , 0.d0, det_alpha_grad_lapl_SVD, 4*size(det_alpha_grad_lapl_SVD,2) )
+  ! !!!
+  if (elec_beta_num /= 0) then
+    call dgemm('N', 'N', 4*elec_beta_num, n_svd_coefs, det_beta_num, 1.d0 &
+      , det_beta_grad_lapl, 4*size(det_beta_grad_lapl,2)          &
+      , psi_svd_beta, size(psi_svd_beta,1)                       &
+      , 0.d0, det_beta_grad_lapl_SVD, 4*size(det_beta_grad_lapl_SVD,2) )
+  endif
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+  !do l = 1, n_svd_coefs
+  !  do mm = 1, 4
+  !  do ee = 1, elec_alpha_num
+  !    tmp = 0.d0
+  !    do ii = 1, det_alpha_num
+  !      tmp = tmp + psi_svd_alpha(ii,l,1) * det_alpha_grad_lapl(mm,ee,ii)
+  !    enddo
+  !    det_alpha_grad_lapl_SVD(mm,ee,l) = tmp
+  !  enddo
+  !  do ee = elec_alpha_num+1, elec_num
+  !    tmp = 0.d0
+  !    do jj = 1, det_beta_num
+  !      tmp = tmp + psi_svd_beta(jj,l,1) * det_beta_grad_lapl(mm,ee,jj)
+  !    enddo
+  !    det_beta_grad_lapl_SVD(mm,ee,l) = tmp
+  !  enddo
+  !  enddo
+  !enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, psidet_grad_lapl_SVD, (4,elec_num) ]
+&BEGIN_PROVIDER [ double precision, pseudo_non_local_SVD, (elec_num) ]
+
+  implicit none
+  BEGIN_DOC
+  ! !!!
+  ! Sigma is diagonal matrix: Sigma = psi_svd_coefs
+  ! !!!
+  ! Gradient of the determinantal part of the wave function after SVD
+  ! Laplacian of determinantal part of the wave function after SVD
+  ! Non-local component of the pseudopotentials after SVD
+  ! !!!
+  ! for each electron:
+  !   for mm = 1, 2, 3 : grad_x Psi, grad_y Psi, grad_z Psi
+  !   for mm = 4: first term (only) of Laplacian Psi
+  ! !!!
+  END_DOC
+
+  integer          :: l, mm, ee
+  integer, save    :: ifirst=0
+  double precision :: tmp
+  if (ifirst == 0) then
+    ifirst = 1
+    psidet_grad_lapl_SVD = 0.d0
+    pseudo_non_local_SVD = 0.d0
+  endif
+
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  do mm = 1, 4
+    ! !!!
+    do ee = 1, elec_alpha_num
+      tmp = 0.d0
+      do l = 1, n_svd_coefs
+        tmp = tmp + det_alpha_grad_lapl_SVD(mm,ee,l) * psi_svd_coefs(l,1) &
+                  * det_beta_value_SVD(l)
+      enddo
+      psidet_grad_lapl_SVD(mm,ee) = tmp
+    enddo
+    ! !!!
+    if (elec_beta_num /= 0) then
+      do ee = elec_alpha_num+1, elec_num
+      tmp = 0.d0
+        do l = 1, n_svd_coefs
+          tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) &
+                    * det_beta_grad_lapl_SVD(mm,ee,l)
+        enddo
+        psidet_grad_lapl_SVD(mm,ee) = tmp
+      enddo
+    endif
+    ! !!!
+  enddo
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+  ! !!!
+  if (do_pseudo) then
+    ! !!!
+    do ee = 1, elec_alpha_num
+      tmp = 0.d0
+      do l = 1, n_svd_coefs
+        tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l)
+      enddo
+      pseudo_non_local_SVD(ee) = tmp * psidet_inv_SVD
+    enddo
+    ! !!!
+    if (elec_beta_num /= 0) then
+      do ee = elec_alpha_num+1, elec_num
+        tmp = 0.d0
+        do l = 1, n_svd_coefs
+          tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) * det_beta_pseudo_SVD(ee,l)
+        enddo
+        pseudo_non_local_SVD(ee) = tmp * psidet_inv_SVD
+      enddo
+    endif
+    ! !!!
+  endif
+  ! !!!
+  ! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
+
+END_PROVIDER
+
+
+
+
+! ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~ ~~ ~
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
+
+
+
+
+!____________________________________________________________________________________________________________________
+
+
+
+ BEGIN_PROVIDER [ integer, n_svd_selected ]
+&BEGIN_PROVIDER [ integer, n_svd_toselect ]
+   implicit none
+   BEGIN_DOC
+   ! !!!
+   ! n_svd_selected : rank of selected space    ( ~ n_svd x n_svd )
+   ! n_svd_toselect : rank space to select from ( ~ na x nb - n_svd x n_svd_toselect )
+   ! !!!
+   END_DOC
+   call get_spindeterminants_n_svd_selected(n_svd_selected)
+   call get_spindeterminants_n_svd_toselect(n_svd_toselect)
+ END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, psi_svd_alpha_numselected, (n_svd_selected , n_states) ]
+&BEGIN_PROVIDER [ integer, psi_svd_beta_numselected , (n_svd_coefs    , n_states) ]
+&BEGIN_PROVIDER [ integer, psi_svd_alpha_numtoselect, (n_svd_toselect, n_states) ]
+&BEGIN_PROVIDER [ integer, psi_svd_beta_numtoselect , (n_svd_toselect, n_states) ]
+   implicit none
+   BEGIN_DOC
+   ! !!!
+   ! pairs of integers indicating the number of the vectors U and V
+   ! !!!
+   END_DOC
+   call get_spindeterminants_psi_svd_alpha_numselected(psi_svd_alpha_numselected)
+   call get_spindeterminants_psi_svd_beta_numselected (psi_svd_beta_numselected )
+   call get_spindeterminants_psi_svd_alpha_numtoselect(psi_svd_alpha_numtoselect)
+   call get_spindeterminants_psi_svd_beta_numtoselect (psi_svd_beta_numtoselect )
+ END_PROVIDER
+
+! BEGIN_PROVIDER [ double precision, psi_svd_alpha_toselect, (det_alpha_num, n_svd_toselect, n_states) ]
+!&BEGIN_PROVIDER [ double precision, psi_svd_beta_toselect , (det_beta_num , n_svd_toselect, n_states) ]
+!   implicit none
+!   BEGIN_DOC
+!   ! !!!
+!   ! pair of | u_i v_j > to select from SVD
+!   ! !!!
+!   END_DOC
+!   integer :: l, i, j
+!   do l = 1, n_svd_toselect
+!     i = psi_svd_alpha_numtoselect(l,1)
+!     j = psi_svd_beta_numtoselect (l,1)
+!     psi_svd_alpha_toselect(:,l,1) = psi_svd_alpha_unique(:,i,1)
+!     psi_svd_beta_toselect (:,l,1) = psi_svd_beta_unique (:,j,1)
+!   enddo
+! END_PROVIDER
+
+!____________________________________________________________________________________________________________________
+
+
+
+
+
+
+
+
 BEGIN_PROVIDER  [ double precision, det_alpha_pseudo_curr, (elec_alpha_num) ]
   implicit none
   BEGIN_DOC
@@ -1856,7 +2341,8 @@ BEGIN_PROVIDER  [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
 !      imo  = mo_list_alpha_curr(j)
 !      do i=1,elec_alpha_num
 !        do k=1,4
-!          det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j)
+!          det_alpha_grad_lapl_curr(k,i) = det_alpha_grad_lapl_curr(k,i) + &
+!                 mo_grad_lapl_alpha(k,i,imo)*slater_matrix_alpha_inv_det(i,j)
 !        enddo
 !      enddo
 !    enddo
@@ -2037,3 +2523,9 @@ END_PROVIDER
   enddo
 END_PROVIDER
 
+
+
+
+
+
+

src/ezfio_interface.irp.f

@@ -48,7 +48,20 @@ data = [ \
 ("simulation_e_trial"              , "double precision" , ""                   ),
 ("simulation_do_run"               , "logical       " , ""                   ),
 ("pseudo_do_pseudo"         , "logical       " , ""                   ),
-
+("spindeterminants_n_svd_coefs_unique", "integer", ""),
+("spindeterminants_n_svd_coefs"       , "integer", ""),
+("spindeterminants_n_svd_selected"    , "integer", ""),
+("spindeterminants_n_svd_toselect"    , "integer", ""),
+("spindeterminants_psi_svd_alpha_unique", "double precision", "(det_alpha_num,n_svd_coefs_unique,n_states)"),
+("spindeterminants_psi_svd_beta_unique" , "double precision", "(det_beta_num,n_svd_coefs_unique,n_states)"),
+("spindeterminants_psi_svd_coefs_unique", "double precision", "(n_svd_coefs_unique,n_states)"),
+("spindeterminants_psi_svd_alpha", "double precision", "(det_alpha_num,n_svd_coefs,n_states)"),
+("spindeterminants_psi_svd_beta" , "double precision", "(det_beta_num,n_svd_coefs,n_states)"),
+("spindeterminants_psi_svd_coefs", "double precision", "(n_svd_coefs,n_states)"),
+("spindeterminants_psi_svd_alpha_numselected" , "integer", "(n_svd_selected,n_states)"),
+("spindeterminants_psi_svd_beta_numselected"  , "integer", "(n_svd_selected,n_states)"),
+("spindeterminants_psi_svd_alpha_numtoselect" , "integer", "(n_svd_toselect,n_states)"),
+("spindeterminants_psi_svd_beta_numtoselect"  , "integer", "(n_svd_toselect,n_states)"),
 ]
 
 data_no_set = [\

src/opt_Jast/opt_jast.py

@@ -0,0 +1,253 @@
+#!/usr/bin/env python3
+
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+from ezfio import ezfio
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+import atexit
+import scipy as sp
+import scipy.optimize
+from math import sqrt
+
+#------------------------------------------------------------------------------
+def make_atom_map():
+    labels    = {}
+    dimension = 0
+    # i: label   of nuclei
+    # k: counter of nuclei
+    for i,k in enumerate(ezfio.nuclei_nucl_label):
+        if k in labels:                
+            labels[k].append(i)            
+        else:
+            labels[k] = [dimension, i]
+            dimension += 1
+    atom_map = [[] for i in range(dimension)]
+    for atom in labels.keys():
+        l = labels[atom]
+        atom_map[l[0]] = l[1:]
+    return atom_map
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_pen():
+    d = ezfio.jastrow_jast_pen
+    return np.array([d[m[0]] for m in atom_map])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_pen(x):
+    #x = np.abs(x)
+    y = list(ezfio.jastrow_jast_pen)
+    for i,m in enumerate(atom_map):
+        for j in m:
+            y[j] = x[i]
+    ezfio.set_jastrow_jast_pen(y)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_b():
+    b = ezfio.get_jastrow_jast_b_up_up()
+    return b
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_b(b):
+    ezfio.set_jastrow_jast_b_up_up(b)
+    ezfio.set_jastrow_jast_b_up_dn(b)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_energy():
+    buffer = subprocess.check_output( 
+            ['qmcchem', 'result', '-e', 'e_loc', EZFIO_file], encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+        return None, None
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_variance():
+    buffer = subprocess.check_output( 
+       ['qmcchem', 'result', '-e', 'e_loc_qmcvar', EZFIO_file], encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, variance, error = [float(x) for x in buffer.split()]
+        return variance, error
+    else:
+        return None, None
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_vmc_params(block_time,total_time):
+    #subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+    #                         '-m', 'VMC',
+    #                         '-l', str(block_time),
+    #                         '--time-step=0.3',
+    #                         '--stop-time=36000',
+    #                         '--norm=1.e-5',
+    #                         '-w', '10',
+    #                         EZFIO_file])
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+                             '-t', str(total_time),
+                             '-l', str(block_time), EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def f(x):
+    # !!!
+    print('f on:')
+    print('nuc param Jast = {}'.format(x[:-1]))
+    print('b   param Jast = {}'.format(x[-1]))
+    # !!!
+    global i_fev
+    i_fev = i_fev + 1
+    # !!!
+    global memo_energy
+    # !!!
+    print ("x = %s"%str(x))
+    h = str(x)
+    if h in memo_energy:
+        return memo_energy[h]
+    # !!!
+    set_params_pen(x[:-1])
+    set_params_b(x[-1])
+    set_params_pen(x)
+    set_vmc_params(block_time_f, total_time_f)
+    # !!!
+    pid = os.fork()
+    if pid == 0:
+        run_qmc()
+        os._exit(os.EX_OK)
+    else:
+        atexit.register(stop_qmc)
+        time.sleep(3.*block_time_f/4.)
+
+        local_thresh = thresh
+        err          = thresh+1.
+        ii           = 1
+        ii_max       = int(total_time_f/block_time_f) - 1
+        while( (err > local_thresh) and (ii<ii_max) ):
+            time.sleep(block_time_f)
+            e, e_err        = get_energy()
+            variance, v_err = get_variance()
+            if e is None or variance is None:
+                continue
+            energy = e + variance # minimise energy + variance
+            err    = sqrt(e_err*e_err+v_err*v_err)
+            print(" energy: %f  %f "%(e, e_err))
+            print(" varian: %f  %f "%(variance, v_err))
+            print(" additi: %f  %f "%(energy, err))
+            print(" ")
+            if (energy-2.*err) > memo_energy['fmin']+thresh:
+               local_thresh = 10.*thresh
+            elif (energy+2.*err) < memo_energy['fmin']-thresh:
+               local_thresh = 10.*thresh
+
+            ii = ii + 1
+            # Check if PID is still running
+            try:
+                os.kill(pid,0)
+            except OSError:
+                print("---")
+                break
+        stop_qmc()
+    # !!!
+    os.wait()
+    memo_energy[h]      = energy + err
+    memo_energy['fmin'] = min(energy, memo_energy['fmin'])
+    # !!!
+    return energy
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    t0 = time.time()
+    # !!!
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    EZFIO_file = "/home/ammar/qp2/src/svdwf/h2o_optJast"
+    # nuclear energy
+    E_toadd = 9.194966082434476
+    # PARAMETERS
+    thresh = 1.e-2
+    # maximum allowed number of function evaluations
+    N_fev = 4
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    # !!!
+    ezfio.set_file(EZFIO_file)
+    print("Today's date:", datetime.now() )
+    print("EZFIO file = {}".format(EZFIO_file))
+    # !!!
+    # map nuclei to a list
+    # for H2O this will give: atom_map = [[0], [1, 2]]
+    atom_map = make_atom_map()
+    n_nucpar = len(atom_map)  # nb of nclear parameters
+    # !!!
+    # x = get_params_pen()
+    x = [1.25816521, 0.23766714]
+    print('initial pen: {}'.format(x))
+    b_par = get_params_b()
+    print('initial b: {}'.format(b_par))
+    x.append(b_par)
+    # !!!
+    i_fev = 0
+    bnds  = [(0.001, 9.99) for _ in range(n_nucpar)]
+    memo_energy = {'fmin': 100000000.}
+    opt = sp.optimize.minimize(f, x, method="Powell", bounds=bnds
+                                , options= {'disp':True, 
+                                            'ftol':thresh, 
+                                            'xtol':0.2, 
+                                            'maxfev':N_fev} )
+    print("x = "+str(opt))
+    set_params_pen(opt['x'])
+    print('number of function evaluations = {}'.format(i_fev))
+    # !!!
+    print('memo_energy: {}'.format(memo_energy))
+    # !!!
+    print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!!

src/opt_Jast/opt_jast_Anthony.py

@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+
+import scipy as sp
+import scipy.optimize
+import numpy as np
+import sys
+import os
+import time
+import subprocess
+from math import sqrt
+
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+
+from ezfio import ezfio
+
+# PARAMETERS
+thresh = 1.e-2
+block_time = 20
+
+def main():
+    if len(sys.argv) != 2:
+      print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
+      sys.exti(1)
+
+    filename = sys.argv[1]
+    ezfio.set_file(filename)
+
+
+    def make_atom_map():
+        labels = {}
+        dimension = 0
+        for i,k in enumerate(ezfio.nuclei_nucl_label):
+            if k in labels:
+                labels[k].append(i)
+            else:
+                labels[k] = [dimension, i]
+                dimension += 1
+        atom_map = [[] for i in range(dimension)]
+        for atom in labels.keys():
+            l = labels[atom]
+            atom_map[l[0]] = l[1:]
+        return atom_map
+    atom_map = make_atom_map()
+
+
+    def get_params_pen():
+        d = ezfio.jastrow_jast_pen
+        return np.array([d[m[0]] for m in atom_map])
+
+
+    def get_energy():
+        buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],
+                                         encoding='UTF-8')
+        if buffer.strip() != "":
+            buffer = buffer.splitlines()[-1]
+            _, energy, error = [float(x) for x in buffer.split()]
+            return energy, error
+        else:
+            return None, None
+
+    def get_variance():
+        buffer = subprocess.check_output(['qmcchem', 'result', '-e',
+                                          'e_loc_qmcvar', filename],
+                                         encoding='UTF-8')
+        if buffer.strip() != "":
+            buffer = buffer.splitlines()[-1]
+            _, variance, error = [float(x) for x in buffer.split()]
+            return variance, error
+        else:
+            return None, None
+
+
+    def set_params_pen(x):
+        x = np.abs(x)
+        y=list(ezfio.jastrow_jast_pen)
+        for i,m in enumerate(atom_map):
+            for j in m:
+                y[j] = x[i]
+        ezfio.set_jastrow_jast_pen(y)
+
+
+    def run_qmc():
+        return subprocess.check_output(['qmcchem', 'run', filename])
+
+
+    def stop_qmc():
+        subprocess.check_output(['qmcchem', 'stop', filename])
+
+
+    def set_vmc_params():
+#        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+#                                 '-m', 'VMC',
+#                                 '-l', str(block_time),
+#                                 '--time-step=0.3',
+#                                 '--stop-time=36000',
+#                                 '--norm=1.e-5',
+#                                 '-w', '10',
+#                                 filename])
+        subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+                                 '-l', str(block_time),
+                                 filename])
+
+    memo_energy = {'fmin': 100000000.}
+    def f(x):
+        print ("x = %s"%str(x))
+        h = str(x)
+        if h in memo_energy:
+            return memo_energy[h]
+        set_params_pen(x)
+        set_vmc_params()
+        pid = os.fork()
+        if pid == 0:
+            run_qmc()
+            os._exit(os.EX_OK)
+        else:
+            import atexit
+            atexit.register(stop_qmc)
+
+            err = thresh+1.
+            time.sleep(3.*block_time/4.)
+            local_thresh = thresh
+            while err > local_thresh:
+                time.sleep(block_time)
+                e, e_err = get_energy()
+                variance, v_err = get_variance()
+                if e is None or variance is None:
+                    continue
+                energy = e + variance
+                err = sqrt(e_err*e_err+v_err*v_err)
+                print(" %f  %f    %f  %f  %f  %f"%(e, e_err, variance, v_err, energy, err))
+                if (energy-2.*err) > memo_energy['fmin']+thresh:
+                   local_thresh = 10.*thresh
+                elif (energy+2.*err) < memo_energy['fmin']-thresh:
+                   local_thresh = 10.*thresh
+
+                # Check if PID is still running
+                try:
+                    os.kill(pid,0)
+                except OSError:
+                    print("---")
+                    break
+            stop_qmc()
+        os.wait()
+        memo_energy[h] = energy + err
+        memo_energy['fmin'] = min(energy, memo_energy['fmin'])
+        return energy
+
+
+    def run():
+        x = get_params_pen()
+        if sum(x) == 0.:
+          jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
+          x += jast_a_up_dn
+        opt = sp.optimize.minimize(f,x,method="Powell",
+                             options= {'disp':True, 'ftol':thresh,'xtol':0.02})
+        print("x = "+str(opt))
+        set_params_pen(opt['x'])
+
+    run()
+
+
+
+if __name__ == '__main__':
+    main()
+
+
+
+
+

src/opt_Jast/opt_jast_freegrad.py

@@ -0,0 +1,233 @@
+
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+from ezfio import ezfio
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+import atexit
+import scipy as sp
+import scipy.optimize
+from math import sqrt
+from modif_powell_imp import my_fmin_powell
+
+#------------------------------------------------------------------------------
+def make_atom_map():
+    labels    = {}
+    dimension = 0
+    # i: label   of nuclei
+    # k: counter of nuclei
+    for i,k in enumerate(ezfio.nuclei_nucl_label):
+        if k in labels:                
+            labels[k].append(i)            
+        else:
+            labels[k] = [dimension, i]
+            dimension += 1
+    atom_map = [[] for i in range(dimension)]
+    for atom in labels.keys():
+        l = labels[atom]
+        atom_map[l[0]] = l[1:]
+    return atom_map
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_pen():
+    d = ezfio.jastrow_jast_pen
+    return np.array([d[m[0]] for m in atom_map])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_pen(x):
+    y = list(ezfio.jastrow_jast_pen)
+    for i,m in enumerate(atom_map):
+        for j in m:
+            y[j] = x[i]
+    ezfio.set_jastrow_jast_pen(y)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_b():
+    b = ezfio.get_jastrow_jast_b_up_up()
+    return b
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_b(b):
+    ezfio.set_jastrow_jast_b_up_up(b)
+    ezfio.set_jastrow_jast_b_up_dn(b)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_energy():
+    buffer = subprocess.check_output( 
+            ['qmcchem', 'result', '-e', 'e_loc', EZFIO_file], encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+        return None, None
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_variance():
+    buffer = subprocess.check_output( 
+       ['qmcchem', 'result', '-e', 'e_loc_qmcvar', EZFIO_file], encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, variance, error = [float(x) for x in buffer.split()]
+        return variance, error
+    else:
+        return None, None
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_vmc_params(block_time,total_time):
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+                             '-t', str(total_time),
+                             '-l', str(block_time), EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def f(x):
+    # !!!
+    global i_fev
+    global memo_energy
+    print(' eval {} of f on:'.format(i_fev))
+    print(' nuc param Jast = {}'.format(x[:-1]))
+    print(' b   param Jast = {}'.format(x[-1]))
+    h = str(x)
+    if h in memo_energy:
+        return memo_energy[h]
+    # !!!
+    i_fev = i_fev + 1
+    # !!!
+    set_params_pen(x[:-1])
+    set_params_b(x[-1])
+    block_time_f = 45
+    total_time_f = 180
+    set_vmc_params(block_time_f, total_time_f)
+    # !!!
+    loc_err = 10.
+    ii      = 0
+    ii_max  = 5 
+    energy  = None
+    err     = None
+    while( thresh < loc_err ):
+        run_qmc()
+        energy, err = get_energy()
+        if( (energy is None) or (err is None) ):
+            continue
+        elif( memo_energy['fmin'] < (energy-2.*err) ):
+            print(" %d energy: %f  %f "%(ii, energy, err))
+            break
+        else:
+            loc_err = err
+            ii      = ii + 1
+            print(" %d energy: %f  %f "%(ii, energy, err))
+            if( ii_max < ii ):
+                break
+    print(" ")
+    # !!!
+    memo_energy[h]      = energy + err
+    memo_energy['fmin'] = min(energy, memo_energy['fmin'])
+    # !!!
+    return energy
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    # !!!
+    t0 = time.time()
+    # !!!
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    EZFIO_file = "/home/aammar/qp2/src/svdwf/h2o_optJast"
+    # PARAMETERS
+    thresh = 1.e-2
+    # maximum allowed number of function evaluations
+    N_fev = 50
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    # !!!
+    ezfio.set_file(EZFIO_file)
+    print(" Today's date:", datetime.now() )
+    print(" EZFIO file = {}".format(EZFIO_file))
+    # !!!
+    # map nuclei to a list
+    # for H2O this will give: atom_map = [[0], [1, 2]]
+    atom_map = make_atom_map()
+    n_par = len(atom_map)  # nb of nclear parameters
+    n_par = n_par + 1      # e-e parameter b
+    # !!!
+    # x = get_params_pen()
+    x = [1.29386006, 0.21362821]
+    print(' initial pen: {}'.format(x))
+    #b_par = get_params_b()
+    b_par = 1.5291090863304375
+    print(' initial   b: {}'.format(b_par))
+    x.append(b_par)
+    # !!!
+    i_fev = 1
+    memo_energy = {'fmin': 100.}
+    # !!!
+    #bnds  = [(0.001, 9.99) for _ in range(n_par)]
+    #opt = sp.optimize.minimize(f, x, method="Newton-CG", bounds=bnds
+    #                            , options= {'disp':True} )
+    # !!!
+    x_min = [ (0.001) for _ in range(n_par) ] 
+    x_max = [ (9.999) for _ in range(n_par) ]
+    opt   = my_fmin_powell( f, x, x_min, x_max
+			     , xtol        = 0.02
+			     , ftol        = thresh 
+		             , maxfev      = N_fev
+			     , full_output = 1
+                             , verbose     = 1 )
+    # !!!
+    print(" x = "+str(opt))
+    #set_params_pen(opt['x'])
+    print(' number of function evaluations = {}'.format(i_fev))
+    # !!!
+    print(' memo_energy: {}'.format(memo_energy))
+    # !!!
+    print(" end after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!!

src/opt_Jast/opt_jast_grad.py

@@ -0,0 +1,301 @@
+#!/usr/bin/env python3
+
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+from ezfio import ezfio
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+import atexit
+import scipy as sp
+import scipy.optimize
+from math import sqrt
+
+#------------------------------------------------------------------------------
+def make_atom_map():
+    labels    = {}
+    dimension = 0
+    # i: label   of nuclei
+    # k: counter of nuclei
+    for i,k in enumerate(ezfio.nuclei_nucl_label):
+        if k in labels:                
+            labels[k].append(i)            
+        else:
+            labels[k] = [dimension, i]
+            dimension += 1
+    atom_map = [[] for i in range(dimension)]
+    for atom in labels.keys():
+        l = labels[atom]
+        atom_map[l[0]] = l[1:]
+    return atom_map
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_pen():
+    d = ezfio.jastrow_jast_pen
+    return np.array([d[m[0]] for m in atom_map])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_pen(x):
+    y = list(ezfio.jastrow_jast_pen)
+    for i,m in enumerate(atom_map):
+        for j in m:
+            y[j] = x[i]
+    ezfio.set_jastrow_jast_pen(y)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_params_b():
+    b = ezfio.get_jastrow_jast_b_up_up()
+    return b
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_params_b(b):
+    ezfio.set_jastrow_jast_b_up_up(b)
+    ezfio.set_jastrow_jast_b_up_dn(b)
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_energy():
+    buffer = subprocess.check_output( 
+            ['qmcchem', 'result', '-e', 'e_loc', EZFIO_file], encoding='UTF-8')
+    buffer = buffer.splitlines()[-1]
+    _, energy, error = [float(x) for x in buffer.split()]
+    return energy, error
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_variance():
+    buffer = subprocess.check_output( 
+       ['qmcchem', 'result', '-e', 'e_loc_qmcvar', EZFIO_file], encoding='UTF-8')
+    buffer = buffer.splitlines()[-1]
+    _, variance, error = [float(x) for x in buffer.split()]
+    return variance, error
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_energy_deriv_nucPar():
+    # !!!
+    buffer = subprocess.check_output(
+             ['qmcchem', 'result', EZFIO_file], encoding='UTF-8')
+    end    = len(buffer)
+    # !!!
+    deriv_nucPar1 = []
+    beg           = buffer.find('E_deriv_nucpar_loc1 : [ ') 
+    beg           = beg +  len( 'E_deriv_nucpar_loc1 : [ ' )
+    der_buf       = buffer[beg:end]
+    der_buf       = der_buf.split( '\n' )
+    for iline in range(1, n_nucpar+1):
+        line = der_buf[iline].split()
+        deriv_nucPar1.append( float(line[2]) )
+    # !!!
+    deriv_nucPar2 = []
+    beg           = buffer.find('E_deriv_nucpar_loc2 : [ ') 
+    beg           = beg +  len( 'E_deriv_nucpar_loc2 : [ ' )
+    der_buf       = buffer[beg:end]
+    der_buf       = der_buf.split( '\n' )
+    for iline in range(1, n_nucpar+1):
+        line          = der_buf[iline].split()
+        deriv_nucPar2.append( float(line[2]) )
+    # !!!
+    deriv_nucPar = []
+    e, _ = get_energy()
+    for i in range(n_nucpar):
+        tmp = 2. * ( deriv_nucPar1[i] - e * deriv_nucPar2[i] )
+        deriv_nucPar.append( tmp )
+    # !!!
+    return deriv_nucPar
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def get_energy_deriv_bPar():
+    # !!!
+    buffer = subprocess.check_output(
+             ['qmcchem', 'result', EZFIO_file], encoding='UTF-8')
+    end    = len(buffer)
+    # !!!
+    beg         = buffer.find('E_deriv_bpar_loc1 : ') 
+    beg         = beg +  len( 'E_deriv_bpar_loc1 : ' )
+    der_buf     = buffer[beg:end]
+    der_buf     = der_buf.split( '\n' )
+    der_buf     = der_buf[0].split()
+    deriv_bPar1 = float(der_buf[0])
+    # !!!
+    beg         = buffer.find('E_deriv_bpar_loc2 : ') 
+    beg         = beg +  len( 'E_deriv_bpar_loc2 : ' )
+    der_buf     = buffer[beg:end]
+    der_buf     = der_buf.split( '\n' )
+    der_buf     = der_buf[0].split()
+    deriv_bPar2 = float(der_buf[0])
+    # !!!
+    e, _       = get_energy()
+    deriv_bPar = deriv_bPar1 - e * deriv_bPar2
+    # !!!
+    return deriv_bPar
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def set_vmc_params(block_time,total_time):
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+                             '-t', str(total_time),
+                             '-l', str(block_time), EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', EZFIO_file])
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def f(x):
+    # !!!
+    print('f on:')
+    print('nuc param Jast = {}'.format(x[:-1]))
+    print('b   param Jast = {}'.format(x[-1]))
+    # !!!
+    global i_fev
+    i_fev = i_fev + 1
+    # !!!
+    set_params_pen(x[:-1])
+    set_params_b(x[-1])
+    block_time_f = 20
+    total_time_f = 40
+    set_vmc_params(block_time_f, total_time_f)
+    # !!!
+    err    = 10.
+    ii     = 0
+    ii_max = 5
+    while( (err > thresh) and (ii<ii_max) ):
+        ii = ii + 1
+        run_qmc()
+        energy, err = get_energy()
+        print(" energy: %f  %f "%(energy, err))
+    print(" ")
+    # !!!
+    global memo_energy
+    h                   = str(x)
+    memo_energy[h]      = energy + err
+    memo_energy['fmin'] = min(energy, memo_energy['fmin'])
+    # !!!
+    return energy
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+def fprime(x):
+    # !!!
+    # used with CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, 
+    # trust-krylov, trust-exact, trust-constr
+    # !!!
+    print('derivative on:')
+    print('nuc param Jast = {}'.format(x[:-1]))
+    print('b   param Jast = {}'.format(x[-1]))
+    # !!!
+    set_params_pen(x[:-1])
+    set_params_b(x[-1])
+    block_time_fp = 20
+    total_time_fp = 30
+    set_vmc_params(block_time_fp, total_time_fp)
+    run_qmc()
+    # !!!
+    fp = get_energy_deriv_nucPar()
+    print('nuc param Jast = {}'.format(fp))
+    fp_bpar = get_energy_deriv_bPar()
+    print('b   param Jast = {}'.format(fp_bpar))
+    fp.append(fp_bpar)
+    print(" ")
+    # !!!
+    return fp
+#------------------------------------------------------------------------------
+##
+###
+##
+#------------------------------------------------------------------------------
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    # !!!
+    t0 = time.time()
+    # !!!
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    EZFIO_file = "/home/ammar/qp2/src/svdwf/h2o_optJast"
+    # nuclear energy
+    E_toadd = 9.194966082434476
+    # PARAMETERS
+    thresh = 1.e-2
+    # maximum allowed number of function evaluations
+    N_fev = 5
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    # !!!
+    ezfio.set_file(EZFIO_file)
+    print("Today's date:", datetime.now() )
+    print("EZFIO file = {}".format(EZFIO_file))
+    # !!!
+    # map nuclei to a list
+    # for H2O this will give: atom_map = [[0], [1, 2]]
+    atom_map = make_atom_map()
+    n_nucpar = len(atom_map)  # nb of nclear parameters
+    # !!!
+    # x = get_params_pen()
+    x = [1.25816521, 0.35]
+    print('initial pen: {}'.format(x))
+    b_par = get_params_b()
+    print('initial b: {}'.format(b_par))
+    x.append(b_par)
+    # !!!
+    #i_fev = 0
+    #bnds  = [(0.001, 9.99) for _ in range(n_nucpar)]
+    #memo_energy = {'fmin': 100000000.}
+    #opt = sp.optimize.minimize(f, x, method="Powell", bounds=bnds
+     #                           , options= {'disp':True, 
+     #                                       'ftol':thresh, 
+     #                                       'xtol':0.2, 
+     #                                       'maxfev':N_fev} )
+    #print("x = "+str(opt))
+    #set_params_pen(opt['x'])
+    #print('number of function evaluations = {}'.format(i_fev))
+    # !!!
+    #print('memo_energy: {}'.format(memo_energy))
+    # !!!
+    #stop_qmc()
+    print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!!

src/prepare_walkers.irp.f

@@ -15,8 +15,8 @@ subroutine draw_init_points
   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
+    do i = 1, elec_alpha_num
       norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
     enddo
   enddo

src/psi.irp.f

@@ -4,10 +4,28 @@
   ! Value of the wave function
   END_DOC
   
+  double precision :: norm_Err, psi_value2
+
+  if (utilise_svd) then
+    psi_value = psidet_value_svd*jast_value
+  else
     psi_value = psidet_value*jast_value
+  endif
+
   if (psi_value == 0.d0) then
     call abrt(irp_here,"Value of the wave function is 0.")
   endif
+
+  !psi_value2 = psidet_value_svd * jast_value
+  !norm_Err   = (psi_value2 - psi_value)**2 /  psi_value**2
+  !if (norm_Err > 1.d-6) then
+  !   print *, 'probleme dans PROVIDER [ psi_value ]:  norm_Err = ', norm_Err
+     !print *, psi_value2
+     !print *, psi_value
+     !print *, irp_here
+     !stop
+  !endif 
+
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, psi_value_inv ]
@@ -32,8 +50,19 @@ BEGIN_PROVIDER [ double precision, psi_lapl, (elec_num_8) ]
   BEGIN_DOC
   ! Laplacian of the wave function
   END_DOC
-  
+  double precision :: psi_lapl2(elec_num)
+  double precision :: norm_Err
   integer          :: i, j
+  if (utilise_svd) then
+    !DIR$ VECTOR ALIGNED
+    !DIR$ LOOP COUNT (100)
+    do j=1,elec_num
+       psi_lapl(j) = jast_value*(psidet_grad_lapl_svd(4,j) + psidet_value_svd*jast_lapl_jast_inv(j) + 2.d0*(&
+           psidet_grad_lapl_svd(1,j)*jast_grad_jast_inv_x(j) +                   &
+           psidet_grad_lapl_svd(2,j)*jast_grad_jast_inv_y(j) +                   &
+           psidet_grad_lapl_svd(3,j)*jast_grad_jast_inv_z(j) ))
+    enddo
+  else
     !DIR$ VECTOR ALIGNED
     !DIR$ LOOP COUNT (100)
     do j=1,elec_num
@@ -42,7 +71,16 @@ BEGIN_PROVIDER [ double precision, psi_lapl, (elec_num_8) ]
           psidet_grad_lapl(2,j)*jast_grad_jast_inv_y(j) +                   &
           psidet_grad_lapl(3,j)*jast_grad_jast_inv_z(j) ))
     enddo
+  endif
 
+  !norm_Err = sum( (psi_lapl2(1:elec_num) - psi_lapl(1:elec_num))**2 ) / sum( psi_lapl(1:elec_num)**2 )
+  !if (norm_Err > 1.d-6) then
+  !   print *, 'probleme dans PROVIDER [ psi_lapl ]:  norm_Err = ', norm_Err
+     !print *, psi_lapl2(1:elec_num) 
+     !print *, psi_lapl(1:elec_num) 
+     !print *, irp_here
+     !stop
+  !endif 
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x, (elec_num_8) ]
@@ -53,7 +91,18 @@ END_PROVIDER
 ! grad(psi)/psi
   END_DOC
   
+  double precision :: psi_grad_psi_inv_x2(elec_num), psi_grad_psi_inv_y2(elec_num), psi_grad_psi_inv_z2(elec_num)
+  double precision :: norm_Err
   integer          :: j
+  if (utilise_svd) then
+    !DIR$ VECTOR ALIGNED
+    !DIR$ LOOP COUNT (100)
+    do j=1,elec_num
+       psi_grad_psi_inv_x(j) = psidet_grad_lapl_svd(1,j)*psidet_inv_svd + jast_grad_jast_inv_x(j)
+       psi_grad_psi_inv_y(j) = psidet_grad_lapl_svd(2,j)*psidet_inv_svd + jast_grad_jast_inv_y(j)
+       psi_grad_psi_inv_z(j) = psidet_grad_lapl_svd(3,j)*psidet_inv_svd + jast_grad_jast_inv_z(j)
+    enddo
+  else
     !DIR$ VECTOR ALIGNED
     !DIR$ LOOP COUNT (100)
     do j=1,elec_num
@@ -61,6 +110,16 @@ END_PROVIDER
       psi_grad_psi_inv_y(j) = psidet_grad_lapl(2,j)*psidet_inv + jast_grad_jast_inv_y(j)
       psi_grad_psi_inv_z(j) = psidet_grad_lapl(3,j)*psidet_inv + jast_grad_jast_inv_z(j)
     enddo
+  endif
+  !norm_Err = sum( (psi_grad_psi_inv_x2(1:elec_num) - psi_grad_psi_inv_x(1:elec_num))**2 ) / sum( psi_grad_psi_inv_x(1:elec_num)**2 )
+  !norm_Err = norm_Err + sum( (psi_grad_psi_inv_y2(1:elec_num) - psi_grad_psi_inv_y(1:elec_num))**2 ) / sum( psi_grad_psi_inv_y(1:elec_num)**2 )
+  !norm_Err = norm_Err + sum( (psi_grad_psi_inv_z2(1:elec_num) - psi_grad_psi_inv_z(1:elec_num))**2 ) / sum( psi_grad_psi_inv_z(1:elec_num)**2 )
+  !if (norm_Err > 1.d-6) then
+  !   print *, 'probleme dans PROVIDER [ psi_grad_psi_inv_xyz ]:  norm_Err = ', norm_Err
+     !print *, irp_here
+     !stop
+  !endif   
+
 
 END_PROVIDER
 

src/psi_SVD.irp.f

@@ -0,0 +1,85 @@
+ BEGIN_PROVIDER [ double precision, psi_value_SVD ]
+  implicit none
+  BEGIN_DOC
+  ! Value of the wave function
+  END_DOC
+  psi_value_SVD = psidet_value_SVD * jast_value
+  if (psi_value_SVD == 0.d0) then
+    call abrt(irp_here,"Value of the wave function is 0.")
+  endif
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, psi_value_inv_SVD ]
+  implicit none
+  BEGIN_DOC
+  ! 1. / psi_value_SVD
+  END_DOC
+  psi_value_inv_SVD = 1.d0 / psi_value_SVD
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, psi_value_inv2_SVD ]
+ implicit none
+ BEGIN_DOC
+ ! 1./(psi_value_SVD)**2
+ END_DOC
+ psi_value_inv2_SVD = psi_value_inv_SVD * psi_value_inv_SVD
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, psi_lapl_SVD, (elec_num_8) ]
+  implicit none
+  BEGIN_DOC
+  ! Laplacian of the wave function
+  END_DOC
+  integer :: i, j
+  !DIR$ VECTOR ALIGNED
+  !DIR$ LOOP COUNT (100)
+  do j = 1, elec_num
+    psi_lapl_SVD(j) = jast_value * (                           &
+        psidet_grad_lapl_SVD(4,j) +                            &
+        psidet_value_SVD * jast_lapl_jast_inv(j) + 2.d0 * (    &
+        psidet_grad_lapl_SVD(1,j) * jast_grad_jast_inv_x(j) +  &
+        psidet_grad_lapl_SVD(2,j) * jast_grad_jast_inv_y(j) +  &
+        psidet_grad_lapl_SVD(3,j) * jast_grad_jast_inv_z(j) ) )
+  enddo
+END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_x_SVD, (elec_num_8) ]
+&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_y_SVD, (elec_num_8) ]
+&BEGIN_PROVIDER [ double precision, psi_grad_psi_inv_z_SVD, (elec_num_8) ]
+  implicit none
+  BEGIN_DOC
+  ! grad(psi_SVD) / psi_SVD
+  END_DOC
+  integer :: j
+  !DIR$ VECTOR ALIGNED
+  !DIR$ LOOP COUNT (100)
+  do j = 1, elec_num
+    psi_grad_psi_inv_x_SVD(j) = psidet_grad_lapl_SVD(1,j) * psidet_inv_SVD + jast_grad_jast_inv_x(j)
+    psi_grad_psi_inv_y_SVD(j) = psidet_grad_lapl_SVD(2,j) * psidet_inv_SVD + jast_grad_jast_inv_y(j)
+    psi_grad_psi_inv_z_SVD(j) = psidet_grad_lapl_SVD(3,j) * psidet_inv_SVD + jast_grad_jast_inv_z(j)
+  enddo 
+END_PROVIDER
+
+
+
+BEGIN_PROVIDER [ double precision, psi_lapl_psi_inv_SVD, (elec_num_8) ]
+  implicit none
+  BEGIN_DOC
+  ! (Laplacian psi_SVD) / psi_SVD
+  END_DOC
+  integer :: i, j
+  !DIR$ VECTOR ALIGNED
+  !DIR$ LOOP COUNT (100)
+  do j = 1, elec_num
+    psi_lapl_psi_inv_SVD(j) = psi_lapl_SVD(j) * psi_value_inv_SVD
+  enddo
+END_PROVIDER
+

src/test_SVD/QMCCHEM_withSVD.py

@@ -0,0 +1,325 @@
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from math import sqrt
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+from scipy.linalg import eig, eigh
+
+# PARAMETERS
+block_time = 20
+eps = 1.
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_energy():
+    buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],  encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+      return None, None
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def run_qmc():
+    return subprocess.check_output(['qmcchem', 'run', filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def stop_qmc():
+    subprocess.check_output(['qmcchem', 'stop', filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def set_vmc_params():
+    #subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
+    #                                 '-m', 'VMC',
+    #                                 '-l', str(20),
+    #                                 '--time-step=0.3',
+    #                                 '--stop-time=36000',
+    #                                 '--norm=1.e-5',
+    #                                 '-w', '10',
+    #                                 filename])
+    subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'None', '-l', str(block_time), filename])
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_svd(n_svd):
+    # !!!
+    Ci_h_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_h_matrix_svd = results.find('Ci_h_matrix_svd : [ ') + len( 'Ci_h_matrix_svd : [ ' ) 
+    end_Ci_h_matrix_svd = len(results)#results.find('E_loc :') 
+    Ci_h_matrix_svd_buf = results[beg_Ci_h_matrix_svd:end_Ci_h_matrix_svd]
+    Ci_h_matrix_svd_buf = Ci_h_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_h_matrix_svd_buf[iline].split()
+        #print(line)
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_svd')
+            stop
+        else:
+          #Ci_h_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+        # !!!
+    # !!!
+    # Ci_h_matrix_svd = np.reshape(Ci_h_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_h_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_svd(n_svd):
+    # !!!
+    Ci_overlap_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_svd = results.find('Ci_overlap_matrix_svd : [ ') + len( 'Ci_overlap_matrix_svd : [ ' ) 
+    end_Ci_overlap_matrix_svd = len(results)#results.find('Ci_h_matrix_svd : [') 
+    Ci_overlap_matrix_svd_buf = results[beg_Ci_overlap_matrix_svd:end_Ci_overlap_matrix_svd]
+    Ci_overlap_matrix_svd_buf = Ci_overlap_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_overlap_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        # !!!
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_svd')
+            stop
+            # !!!
+        # !!!
+        else:
+          #Ci_overlap_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+          # !!!
+      # !!!
+    # !!!
+    #Ci_overlap_matrix_svd = np.reshape(Ci_overlap_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_postsvd(n_svd):
+    #file = open('verif_order.txt','w')
+    # !!!
+    Ci_h_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_h_matrix_postsvd = results.find('Ci_h_matrix_postsvd : [ ') + len( 'Ci_h_matrix_postsvd : [ ' ) 
+    end_Ci_h_matrix_postsvd = len(results)
+    Ci_h_matrix_postsvd_buf = results[beg_Ci_h_matrix_postsvd:end_Ci_h_matrix_postsvd]
+    Ci_h_matrix_postsvd_buf = Ci_h_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_h_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_postsvd')
+            stop
+        else:
+            # !!!
+            kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            #indcrep = kp + (k-1)*n_svd + (lp-1)*n_svd**2 + (l-1)*n_svd**3
+            #file.write( '{:5} {:5} {:5} {:5} {:5} {:5} \n'.format(indc, indc-indcrep, kp, k, lp, l ) )
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_h_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_h_matrix_postsvd[indc-1] = float( line[2] )
+        # !!!
+    #Ci_h_matrix_postsvd = np.reshape(Ci_h_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    #file.close()
+    return(Ci_h_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_postsvd(n_svd):
+    # !!!
+    Ci_overlap_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_postsvd = results.find('Ci_overlap_matrix_postsvd : [ ') + len( 'Ci_overlap_matrix_postsvd : [ ' ) 
+    end_Ci_overlap_matrix_postsvd = len(results)
+    Ci_overlap_matrix_postsvd_buf = results[beg_Ci_overlap_matrix_postsvd:end_Ci_overlap_matrix_postsvd]
+    Ci_overlap_matrix_postsvd_buf = Ci_overlap_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_overlap_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_postsvd')
+            stop
+        else:
+            # !!!
+            kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_overlap_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_overlap_matrix_postsvd[indc-1] = float( line[2] )
+            # !!!
+        # !!!
+    #Ci_overlap_matrix_postsvd = np.reshape(Ci_overlap_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def check_symmetric(a, tol=1e-3):
+    return np.all(np.abs(a-a.T) < tol)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    t0 = time.time()
+    # !!!
+    filename = "/home/aammar/qp2/src/svdwf/h2o_pseudo"
+    ezfio.set_file(filename)
+    # !!!
+    psi_svd_coeff = np.array(ezfio.get_spindeterminants_psi_svd_coefs())
+    U_svd = np.array(ezfio.get_spindeterminants_psi_svd_alpha())
+    V_svd = np.array(ezfio.get_spindeterminants_psi_svd_beta())
+    # !!!
+    print("Today's date:", datetime.now() )
+    print("filename = {}".format(filename))
+    # !!!
+    #set_vmc_params()
+    #run_qmc()
+    #print("start QMC:")
+    #stop_qmc()
+    # !!!
+    results =  subprocess.check_output(['qmcchem', 'result', filename],  encoding='UTF-8')
+    # !!!
+    #file = open('results.txt','a') 
+    #file.write('\n \n \n')
+    #file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    #file.write("Today's date: {}\n".format(datetime.now()) )
+    #file.write("filename = {}\n".format(filename))
+    #file.write('\n')
+    #file.write( results )
+    #file.write('\n')
+    #file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
+    #file.close()
+    # !!!
+    E_loc, ErrEloc = get_energy()
+    print('Eloc = {} +/- {}'.format(E_loc, ErrEloc))
+    # !!!
+    n_svd                 = 10
+    E_toadd               = 6.983610961797779
+    # !!!
+    Ci_h_matrix_svd       = get_Ci_h_matrix_svd(n_svd)
+    #print( 'Ci_h_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_svd)) )
+    Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd(n_svd)
+    #print( 'Ci_overlap_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_svd)) )
+    # !!!
+    aa = Ci_h_matrix_svd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_svd
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, check_finite=True, homogeneous_eigvals=False)
+    print( eigvals_svd + E_toadd )
+    # !!!
+    recouvre_svd = np.abs(psi_svd_coeff @ vr)
+    ind_gssvd = np.argmax(recouvre_svd)
+    print('svd : ', ind_gssvd, eigvals_svd[ind_gssvd] + E_toadd )
+    # !!!
+    Ci_h_matrix_postsvd       = get_Ci_h_matrix_postsvd(n_svd)
+    #print( 'Ci_h_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_postsvd)) )
+    Ci_overlap_matrix_postsvd = get_Ci_overlap_matrix_postsvd(n_svd)
+    #print( 'Ci_overlap_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_postsvd)) )
+    # !!!
+    aa = Ci_h_matrix_postsvd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_postsvd
+    eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, check_finite=True, homogeneous_eigvals=False)
+    print( eigvals_postsvd + E_toadd )
+    # !!!
+    d_postsvd = np.diagflat(psi_svd_coeff)
+    d_postsvd = d_postsvd.reshape( (1,n_svd*n_svd) )
+    recouvre_postsvd = np.abs(d_postsvd @ vr)
+    ind_gspostsvd = np.argmax(recouvre_postsvd)
+    print('post svd :', ind_gspostsvd, eigvals_postsvd[ind_gspostsvd]+E_toadd )
+    # !!!
+    print("end code after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!! 
+# !!!
+
+

src/test_SVD/h2o/QMCCHEM_withSVD.py

@@ -0,0 +1,481 @@
+# !!!
+import sys, os
+QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
+sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
+# !!!
+from ezfio import ezfio
+from math import sqrt
+from datetime import datetime
+import time
+import numpy as np
+import subprocess
+from scipy.linalg import eig, eigh
+from RSVD import powit_RSVD  
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_energy():
+    buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', EZFIO_file],  encoding='UTF-8')
+    if buffer.strip() != "":
+        buffer = buffer.splitlines()[-1]
+        _, energy, error = [float(x) for x in buffer.split()]
+        return energy, error
+    else:
+      return None, None
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_svd():
+    # !!!
+    Ci_h_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_h_matrix_svd = results.find('Ci_h_matrix_svd : [ ') + len( 'Ci_h_matrix_svd : [ ' ) 
+    end_Ci_h_matrix_svd = len(results)
+    Ci_h_matrix_svd_buf = results[beg_Ci_h_matrix_svd:end_Ci_h_matrix_svd]
+    Ci_h_matrix_svd_buf = Ci_h_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_h_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_svd')
+            break
+        else:
+          #Ci_h_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_h_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_h_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+        # !!!
+    # !!!
+    # Ci_h_matrix_svd = np.reshape(Ci_h_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_h_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_svd():
+    # !!!
+    Ci_overlap_matrix_svd = np.zeros( (n_svd,n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_svd = results.find('Ci_overlap_matrix_svd : [ ') + len( 'Ci_overlap_matrix_svd : [ ' ) 
+    end_Ci_overlap_matrix_svd = len(results) 
+    Ci_overlap_matrix_svd_buf = results[beg_Ci_overlap_matrix_svd:end_Ci_overlap_matrix_svd]
+    Ci_overlap_matrix_svd_buf = Ci_overlap_matrix_svd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**2+1):
+        # !!!
+        line = Ci_overlap_matrix_svd_buf[iline].split()
+        indc = int( line[0] )
+        # !!!
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_svd')
+            break
+            # !!!
+        # !!!
+        else:
+          #Ci_overlap_matrix_svd[indc-1] = float( line[2] )
+          irow = indc %  n_svd
+          icol = indc // n_svd
+          if( irow!=0 ):
+              Ci_overlap_matrix_svd[irow-1][icol] = float( line[2] )
+          else:
+              Ci_overlap_matrix_svd[n_svd-1][icol-1] = float( line[2] )
+          # !!!
+      # !!!
+    # !!!
+    #Ci_overlap_matrix_svd = np.reshape(Ci_overlap_matrix_svd, (n_svd, n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_svd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_h_matrix_postsvd():
+    #file = open('verif_order.txt','w')
+    # !!!
+    Ci_h_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_h_matrix_postsvd = results.find('Ci_h_matrix_postsvd : [ ') + len( 'Ci_h_matrix_postsvd : [ ' ) 
+    end_Ci_h_matrix_postsvd = len(results)
+    Ci_h_matrix_postsvd_buf = results[beg_Ci_h_matrix_postsvd:end_Ci_h_matrix_postsvd]
+    Ci_h_matrix_postsvd_buf = Ci_h_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_h_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_h_matrix_postsvd')
+            break
+        else:
+            # !!!
+            kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            #indcrep = kp + (k-1)*n_svd + (lp-1)*n_svd**2 + (l-1)*n_svd**3
+            #file.write( '{:5} {:5} {:5} {:5} {:5} {:5} \n'.format(indc, indc-indcrep, kp, k, lp, l ) )
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_h_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_h_matrix_postsvd[indc-1] = float( line[2] )
+        # !!!
+    #Ci_h_matrix_postsvd = np.reshape(Ci_h_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    #file.close()
+    return(Ci_h_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Ci_overlap_matrix_postsvd():
+    # !!!
+    Ci_overlap_matrix_postsvd = np.zeros( (n_svd*n_svd , n_svd*n_svd) )
+    # !!!
+    beg_Ci_overlap_matrix_postsvd = results.find('Ci_overlap_matrix_postsvd : [ ') + len( 'Ci_overlap_matrix_postsvd : [ ' ) 
+    end_Ci_overlap_matrix_postsvd = len(results)
+    Ci_overlap_matrix_postsvd_buf = results[beg_Ci_overlap_matrix_postsvd:end_Ci_overlap_matrix_postsvd]
+    Ci_overlap_matrix_postsvd_buf = Ci_overlap_matrix_postsvd_buf.split( '\n' )
+    # !!!
+    for iline in range(1, n_svd**4+1):
+        # !!!
+        line = Ci_overlap_matrix_postsvd_buf[iline].split()
+        indc = int( line[0] )
+        errS = float( line[4] ) 
+        #if( errS>eps ):
+            #print( line )
+        if( indc != iline ):
+            print('Error in reading Ci_overlap_matrix_postsvd')
+            break
+        else:
+            # !!!
+            #kp = indc % n_svd
+            if( ( indc % n_svd ) !=0 ):
+                kp = indc % n_svd
+            else:
+                kp = n_svd
+            indc1 = int( ( indc - kp ) / n_svd )
+            k  = indc1 % n_svd + 1
+            indc2 = int( ( indc1 - (k-1) ) / n_svd )
+            lp = indc2 % n_svd + 1
+            l  = int( ( indc2 - (lp-1) ) / n_svd ) + 1
+            # !!!
+            irow = kp + (k-1)*n_svd - 1
+            icol = lp + (l-1)*n_svd - 1
+            Ci_overlap_matrix_postsvd[irow][icol] = float( line[2] )
+            #Ci_overlap_matrix_postsvd[indc-1] = float( line[2] )
+            # !!!
+        # !!!
+    #Ci_overlap_matrix_postsvd = np.reshape(Ci_overlap_matrix_postsvd, (n_svd*n_svd, n_svd*n_svd), order='F')
+    # !!!
+    return(Ci_overlap_matrix_postsvd)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def check_symmetric(a, tol=1e-3):
+    return np.all(np.abs(a-a.T) < tol)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def save_results_to_resultsQMC():
+    file = open('resultsQMC.txt','a') 
+    file.write('\n \n \n')
+    file.write('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n \n')
+    file.write("Today's date: {}\n".format(datetime.now()))
+    file.write("EZFIO file = {}\n".format(EZFIO_file))
+    file.write('\n')
+    file.write( results )
+    file.write('\n')
+    file.write('+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +\n \n')
+    file.close()
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Esvd():
+    # !!!
+    # read CI_SVD matrices
+    Ci_h_matrix_svd       = get_Ci_h_matrix_svd()
+    Ci_overlap_matrix_svd = get_Ci_overlap_matrix_svd()
+    #print( 'Ci_h_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_svd)) )
+    #print( 'Ci_overlap_matrix_svd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_svd)) )
+    # !!!
+    # symmetrise and diagonalise
+    aa = Ci_h_matrix_svd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_svd
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False)
+    #print( eigvals_svd + E_toadd )
+    recouvre_svd = np.abs(psi_svd_coeff @ vr)
+    ind_gssvd = np.argmax(recouvre_svd)
+    # !!!
+    E_svd = eigvals_svd[ind_gssvd] + E_toadd
+    return( E_svd, vr[:,ind_gssvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def get_Epostsvd():
+    # !!!
+    # read CI_postSVD matrices
+    Ci_h_matrix_postsvd       = get_Ci_h_matrix_postsvd()
+    Ci_overlap_matrix_postsvd = get_Ci_overlap_matrix_postsvd()
+    #print( 'Ci_h_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_h_matrix_postsvd)) )
+    #print( 'Ci_overlap_matrix_postsvd is symmetric ? {}' .format(check_symmetric(Ci_overlap_matrix_postsvd)) )
+    # !!!
+    # symmetrise and diagonalise
+    aa = Ci_h_matrix_postsvd
+    aa = 0.5*( aa + aa.T )
+    bb = Ci_overlap_matrix_postsvd
+    eigvals_postsvd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False)
+    #print( eigvals_postsvd + E_toadd )
+    d_postsvd = np.diagflat(psi_svd_coeff)
+    d_postsvd = d_postsvd.reshape( (1,n_svd*n_svd) )
+    recouvre_postsvd = np.abs(d_postsvd @ vr)
+    ind_gspostsvd = np.argmax(recouvre_postsvd)
+    #print(recouvre_postsvd, ind_gspostsvd)
+    # !!!
+    #print( eigvals_postsvd.shape )
+    E_postsvd = eigvals_postsvd[ind_gspostsvd] + E_toadd
+    # !!!
+    return( E_postsvd,  vr[:,ind_gspostsvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+
+
+
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! 
+def SVD_postsvd(sigma_postsvd):
+    # !!!
+    print( 'performing new SVD for the post SVD eigenvector:' )
+    # !!!
+    #sigma_postsvd = sigma_postsvd.reshape( (n_svd,n_svd) )
+    #sigma_postsvd = sigma_postsvd.reshape( n_svd, n_svd, order='F' )
+    #print( 'sigma_postsvd is symmetric ? {}' .format(check_symmetric(sigma_postsvd)) )
+    # !!!
+    sigma_postsvd_mat = np.zeros( (n_svd,n_svd) )
+    for indc in range(1, n_svd**2+1):
+        if( ( indc % n_svd ) !=0 ):
+            kp = indc % n_svd
+        else:
+            kp = n_svd
+        indc1 = int( ( indc - kp ) / n_svd )
+        k     = indc1 % n_svd + 1
+        irow  = kp + (k-1)*n_svd - 1
+        sigma_postsvd_mat[kp-1][k-1] = sigma_postsvd[irow]
+    sigma_postsvd = sigma_postsvd_mat
+    # !!!
+    # construct the new matrix Y
+    Y     = U_svd @ sigma_postsvd @ V_svd.T
+    normY = np.linalg.norm(Y, ord='fro')
+    # !!!
+    # parameters of RSVD
+    rank        = n_svd
+    npow        = 10
+    nb_oversamp = 10
+    # !!!
+    # call RSV
+    U_postSVD, sigma_postsvd_diag, VT_postsvd = powit_RSVD(Y, rank, npow, nb_oversamp)
+    # !!!
+    # check precision
+    Y_SVD = np.dot( U_postSVD , np.dot( np.diag(sigma_postsvd_diag) , VT_postsvd ) )
+    energy = np.sum( np.square(sigma_postsvd_diag) ) / normY**2
+    err_SVD = 100. * np.linalg.norm( Y - Y_SVD, ord="fro") / normY
+    print('energy = {}, error = {}\n'.format(energy, err_SVD))
+    # !!!
+    return(U_postSVD, sigma_postsvd_diag, VT_postsvd)
+    # !!!
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_Hsvd_QP(Hsvd_qp_txt):
+    Hsvd_qp      = np.zeros( (n_svd,n_svd) )
+    Hsvd_qp_file = open(Hsvd_qp_txt, 'r')
+    for line in Hsvd_qp_file:
+        line         = line.split()
+        i            = int(line[0]) - 1
+        j            = int(line[1]) - 1
+        Hsvd_qp[i,j] = float(line[2])
+    return(Hsvd_qp)
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+def get_Esvd_QP(Hsvd_qp):
+    # !!!
+    # symmetrise and diagonalise
+    aa = Hsvd_qp
+    aa = 0.5*( aa + aa.T )
+    bb = np.identity(n_svd)
+    eigvals_svd, vr = eig(aa, bb, left=False, right=True, overwrite_a=True, overwrite_b=True, 
+            check_finite=True, homogeneous_eigvals=False) 
+    recouvre_svd = np.abs(psi_svd_coeff @ vr)
+    ind_gssvd = np.argmax(recouvre_svd)
+    E_svd = eigvals_svd[ind_gssvd] + E_toadd
+    return( E_svd, vr[:,ind_gssvd] )
+# ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ ! ~ !
+
+
+
+
+if __name__ == '__main__':
+    t0 = time.time()
+    # !!!
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    EZFIO_file = "/home/aammar/qp2/src/svdwf/h2o_631g_frez_Jopt_nsvd20"
+    E_toadd    =  9.194966082434476 #6.983610961797779
+    # ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ #
+    # !!!
+    ezfio.set_file(EZFIO_file)
+    n_svd         = ezfio.get_spindeterminants_n_svd_coefs()
+    psi_svd_coeff = np.array(ezfio.get_spindeterminants_psi_svd_coefs())
+    U_svd         = np.array(ezfio.get_spindeterminants_psi_svd_alpha())
+    V_svd         = np.array(ezfio.get_spindeterminants_psi_svd_beta())
+    # !!!
+    #print( U_svd.shape )
+    U_svd = U_svd[0,:,:].T
+    #print( U_svd.shape )
+    V_svd = V_svd[0,:,:].T
+    # !!!
+    print("Today's date:", datetime.now() )
+    print("EZFIO file = {}".format(EZFIO_file))
+    print("nuclear energy = {}".format(E_toadd) )
+    print("n_svd = {} \n".format(n_svd) )
+    # !!!
+    #set_vmc_params()
+    #run_qmc()
+    #print("start QMC:")
+    #stop_qmc()
+    # !!!
+    print( 'getting QMCCHEM results from {}'.format(EZFIO_file) )
+    results =  subprocess.check_output(['qmcchem', 'result', EZFIO_file],  encoding='UTF-8')
+    # !!!
+    Eloc, ErrEloc = get_energy()
+    print('Eloc = {} +/- {}'.format(Eloc, ErrEloc))
+    print('{} <= Eloc <= {}\n'.format(Eloc-ErrEloc,Eloc+ErrEloc))
+    # !!!
+    save_resultsQMC = input( 'save QMC results from {}? (y/n) '.format(EZFIO_file) )
+    if( save_resultsQMC == 'y' ):
+        print('saving in resultsQMC.txt')
+        save_results_to_resultsQMC()
+    print('\n')
+    # !!!
+    read_QPsvd = input( 'read QP Hsvd matrix ? (y/n) ')
+    if( read_QPsvd == 'y' ):
+        Hsvd_qp_txt  = input('name of file with QM H_svd matrix:')
+        Hsvd_qp      = get_Hsvd_QP(Hsvd_qp_txt)
+        E_svd_QP, sigma_svd_QP = get_Esvd_QP(Hsvd_qp)
+        print('QP SVD enegry = {} \n'.format(E_svd_QP) )
+        print('\n')
+    # !!!
+    E_svd, sigma_svd = get_Esvd()
+    print('svd enegry = {} \n '.format(E_svd) )
+    # !!!
+    E_postsvd, sigma_postsvd = get_Epostsvd()
+    print('post svd energy = {} \n'.format(E_postsvd) )
+    # !!!
+    save_to_EZFIO = input( "modify EZFIO (with svd_QP, svd_QMC or postsvd)? "  )
+    # !!!
+    if( save_to_EZFIO == 'svd_QMC' ):
+        ezfio.set_spindeterminants_psi_svd_coefs( sigma_svd )
+        direc_svdcoeff = EZFIO_file + '/spindeterminants/psi_svd_coefs.gz'
+        print( '{} is modified'.format(direc_svdcoeff) )
+        print( 'precedent svd coeff:' )
+        print(psi_svd_coeff)
+        print( 'current svd coeff:' )
+        print(sigma_svd)
+    # !!!
+    elif( save_to_EZFIO == 'svd_QP'):
+        ezfio.set_spindeterminants_psi_svd_coefs( sigma_svd_QP )
+        direc_svdcoeff = EZFIO_file + '/spindeterminants/psi_svd_coefs.gz'
+        print( '{} is modified'.format(direc_svdcoeff) )
+        print( 'precedent svd coeff:' )
+        print(psi_svd_coeff)
+        print( 'current svd coeff:' )
+        print(sigma_svd_QP)
+    # !!!
+    elif( save_to_EZFIO == 'postsvd' ):
+        # SVD first
+        U_postSVD, sigma_postsvd_diag, V_postSVD = SVD_postsvd(sigma_postsvd)
+        V_postSVD = V_postSVD.T
+        # save next in ezfio file   
+        U_postSVD_toEZFIO        = np.zeros( ( U_postSVD.shape[0], U_postSVD.shape[1], 1) )
+        V_postSVD_toEZFIO        = np.zeros( ( V_postSVD.shape[0], V_postSVD.shape[1], 1) )
+        U_postSVD_toEZFIO[:,:,0] = U_postSVD  
+        V_postSVD_toEZFIO[:,:,0] = V_postSVD  
+        #
+        ezfio.set_spindeterminants_psi_svd_alpha( U_postSVD_toEZFIO )        
+        ezfio.set_spindeterminants_psi_svd_coefs( sigma_postsvd_diag )
+        ezfio.set_spindeterminants_psi_svd_beta( V_postSVD_toEZFIO )
+        print(sigma_postsvd_diag)
+    # !!!
+    else:
+        print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+        exit()
+    # !!!
+    print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
+    # !!!
+# !!!
+
+

src/test_SVD/h2o/QR.py

@@ -0,0 +1,10 @@
+# !!!
+import numpy as np
+# !!!
+def QR_fact(X):
+    Q, R     = np.linalg.qr(X, mode="reduced")
+    D        = np.diag( np.sign( np.diag(R) ) )
+    Qunique  = np.dot(Q,D)
+    #Runique = np.dot(D,R)
+    return(Qunique)
+# !!!

src/test_SVD/h2o/RSVD.py

@@ -0,0 +1,20 @@
+# !!!
+import numpy as np
+from QR import QR_fact
+# !!!
+def powit_RSVD(X, new_r, nb_powit, nb_oversamp):
+    # !!!
+    G  = np.random.randn(X.shape[1], new_r+nb_oversamp)
+    Q  = QR_fact( np.dot(X,G) )
+    # !!!
+    for _ in range(nb_powit):
+        Q = QR_fact( np.dot(X.T,Q) )
+        Q = QR_fact( np.dot(X,Q) )
+        # !!!
+    Y = np.dot(Q.T,X)
+    # !!!
+    U, S, VT = np.linalg.svd(Y, full_matrices=0)
+    U = np.dot(Q,U)
+    return U[:,:(new_r)], S[:(new_r)], VT[:(new_r),:]
+    # !!!
+# !!!

src/wf.irp.f

@@ -112,7 +112,11 @@ END_PROVIDER
 
   allocate (psi_det_tmp (N_int,max(det_alpha_num,det_beta_num)))
 
+  if (utilise_svd) then
+    t = -1.d0 
+  else
     t = ci_threshold
+  endif
 
   ! Compute the norm of the alpha and beta determinants
   d_alpha = 0.d0