LCPQ
/
qmcchem
mirror of https://gitlab.com/scemama/qmcchem.git
10
1
Fork 1
Code Releases Activity

More precise Ylm

This commit is contained in:
Anthony Scemama 2020-05-24 03:02:00 +02:00
parent e9706cee7a
commit 9b191eaa7b
3 changed files with 102 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
bin/jast_opt.py
View File

@ -7,6 +7,7 @@ import sys
import os
import time
import subprocess
from math import sqrt
QMCCHEM_PATH=os.environ["QMCCHEM_PATH"]
@ -15,7 +16,7 @@ sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/")
from ezfio import ezfio
# PARAMETERS
thresh = 1.e-3
thresh = 1.e-2
block_time = 20
def main():
@ -59,6 +60,17 @@ def main():
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)
@ -70,7 +82,7 @@ def main():
def run_qmc():
subprocess.check_output(['qmcchem', 'run', filename])
return subprocess.check_output(['qmcchem', 'run', filename])
def stop_qmc():
@ -108,30 +120,40 @@ def main():
atexit.register(stop_qmc)
err = thresh+1.
time.sleep(block_time/2)
time.sleep(3.*block_time/4.)
local_thresh = thresh
while err > local_thresh:
time.sleep(block_time)
energy, err = get_energy()
if energy is not None:
print(" %f %f"%(energy, err))
if (energy-2.*err) > memo_energy['fmin']+thresh:
local_thresh = 10.*thresh
else:
err = local_thresh+1.
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: err = 0.
try:
os.kill(pid,0)
except OSError:
print("---")
break
stop_qmc()
os.wait()
memo_energy[h] = energy
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))

6
src/TOOLS/random.irp.f
View File

@ -50,16 +50,16 @@ BEGIN_PROVIDER [ integer*8, seed, (5) ]
END_DOC
integer :: iargc
integer*8 :: i,j
integer*4 :: clock(12)
integer*4 :: clock(33)
double precision :: r
integer*8 :: pid8
read(current_PID,*) pid8
pid8 = iand( shiftl(pid8, 32), pid8)
do i=1,12
do i=1,size(clock)
clock(i) = i
enddo
call system_clock(count=clock(1))
call random_seed(put=clock)
call random_seed(put=clock)
do i=1,5
call random_number(r)
seed(i) = (r-0.5d0)*huge(1_8)

106
src/pseudo.irp.f
View File

@ -183,8 +183,10 @@ BEGIN_PROVIDER [ double precision, v_pseudo_local, (elec_num) ]
cycle
endif
select case (pseudo_n_k(i,k))
case (-3)
rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)
case (-2)
rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)
rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)
case (-1)
rn = nucl_elec_dist_inv(i,l)
case (0)
@ -193,6 +195,8 @@ BEGIN_PROVIDER [ double precision, v_pseudo_local, (elec_num) ]
rn = r
case (2)
rn = r*r
case (3)
rn = r*r*r
case default
rn = r**pseudo_n_k(i,k)
end select
@ -229,18 +233,22 @@ BEGIN_PROVIDER [ double precision, v_pseudo_non_local, (pseudo_non_loc_dim_8,ele
select case (pseudo_n_kl(i,k,l))
case (0)
rn = 1.d0
case (-1)
rn = nucl_elec_dist_inv(i,j)
case (1)
rn = r
case (-1)
rn = nucl_elec_dist_inv(i,j)
case (2)
rn = r*r
rn = r2
case (-2)
rn = nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j)
case (3)
rn = r2*r
case (-3)
rn = nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j)
case default
rn = r**pseudo_n_kl(i,k,l)
end select
tmp(l) += pseudo_v_kl(i,k,l) * exp(-alpha) * rn
tmp(l) += pseudo_v_kl(i,k,l) * dexp(-alpha) * rn
enddo
enddo
do l=0,pseudo_lmax
@ -273,7 +281,7 @@ BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ]
dr_inv = dble(pseudo_grid_size)/pseudo_grid_rmax
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
endif
PROVIDE pseudo_ylm present_mos mo_pseudo_grid_scaled pseudo_non_loc_dim_count
do j=1,elec_num
tmp = 0.d0
@ -299,7 +307,7 @@ BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ]
do l=kk+1,kk+pseudo_non_loc_dim_count(k)
tmp(l,i) = ( ( mo_pseudo_grid_scaled (l,i,n-1) * w1 - &
mo_pseudo_grid_scaled (l,i,n ) * w0 ) * w2 + &
mo_pseudo_grid_scaled (l,i,n+1) * w1 * w0 ) &
mo_pseudo_grid_scaled (l,i,n+1) * w1 * w0 ) &
* pseudo_ylm(l,j)
enddo
enddo
@ -317,16 +325,20 @@ BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ]
END_PROVIDER
real function ylm(l,m,x,y,z,r_inv)
double precision function ylm(l,m,x_,y_,z_,r_inv)
implicit none
include 'constants.F'
BEGIN_DOC
! Y(l,m) evaluated at x,y,z
END_DOC
integer, intent(in) :: l,m
real, intent(in) :: x,y,z,r_inv
real, intent(in) :: x_,y_,z_,r_inv
double precision :: x, y, z
x = x_
y = y_
z = z_
ylm = 0.5 * one_over_sqpi
ylm = 0.5d0 * one_over_sqpi
select case(l)
case(0)
@ -336,45 +348,55 @@ real function ylm(l,m,x,y,z,r_inv)
ylm = ylm * sq3 * r_inv
select case (m)
case (-1)
ylm = ylm * y
case (0)
ylm = ylm * y
case (0)
ylm = ylm * z
case (1)
case (1)
ylm = ylm * x
end select
case(2)
ylm = ylm * r_inv * r_inv * sqrt(5.)
ylm = ylm * r_inv * r_inv * dsqrt(5.d0)
select case (m)
case(-2)
ylm = ylm * sqrt(3.) * x * y
case(-1)
ylm = ylm * sqrt(3.) * y * z
case(0)
ylm = ylm * 0.5 * (2.*z*z - x*x - y*y)
case(1)
ylm = ylm * sqrt(3.) * z * x
case(2)
ylm = ylm * 0.5 * sqrt(3.) * (x*x - y*y)
case(-2)
ylm = ylm * sq3 * x * y
case(-1)
ylm = ylm * sq3 * y * z
case(0)
! ylm = ylm * 0.5d0 * (2.d0*z*z - x*x - y*y)
! ylm = ylm * 0.5d0 * (z*z - x*x + z*z - y*y)
ylm = ylm * 0.5d0 * ( (z+x)*(z-x) + (z+y)*(z-y) )
case(1)
ylm = ylm * sq3 * z * x
case(2)
! ylm = ylm * 0.5d0 * sq3 * (x*x - y*y)
ylm = ylm * 0.5d0 * sq3 * (x+y)*(x-y)
end select
case(3)
ylm = ylm * r_inv * r_inv * r_inv * sqrt(7.)
ylm = ylm * r_inv * r_inv * r_inv * dsqrt(7.d0)
select case (m)
case(-3)
ylm = ylm * 0.25 * sqrt(10.) * (3.*x*x - y*y)
case(-2)
ylm = ylm * sqrt(15.) * x*y*z
case(-1)
ylm = ylm * 0.25*sqrt(6.) * y * (4.*z*z - x*x -y*y)
case(0)
ylm = ylm * 0.5 * z * (2.*z*z - 3.*(x*x + y*y))
case(1)
ylm = ylm * 0.25*sqrt(6.) * x * (4.*z*z - x*x -y*y)
case(2)
ylm = ylm * 0.5*sqrt(15.) * z * (x*x -y*y)
case(3)
ylm = ylm * 0.25*sqrt(10.) * x * (x*x - 3. * y*y)
case(-3)
! ylm = ylm * 0.25d0 * dsqrt(10.d0) * (3.d0*x*x - y*y)
ylm = ylm * 0.25d0 * dsqrt(10.d0) * (2.d0*x*x + (x+y)*(x-y))
case(-2)
ylm = ylm * dsqrt(15.d0) * x*y*z
case(-1)
! ylm = ylm * 0.25d0*dsqrt(6.d0) * y * (4.d0*z*z - x*x -y*y)
ylm = ylm * 0.25d0*dsqrt(6.d0) * y * (2.d0*z*z + (z+x)*(z-x) + (z+y)*(z-y))
case(0)
! ylm = ylm * 0.5d0 * z * (2.d0*z*z - 3.d0*(x*x + y*y))
! ylm = ylm * 0.5d0 * z * ( (z+x)*(z-x) + (z+y)*(z-y) - 2.d0*(x*x + y*y))
ylm = ylm * 0.5d0 * z * ( (z+x)*(z-x) + (z+y)*(z-y) - 2.d0*(x+y)*(x-y) - 4.d0*y*y)
case(1)
! ylm = ylm * 0.25d0*dsqrt(6.d0) * x * (4.d0*z*z - x*x -y*y)
ylm = ylm * 0.25d0*dsqrt(6.d0) * x * (2.d0*z*z + (z-x)*(z+x) + (z-y)*(z+y))
case(2)
! ylm = ylm * 0.5d0*dsqrt(15.d0) * z * (x*x -y*y)
ylm = ylm * 0.5d0*dsqrt(15.d0) * z * (x+y)*(x-y)
case(3)
! ylm = ylm * 0.25d0*dsqrt(10.d0) * x * (x*x - 3.d0 * y*y)
ylm = ylm * 0.25d0*dsqrt(10.d0) * x * ( (x+y)*(x-y) - 2.d0*y*y)
end select
@ -385,7 +407,7 @@ real function ylm(l,m,x,y,z,r_inv)
end select
end
BEGIN_PROVIDER [ double precision, pseudo_ylm, (pseudo_non_loc_dim_8,elec_num) ]
@ -395,14 +417,14 @@ BEGIN_PROVIDER [ double precision, pseudo_ylm, (pseudo_non_loc_dim_8,elec_num) ]
END_DOC
integer :: i,j,l,m,kk
real, external :: ylm
double precision, external :: ylm
integer, save :: ifirst = 0
if (ifirst == 0) then
pseudo_ylm = 0.d0
ifirst = 1
endif
do j=1,elec_num
kk = 0
do i=1,nucl_num