From 9b191eaa7b4d727239d59215ef7049003d2b37e7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 24 May 2020 03:02:00 +0200 Subject: [PATCH] More precise Ylm --- bin/jast_opt.py | 48 ++++++++++++++----- src/TOOLS/random.irp.f | 6 +-- src/pseudo.irp.f | 106 +++++++++++++++++++++++++---------------- 3 files changed, 102 insertions(+), 58 deletions(-) diff --git a/bin/jast_opt.py b/bin/jast_opt.py index 5bee8b3..c13efad 100755 --- a/bin/jast_opt.py +++ b/bin/jast_opt.py @@ -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)) diff --git a/src/TOOLS/random.irp.f b/src/TOOLS/random.irp.f index aa99839..ea28d9f 100644 --- a/src/TOOLS/random.irp.f +++ b/src/TOOLS/random.irp.f @@ -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) diff --git a/src/pseudo.irp.f b/src/pseudo.irp.f index 156c25a..23273be 100644 --- a/src/pseudo.irp.f +++ b/src/pseudo.irp.f @@ -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