From e9706cee7aae3248ff1cf78ef8380333d66e634e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 11 May 2020 13:49:07 +0200 Subject: [PATCH] Fixed exploding energies with pseudos --- bin/jast_opt.py | 42 +++++++++++++++++------ ocaml/Launcher.ml | 8 +++-- src/SAMPLING/pdmc_step.irp.f | 45 ++++++++++++++++--------- src/SAMPLING/srmc_step.irp.f | 64 ++++++++++++++++++++++++------------ 4 files changed, 110 insertions(+), 49 deletions(-) diff --git a/bin/jast_opt.py b/bin/jast_opt.py index b79b4f5..5bee8b3 100755 --- a/bin/jast_opt.py +++ b/bin/jast_opt.py @@ -8,18 +8,25 @@ import os import time import subprocess -sys.path.insert(0,"/home/scemama/qmcchem/EZFIO/Python/") +QMCCHEM_PATH=os.environ["QMCCHEM_PATH"] + +sys.path.insert(0,QMCCHEM_PATH+"/EZFIO/Python/") from ezfio import ezfio # PARAMETERS thresh = 1.e-3 -block_time = 30 +block_time = 20 def main(): + if len(sys.argv) != 2: + print("Usage: %s "%sys.argv[0]) + sys.exti(1) + filename = sys.argv[1] ezfio.set_file(filename) + def make_atom_map(): labels = {} dimension = 0 @@ -36,10 +43,12 @@ def main(): 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') @@ -50,6 +59,7 @@ def main(): else: return None, None + def set_params_pen(x): x = np.abs(x) y=list(ezfio.jastrow_jast_pen) @@ -58,23 +68,30 @@ def main(): y[j] = x[i] ezfio.set_jastrow_jast_pen(y) + def run_qmc(): 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', - '-m', 'VMC', '-l', str(block_time), - '--time-step=0.3', - '--stop-time=3600', - '--norm=1.e-5', - '-w', '10', filename]) - memo_energy = {} + + memo_energy = {'fmin': 100000000.} def f(x): print ("x = %s"%str(x)) h = str(x) @@ -92,13 +109,16 @@ def main(): err = thresh+1. time.sleep(block_time/2) - while err > thresh: + 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 = thresh+1. + err = local_thresh+1. # Check if PID is still running try: os.kill(pid,0) @@ -106,6 +126,7 @@ def main(): stop_qmc() os.wait() memo_energy[h] = energy + memo_energy['fmin'] = min(energy, memo_energy['fmin']) return energy @@ -115,6 +136,7 @@ def main(): options= {'disp':True, 'ftol':thresh,'xtol':0.02}) print("x = "+str(opt)) set_params_pen(opt['x']) + run() if __name__ == '__main__': diff --git a/ocaml/Launcher.ml b/ocaml/Launcher.ml index 5e84d2b..c456ea5 100644 --- a/ocaml/Launcher.ml +++ b/ocaml/Launcher.ml @@ -6,7 +6,10 @@ type t = let to_string = function | Srun -> "srun" | Bash -> "env" - | MPI -> Lazy.force Qmcchem_config.mpirun + | MPI -> String.concat " " [ Lazy.force Qmcchem_config.mpirun ; + try Sys.getenv "QMCCHEM_MPIRUN_FLAGS" with Not_found -> "" + ] + (* let to_string = function | MPI @@ -16,7 +19,8 @@ let to_string = function | Some p -> p ] | Bash -> "env" -*) + *) + (** Find the launcher for the current job scheduler *) diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index af9b09f..dc62791 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -9,7 +9,7 @@ t = """ &BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ] &BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ] implicit none - BEGIN_DOC + BEGIN_DOC ! PDMC averages of $X. Computed in E_loc_pdmc_block_walk END_DOC $X_pdmc_block_walk = 0.d0 @@ -109,14 +109,14 @@ END_SHELL endif integer :: info -! double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1) +! double precision :: H(0:pdmc_n_diag/2,0:pdmc_n_diag/2), S(0:pdmc_n_diag/2,0:pdmc_n_diag/2), w(0:pdmc_n_diag/2), work(3*pdmc_n_diag+1) ! H = 0.d0 ! S = 0.d0 do while (loop) i_walk = 1 - + if (.not.first_loop) then integer :: i,j,l do l=1,3 @@ -157,7 +157,7 @@ END_SHELL ! 106.d0*E_loc_save(3,i_walk)+19.d0*E_loc_save(4,i_walk))/720.d0 delta = (delta - E_ref)*p - + if (delta >= 0.d0) then pdmc_weight(i_walk) = dexp(-dtime_step*delta) else @@ -183,6 +183,20 @@ END_SHELL E_loc_save(2,i_walk) = E_loc_save(1,i_walk) E_loc_save(1,i_walk) = E_loc + + if (do_print_dmc_data) then + do k=1,walk_num + double precision, external :: qmc_ranf + if (qmc_ranf() < 0.001) then + print *, '--' + do i=1,elec_num + print *, elec_coord_full(i,1:3,k) + enddo + print *, 'w=', pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) + endif + enddo + endif + if (dabs(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag)) > 1.d-15) then dmc_zv_weight = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag)) dmc_zv_weight_half = 1.d0/(pdmc_weight(i_walk)*pdmc_pop_weight_mult(pdmc_n_diag/2)) @@ -205,19 +219,19 @@ t = """ ! $X_pdmc_block_walk += $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) ! $X_2_pdmc_block_walk += $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - + $X_pdmc_block_walk_kahan($D2 3) = $X * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_pdmc_block_walk_kahan($D2 1) $X_pdmc_block_walk_kahan($D2 2) = $X_pdmc_block_walk $D1 + $X_pdmc_block_walk_kahan($D2 3) $X_pdmc_block_walk_kahan($D2 1) = ($X_pdmc_block_walk_kahan($D2 2) - $X_pdmc_block_walk $D1 ) & - $X_pdmc_block_walk_kahan($D2 3) - $X_pdmc_block_walk $D1 = $X_pdmc_block_walk_kahan($D2 2) - - + $X_pdmc_block_walk $D1 = $X_pdmc_block_walk_kahan($D2 2) + + $X_2_pdmc_block_walk_kahan($D2 3) = $X_2 * pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - $X_2_pdmc_block_walk_kahan($D2 1) $X_2_pdmc_block_walk_kahan($D2 2) = $X_2_pdmc_block_walk $D1 + $X_2_pdmc_block_walk_kahan($D2 3) $X_2_pdmc_block_walk_kahan($D2 1) = ($X_2_pdmc_block_walk_kahan($D2 2) - $X_2_pdmc_block_walk $D1 ) & - $X_2_pdmc_block_walk_kahan($D2 3) - $X_2_pdmc_block_walk $D1 = $X_2_pdmc_block_walk_kahan($D2 2) + $X_2_pdmc_block_walk $D1 = $X_2_pdmc_block_walk_kahan($D2 2) endif """ for p in properties: @@ -247,7 +261,6 @@ END_SHELL ! pdmc_pop_weight(:,:) = 1.d0 ! pdmc_pop_weight_mult(:) = 1.d0 ! endif - do k=1,pdmc_n_diag ! Move to the next projection step @@ -266,7 +279,7 @@ END_SHELL endif ! Remove contribution of the old value of the weight at the new ! projection step - + pdmc_pop_weight_mult(k) *= 1.d0/pdmc_pop_weight(pdmc_projection_step(k),k) pdmc_pop_weight(pdmc_projection_step(k),k) = pdmc_weight(i_walk)/dble(walk_num) @@ -288,13 +301,13 @@ END_SHELL cpu2 = cpu1 endif - SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult + SOFT_TOUCH elec_coord_full pdmc_pop_weight_mult first_loop = .False. enddo - + double precision :: factor factor = 1.d0/block_weight SOFT_TOUCH block_weight @@ -337,7 +350,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ] &BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ] implicit none - BEGIN_DOC + BEGIN_DOC ! Number of projection steps for PDMC END_DOC real :: pdmc_projection_time @@ -353,7 +366,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ] implicit none - BEGIN_DOC + BEGIN_DOC ! Population weight of PDMC END_DOC pdmc_pop_weight(:,:) = 1.d0 @@ -361,7 +374,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (0:pdmc_n_diag) ] implicit none - BEGIN_DOC + BEGIN_DOC ! Population weight of PDMC END_DOC pdmc_pop_weight_mult(:) = 1.d0 diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index 79690f7..01b6686 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -1,5 +1,14 @@ ! Providers of *_srmc_block_walk !============================== + +BEGIN_PROVIDER [ logical, do_print_dmc_data ] + implicit none + BEGIN_DOC + ! If true, print in stdout the data to fit a Jastrow + END_DOC + do_print_dmc_data = .False. +END_PROVIDER + BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * @@ -9,7 +18,7 @@ t = """ &BEGIN_PROVIDER [ $T, $X_2_srmc_block_walk $D1 ] &BEGIN_PROVIDER [ $T, $X_2_srmc_block_walk_kahan $D2 ] implicit none - BEGIN_DOC + BEGIN_DOC ! SRMC averages of $X. Computed in E_loc_srmc_block_walk END_DOC $X_srmc_block_walk = 0.d0 @@ -109,7 +118,7 @@ END_SHELL ! Every walker makes a step do i_walk=1,walk_num - + if (.not.first_loop) then integer :: i,j,l do l=1,3 @@ -155,12 +164,12 @@ END_SHELL if (delta >= 0.d0) then srmc_weight(i_walk) = dexp(-dtime_step*delta) else - srmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta) + srmc_weight(i_walk) = max(2.d0-dexp(dtime_step*delta), 0.d0) endif ! Trick to avoid holes in DMC PES. - if (dabs(delta/E_ref) * time_step_sq > p ) then - srmc_weight(i_walk) = 1.d-1 + if (dabs(delta/E_ref) * time_step_sq > p * 0.5d0 ) then + srmc_weight(i_walk) = 0.d0 endif @@ -169,23 +178,23 @@ END_SHELL ! double precision :: delta_old, delta_new ! delta_old = (9.d0*E_loc_save(1,i_walk)+19.d0*E_loc_save(2,i_walk)-& ! 5.d0*E_loc_save(3,i_walk)+E_loc_save(4,i_walk))/24.d0 - E_ref -! -! +! +! ! if (delta_old >= 0.d0) then ! srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(dtime_step*delta_old) ! else ! srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(-dtime_step*delta_old)) ! endif -! +! ! delta_new = (-(E_loc_save_tmp(3,i_walk)-13.d0*E_loc_save_tmp(2,i_walk)& ! -13.d0*E_loc_save_tmp(1,i_walk)+E_loc))/24.d0 - E_ref -! +! ! if (delta_new >= 0.d0) then ! srmc_weight(i_walk) = srmc_weight(i_walk) * dexp(-dtime_step*delta_new) ! else ! srmc_weight(i_walk) = srmc_weight(i_walk) * (2.d0-dexp(dtime_step*delta_new) ) ! endif -! +! ! endif @@ -222,19 +231,19 @@ t = """ ! $X_srmc_block_walk += $X * srmc_pop_weight_mult * srmc_weight(i_walk) ! $X_2_srmc_block_walk += $X_2 * srmc_pop_weight_mult * srmc_weight(i_walk) ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - + $X_srmc_block_walk_kahan($D2 3) = $X * srmc_pop_weight_mult * srmc_weight(i_walk) - $X_srmc_block_walk_kahan($D2 1) $X_srmc_block_walk_kahan($D2 2) = $X_srmc_block_walk $D1 + $X_srmc_block_walk_kahan($D2 3) $X_srmc_block_walk_kahan($D2 1) = ($X_srmc_block_walk_kahan($D2 2) - $X_srmc_block_walk $D1 ) & - $X_srmc_block_walk_kahan($D2 3) - $X_srmc_block_walk $D1 = $X_srmc_block_walk_kahan($D2 2) - - + $X_srmc_block_walk $D1 = $X_srmc_block_walk_kahan($D2 2) + + $X_2_srmc_block_walk_kahan($D2 3) = $X_2 * srmc_pop_weight_mult * srmc_weight(i_walk) - $X_2_srmc_block_walk_kahan($D2 1) $X_2_srmc_block_walk_kahan($D2 2) = $X_2_srmc_block_walk $D1 + $X_2_srmc_block_walk_kahan($D2 3) $X_2_srmc_block_walk_kahan($D2 1) = ($X_2_srmc_block_walk_kahan($D2 2) - $X_2_srmc_block_walk $D1 ) & - $X_2_srmc_block_walk_kahan($D2 3) - $X_2_srmc_block_walk $D1 = $X_2_srmc_block_walk_kahan($D2 2) + $X_2_srmc_block_walk $D1 = $X_2_srmc_block_walk_kahan($D2 2) endif """ for p in properties: @@ -253,7 +262,7 @@ END_SHELL else srmc_weight(i_walk) = 0.d0 endif - + enddo @@ -294,6 +303,19 @@ END_SHELL ! Update the running population weight srmc_pop_weight_mult *= srmc_pop_weight(srmc_projection_step) + if (do_print_dmc_data) then + do k=1,walk_num + double precision, external :: qmc_ranf + if (qmc_ranf() < 0.001) then + print *, '--' + do i=1,elec_num + print *, elec_coord_full(i,1:3,k) + enddo + print *, 'w=', srmc_weight(k) + endif + enddo + endif + ! Reconfiguration integer :: ipos(walk_num) @@ -302,7 +324,7 @@ END_SHELL enddo call dsort(srmc_weight,ipos,walk_num) call reconfigure(ipos,srmc_weight) - + do k=1,walk_num do l=1,3 do i=1,elec_num+1 @@ -315,7 +337,7 @@ END_SHELL psi_value_save_tmp(k) = psi_value_save(k) E_loc_save_tmp(:,k) = E_loc_save(:,k) enddo - + integer :: ipm do k=1,walk_num ipm = ipos(k) @@ -372,7 +394,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, srmc_pop_weight_mult ] implicit none - BEGIN_DOC + BEGIN_DOC ! Population weight of SRMC END_DOC srmc_pop_weight_mult = srmc_pop_weight(srmc_projection) @@ -381,7 +403,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer, srmc_projection ] &BEGIN_PROVIDER [ integer, srmc_projection_step ] implicit none - BEGIN_DOC + BEGIN_DOC ! Number of projection steps for SRMC END_DOC real :: srmc_projection_time @@ -393,7 +415,7 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, srmc_pop_weight, (0:srmc_projection+1) ] implicit none - BEGIN_DOC + BEGIN_DOC ! Population weight of SRMC END_DOC srmc_pop_weight = 1.d0