mirror of
https://gitlab.com/scemama/qmcchem.git
synced 2024-08-29 22:33:41 +02:00
cluster: Jast opt for SCF WV
This commit is contained in:
commit
a93803151e
@ -8,10 +8,13 @@ import os
|
|||||||
import time
|
import time
|
||||||
import subprocess
|
import subprocess
|
||||||
from math import sqrt
|
from math import sqrt
|
||||||
|
#
|
||||||
|
#from bayes_opt import BayesianOptimization
|
||||||
|
|
||||||
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
|
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
|
||||||
|
# this should give: /home/ammar/qmcchem
|
||||||
sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
|
sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
|
||||||
|
# this should add: /home/ammar/qmcchem/EZFIO/Python/
|
||||||
|
|
||||||
from ezfio import ezfio
|
from ezfio import ezfio
|
||||||
|
|
||||||
@ -22,7 +25,7 @@ block_time = 20
|
|||||||
def main():
|
def main():
|
||||||
if len(sys.argv) != 2:
|
if len(sys.argv) != 2:
|
||||||
print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
|
print("Usage: %s <EZFIO_DIRECTORY>"%sys.argv[0])
|
||||||
sys.exti(1)
|
sys.exit(1)
|
||||||
|
|
||||||
filename = sys.argv[1]
|
filename = sys.argv[1]
|
||||||
ezfio.set_file(filename)
|
ezfio.set_file(filename)
|
||||||
@ -47,7 +50,7 @@ def main():
|
|||||||
|
|
||||||
def get_params_pen():
|
def get_params_pen():
|
||||||
d = ezfio.jastrow_jast_pen
|
d = ezfio.jastrow_jast_pen
|
||||||
return np.array([d[m[0]] for m in atom_map])
|
return np.array( [d[m[0]] for m in atom_map])
|
||||||
|
|
||||||
|
|
||||||
def get_energy():
|
def get_energy():
|
||||||
@ -60,6 +63,21 @@ def main():
|
|||||||
else:
|
else:
|
||||||
return None, None
|
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():
|
def get_variance():
|
||||||
buffer = subprocess.check_output(['qmcchem', 'result', '-e',
|
buffer = subprocess.check_output(['qmcchem', 'result', '-e',
|
||||||
'e_loc_qmcvar', filename],
|
'e_loc_qmcvar', filename],
|
||||||
@ -99,28 +117,53 @@ def main():
|
|||||||
# '-w', '10',
|
# '-w', '10',
|
||||||
# filename])
|
# filename])
|
||||||
subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple',
|
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.}
|
memo_energy = {'fmin': 100000000.}
|
||||||
def f(x):
|
def f(x):
|
||||||
print ("x = %s"%str(x))
|
print ("x = %s"%str(x))
|
||||||
|
# !!!
|
||||||
h = str(x)
|
h = str(x)
|
||||||
if h in memo_energy:
|
if h in memo_energy:
|
||||||
return memo_energy[h]
|
return memo_energy[h]
|
||||||
|
# !!!
|
||||||
set_params_pen(x)
|
set_params_pen(x)
|
||||||
set_vmc_params()
|
set_vmc_params()
|
||||||
|
# !!!
|
||||||
pid = os.fork()
|
pid = os.fork()
|
||||||
if pid == 0:
|
if pid == 0:
|
||||||
|
# In child process
|
||||||
run_qmc()
|
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)
|
os._exit(os.EX_OK)
|
||||||
else:
|
else:
|
||||||
|
# In parent process
|
||||||
import atexit
|
import atexit
|
||||||
atexit.register(stop_qmc)
|
atexit.register(stop_qmc)
|
||||||
|
# !!!
|
||||||
err = thresh+1.
|
|
||||||
time.sleep(3.*block_time/4.)
|
time.sleep(3.*block_time/4.)
|
||||||
|
# !!!
|
||||||
|
err = thresh+1.
|
||||||
local_thresh = thresh
|
local_thresh = thresh
|
||||||
while err > local_thresh:
|
while err > local_thresh:
|
||||||
time.sleep(block_time)
|
time.sleep(block_time)
|
||||||
@ -128,14 +171,20 @@ def main():
|
|||||||
variance, v_err = get_variance()
|
variance, v_err = get_variance()
|
||||||
if e is None or variance is None:
|
if e is None or variance is None:
|
||||||
continue
|
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:
|
if (energy-2.*err) > memo_energy['fmin']+thresh:
|
||||||
local_thresh = 10.*thresh
|
local_thresh = 10.*thresh
|
||||||
elif (energy+2.*err) < memo_energy['fmin']-thresh:
|
elif (energy+2.*err) < memo_energy['fmin']-thresh:
|
||||||
local_thresh = 10.*thresh
|
local_thresh = 10.*thresh
|
||||||
|
# !!!
|
||||||
# Check if PID is still running
|
# Check if PID is still running
|
||||||
try:
|
try:
|
||||||
os.kill(pid,0)
|
os.kill(pid,0)
|
||||||
@ -144,22 +193,34 @@ def main():
|
|||||||
break
|
break
|
||||||
stop_qmc()
|
stop_qmc()
|
||||||
os.wait()
|
os.wait()
|
||||||
memo_energy[h] = energy + err
|
# !!!
|
||||||
|
#memo_energy[h] = energy + err
|
||||||
|
memo_energy[h] = energy
|
||||||
|
# !!!
|
||||||
memo_energy['fmin'] = min(energy, memo_energy['fmin'])
|
memo_energy['fmin'] = min(energy, memo_energy['fmin'])
|
||||||
return energy
|
return energy
|
||||||
|
# !!!
|
||||||
|
|
||||||
def run():
|
def run():
|
||||||
x = get_params_pen()
|
x = get_params_pen()
|
||||||
|
# !!!
|
||||||
if sum(x) == 0.:
|
if sum(x) == 0.:
|
||||||
jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
|
jast_a_up_dn = ezfio.jastrow_jast_a_up_dn
|
||||||
x += 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))
|
print("x = "+str(opt))
|
||||||
set_params_pen(opt['x'])
|
set_params_pen(opt['x'])
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
run()
|
run()
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
if __name__ == '__main__':
|
if __name__ == '__main__':
|
||||||
main()
|
main()
|
||||||
|
@ -4,7 +4,7 @@
|
|||||||
BEGIN_PROVIDER [ double precision , jast_elec_Simple_value, (elec_num_8) ]
|
BEGIN_PROVIDER [ double precision , jast_elec_Simple_value, (elec_num_8) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
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
|
END_DOC
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
double precision :: a, b, rij, tmp
|
double precision :: a, b, rij, tmp
|
||||||
@ -20,10 +20,11 @@ implicit none
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
a = 0.5*jast_a_up_up
|
a = 0.5d0*jast_a_up_up
|
||||||
b = jast_b_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)
|
!DIR$ LOOP COUNT (50)
|
||||||
do i=j+1,elec_alpha_num
|
do i=j+1,elec_alpha_num
|
||||||
rij = elec_dist(i,j)
|
rij = elec_dist(i,j)
|
||||||
@ -33,7 +34,8 @@ implicit none
|
|||||||
enddo
|
enddo
|
||||||
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)
|
!DIR$ LOOP COUNT (50)
|
||||||
do i=j+1,elec_num
|
do i=j+1,elec_num
|
||||||
rij = elec_dist(i,j)
|
rij = elec_dist(i,j)
|
||||||
@ -43,9 +45,11 @@ implicit none
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
a = 0.5*jast_a_up_dn
|
|
||||||
|
a = 0.5d0*jast_a_up_dn
|
||||||
b = jast_b_up_dn
|
b = jast_b_up_dn
|
||||||
|
|
||||||
|
! anti-parallele spin
|
||||||
do j=1,elec_alpha_num
|
do j=1,elec_alpha_num
|
||||||
!DIR$ LOOP COUNT (50)
|
!DIR$ LOOP COUNT (50)
|
||||||
do i=elec_alpha_num+1,elec_num
|
do i=elec_alpha_num+1,elec_num
|
||||||
@ -190,6 +194,31 @@ BEGIN_PROVIDER [ double precision , jast_elec_Simple_lapl, (elec_num_8) ]
|
|||||||
enddo
|
enddo
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -13,6 +13,9 @@ program qmcchem_info
|
|||||||
endif
|
endif
|
||||||
print *, 'Number of determinants : ', det_num
|
print *, 'Number of determinants : ', det_num
|
||||||
print *, 'Number of unique alpha/beta determinants : ', det_alpha_num, det_beta_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 *, 'Closed-shell MOs : ', mo_closed_num
|
||||||
print *, 'Number of MOs in determinants : ', num_present_mos
|
print *, 'Number of MOs in determinants : ', num_present_mos
|
||||||
! print *, 'Det alpha norm:'
|
! print *, 'Det alpha norm:'
|
||||||
|
51
src/MAIN/test_debug.irp.f
Normal file
51
src/MAIN/test_debug.irp.f
Normal file
@ -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
|
@ -38,7 +38,7 @@ BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ]
|
|||||||
do l=1,elec_alpha_num
|
do l=1,elec_alpha_num
|
||||||
T += det_alpha_grad_lapl(4,l,i)*det_beta_value (j)
|
T += det_alpha_grad_lapl(4,l,i)*det_beta_value (j)
|
||||||
enddo
|
enddo
|
||||||
do l=elec_beta_num+1,elec_num
|
do l=elec_alpha_num+1,elec_num
|
||||||
T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i)
|
T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i)
|
||||||
enddo
|
enddo
|
||||||
ci_h_psidet(k) = -0.5d0*T + E_pot * det_alpha_value(i)*det_beta_value (j)
|
ci_h_psidet(k) = -0.5d0*T + E_pot * det_alpha_value(i)*det_beta_value (j)
|
||||||
@ -98,7 +98,7 @@ BEGIN_PROVIDER [ double precision, ci_h_matrix, (size_ci_h_matrix) ]
|
|||||||
do e=1,elec_alpha_num
|
do e=1,elec_alpha_num
|
||||||
g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
|
g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
|
||||||
enddo
|
enddo
|
||||||
do e=elec_beta_num+1,elec_num
|
do e=elec_alpha_num+1,elec_num
|
||||||
g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
|
g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
|
||||||
enddo
|
enddo
|
||||||
T = g
|
T = g
|
||||||
@ -169,7 +169,7 @@ BEGIN_PROVIDER [ double precision, ci_h_matrix_diag, (size_ci_h_matrix_diag) ]
|
|||||||
do e=1,elec_alpha_num
|
do e=1,elec_alpha_num
|
||||||
g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
|
g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
|
||||||
enddo
|
enddo
|
||||||
do e=elec_beta_num+1,elec_num
|
do e=elec_alpha_num+1,elec_num
|
||||||
g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
|
g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
|
||||||
enddo
|
enddo
|
||||||
T = g
|
T = g
|
||||||
|
@ -9,15 +9,17 @@ BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc1, (size_E_deriv_nucPar_lo
|
|||||||
|
|
||||||
integer :: i
|
integer :: i
|
||||||
|
|
||||||
do i=1,nucl_num
|
do i = 1, nucl_num
|
||||||
E_deriv_nucPar_loc1(i) = E_loc*jast_elec_Simple_deriv_nucPar(i)
|
E_deriv_nucPar_loc1(i) = E_loc*jast_elec_Simple_deriv_nucPar(i)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
E_deriv_nucPar_loc1_min = min(E_deriv_nucPar_loc1_min,minval(E_deriv_nucPar_loc1))
|
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))
|
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
|
SOFT_TOUCH E_deriv_nucPar_loc1_min E_deriv_nucPar_loc1_max
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc2, (size_E_deriv_nucPar_loc2) ]
|
BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc2, (size_E_deriv_nucPar_loc2) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
@ -37,6 +39,38 @@ BEGIN_PROVIDER [ double precision, E_deriv_nucPar_loc2, (size_E_deriv_nucPar_lo
|
|||||||
SOFT_TOUCH E_deriv_nucPar_loc2_min E_deriv_nucPar_loc2_max
|
SOFT_TOUCH E_deriv_nucPar_loc2_min E_deriv_nucPar_loc2_max
|
||||||
END_PROVIDER
|
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) ]
|
BEGIN_PROVIDER [ double precision, J_deriv_nucPar_verif, (size_J_deriv_nucPar_verif) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
@ -45,80 +79,157 @@ BEGIN_PROVIDER [ double precision, J_deriv_nucPar_verif, (size_J_deriv_nucPar_v
|
|||||||
! Dimensions : nucl_num
|
! Dimensions : nucl_num
|
||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
integer :: i,j
|
integer :: i, j
|
||||||
double precision :: sqt_eps = 1d-3, der1, der2, pos1, pos2, pos0
|
double precision :: eps = 1d-7, der1, der2, pos0
|
||||||
|
|
||||||
do j=1,nucl_num
|
do j = 1, nucl_num
|
||||||
!!
|
!!!
|
||||||
pos0 = sqt_eps * jast_pen(j)
|
pos0 = eps * jast_pen(j)
|
||||||
!!
|
!!!
|
||||||
jast_pen(j) = jast_pen(j) + pos0
|
jast_pen(j) = jast_pen(j) + pos0
|
||||||
TOUCH jast_pen
|
TOUCH jast_pen
|
||||||
pos1 = jast_pen(j)
|
der1 = 0.d0
|
||||||
!DIR$ LOOP COUNT (100)
|
!DIR$ LOOP COUNT (100)
|
||||||
do i=1,elec_num
|
do i = 1, elec_num
|
||||||
der1 += jast_elec_Simple_value(i)
|
der1 += jast_elec_Simple_value(i)
|
||||||
enddo
|
end do
|
||||||
!!
|
!!!
|
||||||
jast_pen(j) = jast_pen(j) - 2.d0 * pos0
|
jast_pen(j) = jast_pen(j) - 2.d0 * pos0
|
||||||
TOUCH jast_pen
|
TOUCH jast_pen
|
||||||
pos2 = jast_pen(j)
|
der2 = 0.d0
|
||||||
!DIR$ LOOP COUNT (100)
|
!DIR$ LOOP COUNT (100)
|
||||||
do i=1,elec_num
|
do i = 1, elec_num
|
||||||
der2 += jast_elec_Simple_value(i)
|
der2 += jast_elec_Simple_value(i)
|
||||||
enddo
|
end do
|
||||||
!!
|
!!!
|
||||||
J_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0
|
J_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0
|
||||||
!!
|
!!!
|
||||||
enddo
|
end do
|
||||||
|
|
||||||
J_deriv_nucPar_verif_min = min(J_deriv_nucPar_verif_min,minval(J_deriv_nucPar_verif))
|
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))
|
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
|
SOFT_TOUCH J_deriv_nucPar_verif_min J_deriv_nucPar_verif_max
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, E_deriv_nucPar_verif, (size_E_deriv_nucPar_verif) ]
|
BEGIN_PROVIDER [ double precision, J_deriv_nucPar_ex, (size_J_deriv_nucPar_ex) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Local energy variation with respect to nuclear parameters: verification
|
! Jastrow variation with respect to nuclear parameters
|
||||||
!
|
!
|
||||||
! Dimensions : nucl_num
|
! Dimensions : nucl_num
|
||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
integer :: j
|
integer :: j
|
||||||
double precision :: sqt_eps = 1d-3, der1, der2, pos1, pos2, pos0
|
|
||||||
|
|
||||||
do j=1,nucl_num
|
do j = 1, nucl_num
|
||||||
!!
|
J_deriv_nucPar_ex(j) = jast_elec_Simple_deriv_nucPar(j)
|
||||||
pos0 = sqt_eps * jast_pen(j)
|
end do
|
||||||
!!
|
|
||||||
jast_pen(j) = jast_pen(j) + pos0
|
|
||||||
TOUCH jast_pen
|
|
||||||
pos1 = jast_pen(j)
|
|
||||||
der1 = E_loc
|
|
||||||
!!DIR$ LOOP COUNT (100)
|
|
||||||
!do i=1,elec_num
|
|
||||||
! der1 += jast_elec_Simple_value(i)
|
|
||||||
!enddo
|
|
||||||
!!
|
|
||||||
jast_pen(j) = jast_pen(j) - 2.d0 * pos0
|
|
||||||
TOUCH jast_pen
|
|
||||||
pos2 = jast_pen(j)
|
|
||||||
der2 = E_loc
|
|
||||||
!!DIR$ LOOP COUNT (100)
|
|
||||||
!do i=1,elec_num
|
|
||||||
! der2 += jast_elec_Simple_value(i)
|
|
||||||
!enddo
|
|
||||||
!!
|
|
||||||
E_deriv_nucPar_verif(j) = 0.5d0 * (der1 - der2) / pos0
|
|
||||||
!!
|
|
||||||
enddo
|
|
||||||
|
|
||||||
E_deriv_nucPar_verif_min = min(E_deriv_nucPar_verif_min,minval(E_deriv_nucPar_verif))
|
J_deriv_nucPar_ex_min = min(J_deriv_nucPar_ex_min,minval(J_deriv_nucPar_ex))
|
||||||
E_deriv_nucPar_verif_max = max(E_deriv_nucPar_verif_max,maxval(E_deriv_nucPar_verif))
|
J_deriv_nucPar_ex_max = max(J_deriv_nucPar_ex_max,maxval(J_deriv_nucPar_ex))
|
||||||
SOFT_TOUCH E_deriv_nucPar_verif_min E_deriv_nucPar_verif_max
|
SOFT_TOUCH J_deriv_nucPar_ex_min J_deriv_nucPar_ex_max
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
||||||
|
@ -279,5 +279,3 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
|
|||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
262
src/QMC_SVD/v0/test_SVD.py
Executable file
262
src/QMC_SVD/v0/test_SVD.py
Executable file
@ -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.) )
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
|
|
10
src/QMC_SVD/v1/QR.py
Normal file
10
src/QMC_SVD/v1/QR.py
Normal file
@ -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)
|
||||||
|
# !!!
|
20
src/QMC_SVD/v1/RSVD.py
Normal file
20
src/QMC_SVD/v1/RSVD.py
Normal file
@ -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),:]
|
||||||
|
# !!!
|
||||||
|
# !!!
|
513
src/QMC_SVD/v1/diag_after_QMCCHEM.py
Executable file
513
src/QMC_SVD/v1/diag_after_QMCCHEM.py
Executable file
@ -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.) )
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
|
|
75
src/QMC_SVD/v1/save_svd.py
Normal file
75
src/QMC_SVD/v1/save_svd.py
Normal file
@ -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.) )
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
|
|
@ -1188,6 +1188,7 @@ END_PROVIDER
|
|||||||
|
|
||||||
ddet = 0.d0
|
ddet = 0.d0
|
||||||
|
|
||||||
|
! n_to_do < 2*elec_alpha_num
|
||||||
if (n_to_do < shiftl(elec_alpha_num,1)) then
|
if (n_to_do < shiftl(elec_alpha_num,1)) then
|
||||||
|
|
||||||
do while ( n_to_do > 0 )
|
do while ( n_to_do > 0 )
|
||||||
@ -1251,6 +1252,7 @@ END_PROVIDER
|
|||||||
ASSERT (ddet /= 0.d0)
|
ASSERT (ddet /= 0.d0)
|
||||||
|
|
||||||
det_alpha_value_curr = ddet
|
det_alpha_value_curr = ddet
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
|
BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
|
||||||
@ -1509,8 +1511,8 @@ END_PROVIDER
|
|||||||
DaC = 0.d0
|
DaC = 0.d0
|
||||||
CDb = 0.d0
|
CDb = 0.d0
|
||||||
|
|
||||||
if (det_num < shiftr(det_alpha_num*det_beta_num,2)) then
|
!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))
|
det_num4 = iand(det_num,not(3))
|
||||||
!DIR$ VECTOR ALIGNED
|
!DIR$ VECTOR ALIGNED
|
||||||
do k=1,det_num4,4
|
do k=1,det_num4,4
|
||||||
@ -1528,7 +1530,6 @@ END_PROVIDER
|
|||||||
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
|
CDb(i2) = CDb(i2) + det_coef_matrix_values(k+1)*f
|
||||||
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
|
CDb(i3) = CDb(i3) + det_coef_matrix_values(k+2)*f
|
||||||
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*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
|
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) &
|
DaC(j1) = DaC(j1) + det_coef_matrix_values(k)*det_alpha_value(i1) &
|
||||||
+ det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
|
+ det_coef_matrix_values(k+1)*det_alpha_value(i1+1) &
|
||||||
@ -1551,43 +1552,61 @@ END_PROVIDER
|
|||||||
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
|
CDb(i4) = CDb(i4) + det_coef_matrix_values(k+3)*det_beta_value (j4)
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
do k=det_num4+1,det_num
|
do k=det_num4+1,det_num
|
||||||
i = det_coef_matrix_rows(k)
|
i = det_coef_matrix_rows(k)
|
||||||
j = det_coef_matrix_columns(k)
|
j = det_coef_matrix_columns(k)
|
||||||
DaC(j) = DaC(j) + det_coef_matrix_values(k)*det_alpha_value(i)
|
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)
|
CDb(i) = CDb(i) + det_coef_matrix_values(k)*det_beta_value (j)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
! DaC(det_beta_num) = Trans( det_coef_matrix_dense(det_alpha_num,det_beta_num) ) @ det_alpha_value(det_alpha_num)
|
! 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, &
|
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)
|
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)
|
! 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, &
|
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)
|
size(det_coef_matrix_dense,1), det_beta_value, 1, 0.d0, CDb, 1)
|
||||||
|
|
||||||
endif
|
endif
|
||||||
|
!endif
|
||||||
|
|
||||||
|
|
||||||
! Value
|
! Value
|
||||||
! -----
|
! -----
|
||||||
|
|
||||||
psidet_value = 0.d0
|
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
|
do j=1,det_beta_num
|
||||||
psidet_value = psidet_value + det_beta_value(j) * DaC(j)
|
psidet_value = psidet_value + det_beta_value(j) * DaC(j)
|
||||||
enddo
|
enddo
|
||||||
|
endif
|
||||||
|
|
||||||
if (psidet_value == 0.d0) then
|
if (psidet_value == 0.d0) then
|
||||||
call abrt(irp_here,'Determinantal component of the wave function is zero.')
|
call abrt(irp_here,'Determinantal component of the wave function is zero.')
|
||||||
endif
|
endif
|
||||||
psidet_inv = 1.d0/psidet_value
|
psidet_inv = 1.d0/psidet_value
|
||||||
|
|
||||||
|
|
||||||
! Gradients
|
! Gradients
|
||||||
! ---------
|
! ---------
|
||||||
|
if(det_num .eq. 1) then
|
||||||
! psidet_grad_lapl(4,elec_alpha_num) = det_alpha_grad_lapl(4,elec_alpha_num,det_alpha_num) @ CDb(det_alpha_num)
|
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, &
|
call dgemv('N',elec_alpha_num*4,det_alpha_num,1.d0, &
|
||||||
det_alpha_grad_lapl, &
|
det_alpha_grad_lapl, &
|
||||||
size(det_alpha_grad_lapl,1)*size(det_alpha_grad_lapl,2), &
|
size(det_alpha_grad_lapl,1)*size(det_alpha_grad_lapl,2), &
|
||||||
@ -1598,6 +1617,8 @@ END_PROVIDER
|
|||||||
size(det_beta_grad_lapl,1)*size(det_beta_grad_lapl,2), &
|
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)
|
DaC, 1, 0.d0, psidet_grad_lapl(1,elec_alpha_num+1), 1)
|
||||||
endif
|
endif
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
if (do_pseudo) then
|
if (do_pseudo) then
|
||||||
call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv, &
|
call dgemv('N',elec_alpha_num,det_alpha_num,psidet_inv, &
|
||||||
@ -1750,7 +1771,8 @@ END_PROVIDER
|
|||||||
! !!!
|
! !!!
|
||||||
psidet_value_SVD = 0.d0
|
psidet_value_SVD = 0.d0
|
||||||
do l = 1, n_svd_coefs
|
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)
|
psidet_value_SVD = psidet_value_SVD + det_beta_value_SVD(l) &
|
||||||
|
* psi_svd_coefs(l,1) * det_alpha_value_SVD(l)
|
||||||
enddo
|
enddo
|
||||||
if (psidet_value_SVD == 0.d0) then
|
if (psidet_value_SVD == 0.d0) then
|
||||||
call abrt(irp_here,'Determinantal component of the SVD wave function is zero.')
|
call abrt(irp_here,'Determinantal component of the SVD wave function is zero.')
|
||||||
@ -1763,14 +1785,14 @@ END_PROVIDER
|
|||||||
! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
|
! -~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~-~
|
||||||
! !!!
|
! !!!
|
||||||
if (do_pseudo) then
|
if (do_pseudo) then
|
||||||
! det_alpha_pseudo_SVD = det_alpha_pseudo @ psi_svd_alpha / psidet_inv_SVD
|
! det_alpha_pseudo_SVD = det_alpha_pseudo @ psi_svd_alpha * psidet_inv_SVD
|
||||||
call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, 1.d0 &
|
call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, psidet_inv_SVD &
|
||||||
, det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) &
|
, det_alpha_pseudo, size(det_alpha_pseudo,1), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) &
|
||||||
, 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) )
|
, 0.d0, det_alpha_pseudo_SVD, size(det_alpha_pseudo_SVD,1) )
|
||||||
! !!!
|
! !!!
|
||||||
if (elec_beta_num /= 0) then
|
if (elec_beta_num /= 0) then
|
||||||
! det_beta_pseudo_SVD = det_beta_pseudo @ psi_svd_beta / psidet_inv_SVD
|
! det_beta_pseudo_SVD = det_beta_pseudo @ psi_svd_beta * psidet_inv_SVD
|
||||||
call dgemm('N', 'N', elec_beta_num, n_svd_coefs , det_beta_num, 1.d0 &
|
call dgemm('N', 'N', 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), size(psi_svd_beta,1) &
|
, det_beta_pseudo(elec_alpha_num+1,1), size(det_beta_pseudo,1), psi_svd_beta(:,:,1), size(psi_svd_beta,1) &
|
||||||
, 0.d0, det_beta_pseudo_SVD(elec_alpha_num+1,1), size(det_beta_pseudo_SVD,1) )
|
, 0.d0, det_beta_pseudo_SVD(elec_alpha_num+1,1), size(det_beta_pseudo_SVD,1) )
|
||||||
endif
|
endif
|
||||||
@ -1832,13 +1854,15 @@ END_PROVIDER
|
|||||||
! !!!
|
! !!!
|
||||||
do mm = 1, 4
|
do mm = 1, 4
|
||||||
! !!!
|
! !!!
|
||||||
call dgemm('N', 'N', elec_alpha_num, n_svd_coefs , det_alpha_num, 1.d0 &
|
call dgemm('N', 'N', elec_alpha_num, n_svd_coefs, det_alpha_num, 1.d0 &
|
||||||
, det_alpha_grad_lapl(mm,:,:), size(det_alpha_grad_lapl,2), psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) &
|
, det_alpha_grad_lapl(mm,:,:), size(det_alpha_grad_lapl,2) &
|
||||||
|
, psi_svd_alpha(:,:,1), size(psi_svd_alpha,1) &
|
||||||
, 0.d0, det_alpha_grad_lapl_SVD(mm,:,:), size(det_alpha_grad_lapl_SVD,2) )
|
, 0.d0, det_alpha_grad_lapl_SVD(mm,:,:), size(det_alpha_grad_lapl_SVD,2) )
|
||||||
! !!!
|
! !!!
|
||||||
if (elec_beta_num /= 0) then
|
if (elec_beta_num /= 0) then
|
||||||
call dgemm('N', 'N', elec_beta_num, n_svd_coefs , det_beta_num, 1.d0 &
|
call dgemm('N', 'N', elec_beta_num, n_svd_coefs, det_beta_num, 1.d0 &
|
||||||
, det_beta_grad_lapl(mm,:,:), size(det_beta_grad_lapl,2), psi_svd_beta(:,:,1), size(psi_svd_beta,1) &
|
, det_beta_grad_lapl(mm,:,:), size(det_beta_grad_lapl,2) &
|
||||||
|
, psi_svd_beta(:,:,1), size(psi_svd_beta,1) &
|
||||||
, 0.d0, det_beta_grad_lapl_SVD(mm,:,:), size(det_beta_grad_lapl_SVD,2) )
|
, 0.d0, det_beta_grad_lapl_SVD(mm,:,:), size(det_beta_grad_lapl_SVD,2) )
|
||||||
endif
|
endif
|
||||||
! !!!
|
! !!!
|
||||||
@ -1905,7 +1929,8 @@ END_PROVIDER
|
|||||||
do ee = 1, elec_alpha_num
|
do ee = 1, elec_alpha_num
|
||||||
tmp = 0.d0
|
tmp = 0.d0
|
||||||
do l = 1, n_svd_coefs
|
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)
|
tmp = tmp + det_alpha_grad_lapl_SVD(mm,ee,l) * psi_svd_coefs(l,1) &
|
||||||
|
* det_beta_value_SVD(l)
|
||||||
enddo
|
enddo
|
||||||
psidet_grad_lapl_SVD(mm,ee) = tmp
|
psidet_grad_lapl_SVD(mm,ee) = tmp
|
||||||
enddo
|
enddo
|
||||||
@ -1914,7 +1939,8 @@ END_PROVIDER
|
|||||||
do ee = elec_alpha_num+1, elec_num
|
do ee = elec_alpha_num+1, elec_num
|
||||||
tmp = 0.d0
|
tmp = 0.d0
|
||||||
do l = 1, n_svd_coefs
|
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)
|
tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) &
|
||||||
|
* det_beta_grad_lapl_SVD(mm,ee,l)
|
||||||
enddo
|
enddo
|
||||||
psidet_grad_lapl_SVD(mm,ee) = tmp
|
psidet_grad_lapl_SVD(mm,ee) = tmp
|
||||||
enddo
|
enddo
|
||||||
@ -1934,7 +1960,7 @@ END_PROVIDER
|
|||||||
do l = 1, n_svd_coefs
|
do l = 1, n_svd_coefs
|
||||||
tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l)
|
tmp = tmp + det_alpha_pseudo_SVD(ee,l) * psi_svd_coefs(l,1) * det_beta_value_SVD(l)
|
||||||
enddo
|
enddo
|
||||||
pseudo_non_local_SVD(ee) = tmp * psi_value_inv_svd
|
pseudo_non_local_SVD(ee) = tmp * psidet_inv_SVD
|
||||||
enddo
|
enddo
|
||||||
! !!!
|
! !!!
|
||||||
if (elec_beta_num /= 0) then
|
if (elec_beta_num /= 0) then
|
||||||
@ -1943,7 +1969,7 @@ END_PROVIDER
|
|||||||
do l = 1, n_svd_coefs
|
do l = 1, n_svd_coefs
|
||||||
tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) * det_beta_pseudo_SVD(ee,l)
|
tmp = tmp + det_alpha_value_SVD(l) * psi_svd_coefs(l,1) * det_beta_pseudo_SVD(ee,l)
|
||||||
enddo
|
enddo
|
||||||
pseudo_non_local_SVD(ee) = tmp * psi_value_inv_svd
|
pseudo_non_local_SVD(ee) = tmp * psidet_inv_SVD
|
||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
! !!!
|
! !!!
|
||||||
@ -2052,7 +2078,8 @@ BEGIN_PROVIDER [ double precision, det_alpha_grad_lapl_curr, (4,elec_alpha_num)
|
|||||||
! imo = mo_list_alpha_curr(j)
|
! imo = mo_list_alpha_curr(j)
|
||||||
! do i=1,elec_alpha_num
|
! do i=1,elec_alpha_num
|
||||||
! do k=1,4
|
! 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
|
! enddo
|
||||||
! enddo
|
! enddo
|
||||||
|
253
src/opt_Jast/opt_jast.py
Normal file
253
src/opt_Jast/opt_jast.py
Normal file
@ -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.) )
|
||||||
|
# !!!
|
171
src/opt_Jast/opt_jast_Anthony.py
Normal file
171
src/opt_Jast/opt_jast_Anthony.py
Normal file
@ -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()
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
225
src/opt_Jast/opt_jast_freegrad.py
Normal file
225
src/opt_Jast/opt_jast_freegrad.py
Normal file
@ -0,0 +1,225 @@
|
|||||||
|
#!/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')
|
||||||
|
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 = 30
|
||||||
|
total_time_f = 65
|
||||||
|
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/ammar/qp2/src/svdwf/h2o_optJast"
|
||||||
|
# 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_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
|
||||||
|
bnds = [(0.001, 9.99) for _ in range(n_par)]
|
||||||
|
memo_energy = {'fmin': 100.}
|
||||||
|
opt = sp.optimize.minimize(f, x, method="Powell", bounds=bnds
|
||||||
|
, options= {'disp':True,
|
||||||
|
'ftol':0.2,
|
||||||
|
'xtol':0.2,
|
||||||
|
'maxfev':5} )
|
||||||
|
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.) )
|
||||||
|
# !!!
|
301
src/opt_Jast/opt_jast_grad.py
Normal file
301
src/opt_Jast/opt_jast_grad.py
Normal file
@ -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.) )
|
||||||
|
# !!!
|
@ -15,8 +15,8 @@ subroutine draw_init_points
|
|||||||
allocate (do_elec(elec_num), xmin(3,elec_num))
|
allocate (do_elec(elec_num), xmin(3,elec_num))
|
||||||
xmin = -huge(1.)
|
xmin = -huge(1.)
|
||||||
norm = 0.
|
norm = 0.
|
||||||
do i=1,elec_alpha_num
|
do j = 1, ao_num
|
||||||
do j=1,ao_num
|
do i = 1, elec_alpha_num
|
||||||
norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
|
norm += mo_coef_transp(i,j)*mo_coef_transp(i,j)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
325
src/test_SVD/QMCCHEM_withSVD.py
Executable file
325
src/test_SVD/QMCCHEM_withSVD.py
Executable file
@ -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.) )
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
|
|
510
src/test_SVD/h2o/QMCCHEM_withSVD.py
Executable file
510
src/test_SVD/h2o/QMCCHEM_withSVD.py
Executable file
@ -0,0 +1,510 @@
|
|||||||
|
# !!!
|
||||||
|
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)
|
||||||
|
# !!!
|
||||||
|
#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_J_art"
|
||||||
|
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 )
|
||||||
|
# !!!
|
||||||
|
else:
|
||||||
|
print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
|
||||||
|
exit()
|
||||||
|
# !!!
|
||||||
|
print("end after {:.3f} minutes".format((time.time()-t0)/60.) )
|
||||||
|
# !!!
|
||||||
|
# !!!
|
||||||
|
|
||||||
|
|
10
src/test_SVD/h2o/QR.py
Normal file
10
src/test_SVD/h2o/QR.py
Normal file
@ -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)
|
||||||
|
# !!!
|
20
src/test_SVD/h2o/RSVD.py
Normal file
20
src/test_SVD/h2o/RSVD.py
Normal file
@ -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),:]
|
||||||
|
# !!!
|
||||||
|
# !!!
|
Loading…
Reference in New Issue
Block a user