LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-12-21 11:53:30 +01:00
Code Releases Activity

Fixed exploding energies with pseudos

Browse Source
This commit is contained in:
Anthony Scemama 2020-05-11 13:49:07 +02:00
parent f41412658e
commit e9706cee7a
4 changed files with 110 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
bin/jast_opt.py
View File

@ -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 <EZFIO_DIRECTORY>"%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__':

8
ocaml/Launcher.ml
View File

@ -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 *)

45
src/SAMPLING/pdmc_step.irp.f
View File

@ -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

64
src/SAMPLING/srmc_step.irp.f
View File

@ -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