LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-12-22 04:13:31 +01:00
Code Releases Activity

Added 1body jastrow mu erf

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-26 19:48:16 +01:00
parent 557ef562cb
commit 1b12940f93
8 changed files with 183 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
bin/jast_opt.py
View File

@ -29,6 +29,8 @@ def main():
filename = sys.argv[1] filename = sys.argv[1]
ezfio.set_file(filename) ezfio.set_file(filename)
jast_type = ezfio.jastrow_jast_type
print (jast_type)
def make_atom_map(): def make_atom_map():
labels = {} labels = {}
@ -47,18 +49,62 @@ def main():
atom_map = make_atom_map() atom_map = make_atom_map()
if jast_type == 'Simple':
def get_params_b(): def get_params_b():
return ezfio.jastrow_jast_b_up_dn return ezfio.jastrow_jast_b_up_dn
def set_params_b(x):
x = np.abs(x)
ezfio.set_jastrow_jast_b_up_up(x)
ezfio.set_jastrow_jast_b_up_dn(x)
elif jast_type == 'Mu':
def get_params_b():
return ezfio.jastrow_mu_erf
def set_params_b(x):
x = np.abs(x)
ezfio.set_jastrow_mu_erf(x)
if jast_type == 'Simple':
def get_params_pen(): def get_params_pen():
d = ezfio.jastrow_jast_pen d = ezfio.jastrow_jast_pen
print(atom_map)
for m in atom_map:
print(m[0])
print (d[m[0]])
sys.stdout.flush()
return np.array([d[m[0]] for m in atom_map]) 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)
elif jast_type == 'Mu':
def get_params_pen():
d = ezfio.jastrow_jast_1bgauss_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_1bgauss_pen)
for i,m in enumerate(atom_map):
for j in m:
y[j] = x[i]
ezfio.set_jastrow_jast_1bgauss_pen(y)
def get_norm():
return 1.,0.
# buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'psi_norm', filename],
# encoding='UTF-8')
# if buffer.strip() != "":
# buffer = buffer.splitlines()[-1]
# _, energy, error = [float(x) for x in buffer.split()]
# else:
# return None, None
def get_energy(): def get_energy():
buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename], buffer = subprocess.check_output(['qmcchem', 'result', '-e', 'e_loc', filename],
@ -82,19 +128,6 @@ def main():
return None, None return None, None
def set_params_b(x):
x = np.abs(x)
ezfio.set_jastrow_jast_b_up_up(x)
ezfio.set_jastrow_jast_b_up_dn(x)
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(): def run_qmc():
return subprocess.check_output(['qmcchem', 'run', filename]) return subprocess.check_output(['qmcchem', 'run', filename])
@ -105,7 +138,7 @@ def main():
def set_vmc_params(): def set_vmc_params():
# subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple', # subprocess.check_output(['qmcchem', 'edit', '-c', '-j', jast_type,
# '-m', 'VMC', # '-m', 'VMC',
# '-l', str(block_time), # '-l', str(block_time),
# '--time-step=0.3', # '--time-step=0.3',
@ -113,7 +146,7 @@ def main():
# '--norm=1.e-5', # '--norm=1.e-5',
# '-w', '10', # '-w', '10',
# filename]) # filename])
subprocess.check_output(['qmcchem', 'edit', '-c', '-j', 'Simple', subprocess.check_output(['qmcchem', 'edit', '-c', '-j', jast_type,
'-l', str(block_time), '-l', str(block_time),
filename]) filename])
@ -145,6 +178,9 @@ 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
# norm, _ = get_norm()
# e, e_err = e, e_err/norm
# variance, v_err =variance/norm, v_err/norm
energy = e #+ variance energy = e #+ variance
err = e_err #sqrt(e_err*e_err+v_err*v_err) err = e_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)) print(" %f %f %f %f %f %f"%(e, e_err, variance, v_err, energy, err))
@ -173,7 +209,9 @@ def main():
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", # opt = sp.optimize.minimize(f,x,method="Powell",
# options= {'disp':True, 'ftol':thresh,'xtol':0.02})
opt = sp.optimize.minimize(f,x,method="Nelder-Mead",
options= {'disp':True, 'ftol':thresh,'xtol':0.02}) options= {'disp':True, 'ftol':thresh,'xtol':0.02})
print("x = "+str(opt)) print("x = "+str(opt))
sys.stdout.flush() sys.stdout.flush()

11
src/JASTROW/jastrow_1b.irp.f
View File

@ -49,8 +49,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_value, (elec_num_8) ]
elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow
! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ] ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
! mu = jast_mu_erf mu = jast_mu_erf
mu = mu_erf
mu_pi = 1.d0 / ( dsqpi * mu ) mu_pi = 1.d0 / ( dsqpi * mu )
!DIR$ LOOP COUNT (100) !DIR$ LOOP COUNT (100)
do j = 1, nucl_num do j = 1, nucl_num
@ -151,8 +150,7 @@ END_PROVIDER
elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow
! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ] ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
! mu = jast_mu_erf mu = jast_mu_erf
mu = mu_erf
!DIR$ LOOP COUNT (100) !DIR$ LOOP COUNT (100)
do j = 1, nucl_num do j = 1, nucl_num
rij = nucl_elec_dist(j,i) rij = nucl_elec_dist(j,i)
@ -248,8 +246,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_lapl, (elec_num_8) ]
elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow elseif( jast_1b_type .eq. 4 ) then ! add 1body-RSDFT Jastrow
! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ] ! J(i) = - \sum_A [ -z_A r_iA erfc(mu*r_iA) + z_A exp(-(mu*r_iA)^2)/(mu*sqt_pi) ]
! mu = jast_mu_erf mu = jast_mu_erf
mu = mu_erf
mu_pi = mu / dsqpi mu_pi = mu / dsqpi
!DIR$ LOOP COUNT (100) !DIR$ LOOP COUNT (100)
do j = 1, nucl_num do j = 1, nucl_num
@ -278,7 +275,7 @@ BEGIN_PROVIDER [ double precision, jast_1b_lapl, (elec_num_8) ]
a = jast_1bGauss_pen(j) a = jast_1bGauss_pen(j)
rij = nucl_elec_dist(j,i) rij = nucl_elec_dist(j,i)
c = a * rij * rij c = a * rij * rij
tmp = 2.d0 * a * dexp(-c) * (3.d0-2.d0*a*c) tmp = 2.d0 * a * dexp(-c) * (3.d0-2.d0*c)
jast_1b_lapl(i) -= tmp jast_1b_lapl(i) -= tmp
enddo enddo

27
src/JASTROW/jastrow_mu.irp.f
View File

@ -13,10 +13,10 @@ implicit none
double precision :: a, b, rij, tmp double precision :: a, b, rij, tmp
include '../constants.F' include '../constants.F'
double precision :: mu double precision :: mu
mu = mu_erf mu = jast_mu_erf
do i=1,elec_num do i=1,elec_num
jast_elec_Mu_value(i) = 0.d0 jast_elec_Mu_value(i) = jast_1b_value(i)
enddo enddo
do j=1,elec_num do j=1,elec_num
@ -25,10 +25,9 @@ implicit none
if(j==i)cycle if(j==i)cycle
rij = elec_dist(i,j) rij = elec_dist(i,j)
tmp = 0.5d0 * rij * (1.d0 - derf(mu*rij)) - 0.5d0/(dsqpi*mu) * dexp(-mu*mu*rij*rij) tmp = 0.5d0 * rij * (1.d0 - derf(mu*rij)) - 0.5d0/(dsqpi*mu) * dexp(-mu*mu*rij*rij)
jast_elec_Mu_value(i) += tmp jast_elec_Mu_value(j) += 0.5d0*tmp ! symmetrization
enddo enddo
enddo enddo
jast_elec_Mu_value = jast_elec_Mu_value * 0.5d0 ! symmetrization
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_x, (elec_num_8) ] BEGIN_PROVIDER [ double precision , jast_elec_Mu_grad_x, (elec_num_8) ]
@ -44,13 +43,12 @@ END_PROVIDER
double precision :: a, b, rij, tmp, x, y, z double precision :: a, b, rij, tmp, x, y, z
include '../constants.F' include '../constants.F'
double precision :: mu double precision :: mu
mu = mu_erf mu = jast_mu_erf
do i=1,elec_num do i=1,elec_num
jast_elec_Mu_grad_x(i) = 0.d0 jast_elec_Mu_grad_x(i) = jast_1b_grad_x(i)
jast_elec_Mu_grad_y(i) = 0.d0 jast_elec_Mu_grad_y(i) = jast_1b_grad_y(i)
jast_elec_Mu_grad_z(i) = 0.d0 jast_elec_Mu_grad_z(i) = jast_1b_grad_z(i)
!DIR$ LOOP COUNT (100)
enddo enddo
! (grad of J(r12) with respect to xi, yi, zi) ! (grad of J(r12) with respect to xi, yi, zi)
@ -76,10 +74,10 @@ BEGIN_PROVIDER [ double precision , jast_elec_Mu_lapl, (elec_num_8) ]
double precision :: a, b, rij, tmp, x, y, z double precision :: a, b, rij, tmp, x, y, z
include '../constants.F' include '../constants.F'
double precision :: mu, x_ij, y_ij, z_ij, rij_inv double precision :: mu, x_ij, y_ij, z_ij, rij_inv
mu = mu_erf mu = jast_mu_erf
do i=1,elec_num do i=1,elec_num
jast_elec_Mu_lapl(i) = 0.d0 jast_elec_Mu_lapl(i) = jast_1b_lapl(i)
enddo enddo
do i=1, elec_num do i=1, elec_num
@ -96,10 +94,6 @@ BEGIN_PROVIDER [ double precision , jast_elec_Mu_lapl, (elec_num_8) ]
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, mu_erf ]
implicit none
mu_erf = 0.5d0
END_PROVIDER
BEGIN_PROVIDER [double precision, grad_j_mu_x,(elec_num, elec_num)] BEGIN_PROVIDER [double precision, grad_j_mu_x,(elec_num, elec_num)]
&BEGIN_PROVIDER [double precision, grad_j_mu_y,(elec_num, elec_num)] &BEGIN_PROVIDER [double precision, grad_j_mu_y,(elec_num, elec_num)]
@ -110,7 +104,7 @@ END_PROVIDER
END_DOC END_DOC
integer :: i,j integer :: i,j
double precision :: rij, mu,scal double precision :: rij, mu,scal
mu = mu_erf mu = jast_mu_erf
grad_j_mu_x = 0.d0 grad_j_mu_x = 0.d0
grad_j_mu_y = 0.d0 grad_j_mu_y = 0.d0
grad_j_mu_z = 0.d0 grad_j_mu_z = 0.d0
@ -126,3 +120,4 @@ END_PROVIDER
enddo enddo
END_PROVIDER END_PROVIDER

14
src/JASTROW/jastrow_param.irp.f
View File

@ -156,16 +156,24 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ real, jast_mu_erf ]
implicit none
include '../types.F'
jast_mu_erf = 0.5
call get_jastrow_mu_erf(jast_mu_erf)
END_PROVIDER
BEGIN_PROVIDER [ integer, jast_1b_type ] BEGIN_PROVIDER [ integer, jast_1b_type ]
implicit none implicit none
include '../types.F' include '../types.F'
jast_1b_type = 0 ! no 1body Jastrow ! jast_1b_type = 0 ! no 1body Jastrow
!jast_1b_type = 2 ! add 1body-Tanh Jastrow !jast_1b_type = 2 ! add 1body-Tanh Jastrow
!jast_1b_type = 3 ! add 1body-Simple Jastrow !jast_1b_type = 3 ! add 1body-Simple Jastrow
!jast_1b_type = 4 ! add 1body-RSDFT Jastrow !jast_1b_type = 4 ! add 1body-RSDFT Jastrow
!jast_1b_type = 5 ! add 1body-erf Jastrow !jast_1b_type = 5 ! add 1body-erf Jastrow
!jast_1b_type = 6 ! add 1body-Gauss Jastrow jast_1b_type = 6 ! add 1body-Gauss Jastrow
call get_jastrow_jast_1b_type(jast_1b_type) ! call get_jastrow_jast_1b_type(jast_1b_type)
END_PROVIDER END_PROVIDER
! useful if jast_1b_type = 2 ! useful if jast_1b_type = 2

5
src/Makefile
View File

@ -5,3 +5,8 @@ export
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile
$(IRPF90) $(IRPF90)
IRPF90_temp/irp_touches.irp.o: IRPF90_temp/irp_touches.irp.F90
$(FC) -O2 -c -g IRPF90_temp/irp_touches.irp.F90 -o IRPF90_temp/irp_touches.irp.o

68
src/PROPERTIES/properties_ci.irp.f
View File

@ -331,3 +331,71 @@ return
ci_dress_max = max(ci_dress_max,maxval(ci_dress)) ci_dress_max = max(ci_dress_max,maxval(ci_dress))
SOFT_TOUCH ci_dress_min ci_dress_max SOFT_TOUCH ci_dress_min ci_dress_max
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_dress_mu_opt ]
BEGIN_DOC
! Use for optimizing mu
END_DOC
implicit none
integer :: i, j, k, l
double precision :: T, dij, f, E_noJ, dE
! E_noJ = -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl
! ci_dress_mu_opt = (E_loc - energy_mu) * psi_value_inv * jast_value_inv * &
! det_alpha_value(1) * det_beta_value(1)
! energy_mu = H_mu \Phi / \Phi
dE = (E_loc - energy_mu) * psi_value_inv * jast_value_inv
! dE = (E_loc - E_noJ) * psi_value_inv * jast_value_inv
k = 1
i = det_coef_matrix_rows( k)
j = det_coef_matrix_columns(k)
f = det_alpha_value(i) * det_beta_value(j)
ci_dress_mu_opt = dE * f
ci_dress_mu_opt_min = min(ci_dress_mu_opt_min, ci_dress_mu_opt)
ci_dress_mu_opt_max = max(ci_dress_mu_opt_max, ci_dress_mu_opt)
SOFT_TOUCH ci_dress_mu_opt_min ci_dress_mu_opt_max
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_dress_Htilde, (size_ci_dress) ]
implicit none
BEGIN_DOC
! < det(i) e^{-J} |H| Psi >
!
! Dimensions : det_num
END_DOC
integer :: i, j, k, l
double precision :: T, h_psidet, dij, f, E_noJ, dE
h_psidet = -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl
E_noJ = h_psidet
dE = E_loc - E_noJ
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
f = det_alpha_value(i)*det_beta_value(j) * psi_value_inv * jast_value_inv
ci_dress(k) = dE * f
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ci_dress_H, (size_ci_dress) ]
implicit none
BEGIN_DOC
! < det(i) e^{-J} |H| Psi >
!
! Dimensions : det_num
END_DOC
integer :: i, j, k, l
double precision :: T, h_psidet, dij, f, E_noJ, dE
E_noJ= -0.5d0*psidet_lapl*psidet_inv + E_pot + E_nucl
do k=1,det_num
i = det_coef_matrix_rows(k)
j = det_coef_matrix_columns(k)
f = det_alpha_value(i)*det_beta_value(j) * psi_value_inv * jast_value_inv
ci_dress(k) = E_noJ * f
enddo
END_PROVIDER

34
src/PROPERTIES/properties_mu.irp.f
View File

@ -57,8 +57,7 @@ END_PROVIDER
integer :: i,j integer :: i,j
double precision :: rij, mu double precision :: rij, mu
mu = mu_erf mu = jast_mu_erf
! mu = jast_mu_erf
Eff_pot_mu_elec = 0.d0 Eff_pot_mu_elec = 0.d0
@ -85,18 +84,21 @@ END_PROVIDER
enddo enddo
! 1-body Jastrow ! 1-body Jastrow
! if( jast_1b_type .gt. 0 ) then if( jast_1b_type .gt. 0 ) then
! do i = 1, elec_num do i = 1, elec_num
! Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_lapl(i) Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_lapl(i)
! Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_grad_sq(i) Eff_pot_mu_elec(i) -= 0.5d0 * jast_1b_grad_sq(i)
! do j = 1, elec_num do j = 1, elec_num
! if(i==j) cycle if(i==j) cycle
! Eff_pot_mu_elec(i) -= ( jast_elec_Mu_grad_x(i) * ( jast_1b_grad_x(i) - jast_1b_grad_x(j) ) & ! + sign for i <--> j
! + jast_elec_Mu_grad_y(i) * ( jast_1b_grad_y(i) - jast_1b_grad_y(j) ) & ! 0.5d0 for double counting
! + jast_elec_Mu_grad_z(i) * ( jast_1b_grad_z(i) - jast_1b_grad_z(j) ) ) Eff_pot_mu_elec(i) += 0.5d0 * &
! enddo ( grad_j_mu_x(j,i) * ( jast_1b_grad_x(i) - jast_1b_grad_x(j) ) &
! enddo + grad_j_mu_y(j,i) * ( jast_1b_grad_y(i) - jast_1b_grad_y(j) ) &
! endif + grad_j_mu_z(j,i) * ( jast_1b_grad_z(i) - jast_1b_grad_z(j) ) )
enddo
enddo
endif
END_PROVIDER END_PROVIDER
@ -143,8 +145,7 @@ BEGIN_PROVIDER [double precision, Eff_pot_deriv_mu_elec, (elec_num) ]
integer :: i, j integer :: i, j
double precision :: rij, mu double precision :: rij, mu
! mu = jast_mu_erf mu = jast_mu_erf
mu = mu_erf
Eff_pot_deriv_mu_elec = 0.d0 Eff_pot_deriv_mu_elec = 0.d0
! 2body-Jastrow: (eq A4) ! 2body-Jastrow: (eq A4)
@ -237,3 +238,4 @@ BEGIN_PROVIDER [ double precision, ci_dress_mu, (size_ci_dress_mu) ]
SOFT_TOUCH ci_dress_mu_min ci_dress_mu_max SOFT_TOUCH ci_dress_mu_min ci_dress_mu_max
END_PROVIDER END_PROVIDER