From 4edcec78db76021a0c0564809ce8f2efaed8afdb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 2 May 2018 21:26:22 +0200 Subject: [PATCH] python2 --- build.ninja | 2 +- install/build.ninja | 8 +- install/scripts/install_f77_zmq.sh | 2 +- install/scripts/install_zmq.sh | 15 +- ocaml/Qmcchem_dataserver.ml | 12 +- ocaml/Qmcchem_result.ml | 40 -- ocaml/compile.sh | 2 +- ocaml/ninja_ocaml.py | 4 +- ocaml/qmcchem.ml | 1 + scripts/compile_irpf90.sh | 2 +- scripts/compile_ocaml_dep.sh | 2 +- scripts/create_properties_ezfio.py | 2 +- scripts/create_properties_python.py | 2 +- src/PROPERTIES/properties.irp.f | 4 +- src/PROPERTIES/properties_energy.irp.f | 144 +----- src/PROPERTIES/properties_general.irp.f | 16 - src/SAMPLING/block.irp.f | 2 +- src/SAMPLING/dmc_step.irp.f | 8 +- src/SAMPLING/fkmc_step.irp.f | 8 +- src/SAMPLING/pdmc_step.irp.f | 627 +++++++++++------------- src/SAMPLING/srmc_step.irp.f | 23 +- src/SAMPLING/vmc_step.irp.f | 11 +- src/TOOLS/random.irp.f | 2 +- src/ZMQ/qmc.irp.f | 2 +- src/ezfio_interface.irp.f | 2 +- src/properties.py | 2 +- 26 files changed, 369 insertions(+), 576 deletions(-) diff --git a/build.ninja b/build.ninja index 41243c7..6cbb168 100644 --- a/build.ninja +++ b/build.ninja @@ -4,7 +4,7 @@ rule compile_ezfio pool = console rule build_properties_config - command = bash -c "source qmcchemrc ; exec python ./scripts/create_properties_ezfio.py" + command = bash -c "source qmcchemrc ; exec ./scripts/create_properties_ezfio.py" pool = console rule compile_irpf90 diff --git a/install/build.ninja b/install/build.ninja index 34ab4f7..2acb771 100644 --- a/install/build.ninja +++ b/install/build.ninja @@ -4,11 +4,11 @@ ###### URL_OPAM ="https://raw.github.com/ocaml/opam/master/shell/opam_installer.sh" -URL_IRPF90="https://github.com/scemama/irpf90/archive/v1.6.7.tar.gz" -URL_EZFIO ="https://github.com/scemama/EZFIO/archive/v1.3.1.tar.gz" +URL_IRPF90="https://github.com/scemama/irpf90/archive/v1.7.2.tar.gz" +URL_EZFIO ="https://github.com/scemama/EZFIO/archive/v1.3.2.tar.gz" -URL_ZMQ ="http://download.zeromq.org/zeromq-4.1.4.tar.gz" -URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/4.1.4.tar.gz" +URL_ZMQ ="http://github.com/zeromq/libzmq/releases/download/v4.2.5/zeromq-4.2.5.tar.gz" +URL_F77ZMQ="https://github.com/scemama/f77_zmq/archive/v4.2.5.tar.gz" # Rules ####### diff --git a/install/scripts/install_f77_zmq.sh b/install/scripts/install_f77_zmq.sh index a996788..bcc1732 100755 --- a/install/scripts/install_f77_zmq.sh +++ b/install/scripts/install_f77_zmq.sh @@ -15,7 +15,7 @@ function _install() cd - rm -f -- "${QMCCHEM_PATH}"/src/ZMQ/f77_zmq.h "${QMCCHEM_PATH}"/lib/libf77zmq.a "${QMCCHEM_PATH}"/lib/libf77zmq.so cp "${BUILD}"/libf77zmq.{a,so} ../lib/ - cp "${BUILD}"/f77_zmq.h ../src/ZMQ/ + cp "${BUILD}"/f77_zmq_free.h ../src/ZMQ/ f77_zmq.h return 0 } diff --git a/install/scripts/install_zmq.sh b/install/scripts/install_zmq.sh index 6c7609c..be26568 100755 --- a/install/scripts/install_zmq.sh +++ b/install/scripts/install_zmq.sh @@ -10,18 +10,11 @@ function _install() set -e set -u cd "${BUILD}" - ./configure --without-libsodium - make -j 8 - cd - rm -f -- "${QMCCHEM_PATH}"/lib/libzmq.{a,so,so.$LIBVERSION} - cp "${BUILD}"/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ - cp "${BUILD}"/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION -# cp "${BUILD}"/src/.libs/libzmq.a "${QMCCHEM_PATH}"/lib/ -# cp "${BUILD}"/src/.libs/libzmq.so "${QMCCHEM_PATH}"/lib/libzmq.so.$LIBVERSION - cp "${BUILD}"/include/{zmq,zmq_utils}.h "${QMCCHEM_PATH}"/lib/ - cd "${QMCCHEM_PATH}"/lib - ln libzmq.so.$LIBVERSION libzmq.so || cp libzmq.so.$LIBVERSION libzmq.so - cd - + ./configure --without-libsodium --enable-libunwind=no --prefix="${QMCCHEM_PATH}" + make -j 8 + make install + cd - return 0 } diff --git a/ocaml/Qmcchem_dataserver.ml b/ocaml/Qmcchem_dataserver.ml index e4ce5c1..b870de8 100644 --- a/ocaml/Qmcchem_dataserver.ml +++ b/ocaml/Qmcchem_dataserver.ml @@ -59,12 +59,10 @@ let run ?(daemon=true) ezfio_filename = -(* (** Maximum size of the blocks file before compressing *) let max_file_size = ref ( Byte_units.create `Kilobytes 64.) in -*) let hostname = @@ -379,18 +377,20 @@ let run ?(daemon=true) ezfio_filename = same host name, but different PIDs. The result is merging all the CPU cores of the compute nodes. Happens when [max_file_size] is reached. *) - (* + let compress_block_file filename = let t0 = Time.now () in Out_channel.close !block_channel; Unix.rename ~src:block_channel_filename_locked ~dst:block_channel_filename_tmp; + (* Random_variable.compress_files (); + *) send_log "status" 0 t0 "Compressed block file"; block_channel := Out_channel.create ~append:true block_channel_filename_locked in - *) + @@ -785,7 +785,6 @@ let run ?(daemon=true) ezfio_filename = | _ -> begin Out_channel.flush !block_channel ; - (* let file_size = (Unix.stat block_channel_filename_locked).Unix.st_size |> Float.of_int64 @@ -796,7 +795,6 @@ let run ?(daemon=true) ezfio_filename = compress_block_file (); max_file_size := Byte_units.scale file_size 1.2; end - *) end end done; @@ -819,9 +817,7 @@ let run ?(daemon=true) ezfio_filename = let finalize ~t0 = print_string "Finalizing..."; change_status Status.Stopped; - (* compress_block_file (); - *) send_log "status" 0 t0 "Done"; close_debug_socket (); ZMQ.Context.terminate zmq_context; diff --git a/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml index c8d5b13..ce02468 100644 --- a/ocaml/Qmcchem_result.ml +++ b/ocaml/Qmcchem_result.ml @@ -4,17 +4,6 @@ open Qptypes (** Display a table that can be plotted by gnuplot *) let display_table ~range property = let p = - (*---*) - if (property = "xxx") then - let open Random_variable in - (of_raw_data ~range Property.E_ref_zv) - /! (of_raw_data ~range Property.Psi_no_cusp_weight) - -! ( (of_raw_data ~range Property.E_loc_zv) - /! (of_raw_data ~range Property.Weight_zv) ) - /! (of_raw_data ~range Property.Psi_no_cusp_weight2) - +@ Input.Trial_wf_energy.to_float (Input.Trial_wf_energy. read () ) - else - (*---*) Property.of_string property |> Random_variable.of_raw_data ~range in @@ -157,35 +146,6 @@ let display_summary ~range = in List.iter properties ~f:print_property ; - (*---*) - let open Random_variable in - let p = (of_raw_data ~range Property.E_ref_zv) - /! (of_raw_data ~range Property.Psi_no_cusp_weight) - -! ( (of_raw_data ~range Property.E_loc_zv) - /! (of_raw_data ~range Property.Weight_zv) ) - /! (of_raw_data ~range Property.Psi_no_cusp_weight2) - +@ Input.Trial_wf_energy.to_float (Input.Trial_wf_energy. read () ) - in - Printf.printf "%20s : %s\n" - ("E_loc_zv(+)") - (Random_variable.to_string p); - (*---*) - -(*--- - let open Random_variable in - let e_trial = - Input.Trial_wf_energy.read () |> Input.Trial_wf_energy.to_float - in - let p = ( (of_raw_data ~range Property.E_pdmc) - /! (of_raw_data ~range Property.W_pdmc) ) - -! ( (of_raw_data ~range Property.E_loc) -@ e_trial ) - in - Printf.printf "%20s : %s %f\n" - ("E_pdmc(zv)") - (Random_variable.to_string p) - e_trial; ----*) - let cpu = Random_variable.of_raw_data ~range Property.Cpu |> Random_variable.sum diff --git a/ocaml/compile.sh b/ocaml/compile.sh index 8ae5c7c..9ac0bd8 100755 --- a/ocaml/compile.sh +++ b/ocaml/compile.sh @@ -25,7 +25,7 @@ if [[ ${MD5} != ${REF} ]] then echo ${MD5} > ${LSMD5_FILE} echo Finding dependencies in OCaml files - python ./ninja_ocaml.py + python2 ./ninja_ocaml.py fi ninja ${@} diff --git a/ocaml/ninja_ocaml.py b/ocaml/ninja_ocaml.py index 775b670..fde91e3 100755 --- a/ocaml/ninja_ocaml.py +++ b/ocaml/ninja_ocaml.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python2 # # Copyright 2015 Anthony Scemama # @@ -218,7 +218,7 @@ GENERATED_NINJA=generated.ninja # Name of the auto-generated ninja file rule create_generated - command = python ./ninja_ocaml.py + command = python2 ./ninja_ocaml.py description = Finding dependencies between modules rule run_ninja diff --git a/ocaml/qmcchem.ml b/ocaml/qmcchem.ml index ad341b3..59a1563 100644 --- a/ocaml/qmcchem.ml +++ b/ocaml/qmcchem.ml @@ -8,6 +8,7 @@ let command = "info" , Qmcchem_info.command ; "md5" , Qmcchem_md5.command ; "result", Qmcchem_result.command ; + "result-pdmc", Qmcchem_result_pdmc.command ; "run" , Qmcchem_run.command ; "stop" , Qmcchem_stop.command ; ] diff --git a/scripts/compile_irpf90.sh b/scripts/compile_irpf90.sh index 2b790d5..aceffce 100755 --- a/scripts/compile_irpf90.sh +++ b/scripts/compile_irpf90.sh @@ -33,7 +33,7 @@ then IRPF90_FLAGS="${IRPF90_FLAGS} --ninja" # Check IRPF90 version - if [[ $( ${IRPF90} -v | python -c "import sys ; print float(sys.stdin.read().rsplit('.',1)[0]) >= 1.6") == False ]] + if [[ $( ${IRPF90} -v | python2 -c "import sys ; print float(sys.stdin.read().rsplit('.',1)[0]) >= 1.6") == False ]] then echo "IRPF90 version >= 1.6 required" exit -1 diff --git a/scripts/compile_ocaml_dep.sh b/scripts/compile_ocaml_dep.sh index ce7f610..84fe181 100755 --- a/scripts/compile_ocaml_dep.sh +++ b/scripts/compile_ocaml_dep.sh @@ -25,7 +25,7 @@ if [[ ${MD5} != ${REF} ]] then echo ${MD5} > ${LSMD5_FILE} echo Finding dependencies in OCaml files - python ./ninja_ocaml.py || exit -1 + python2 ./ninja_ocaml.py || exit -1 fi diff --git a/scripts/create_properties_ezfio.py b/scripts/create_properties_ezfio.py index d2fda6c..ed82d8a 100755 --- a/scripts/create_properties_ezfio.py +++ b/scripts/create_properties_ezfio.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python2 # # Creates the properties.config file in the EZFIO directory. This is # done by reading all the properties written in the src/PROPERTIES diff --git a/scripts/create_properties_python.py b/scripts/create_properties_python.py index f763920..02e62f2 100755 --- a/scripts/create_properties_python.py +++ b/scripts/create_properties_python.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/env python2 import string import os diff --git a/src/PROPERTIES/properties.irp.f b/src/PROPERTIES/properties.irp.f index 93a4658..f50473f 100644 --- a/src/PROPERTIES/properties.irp.f +++ b/src/PROPERTIES/properties.irp.f @@ -1,4 +1,4 @@ -BEGIN_SHELL [ /usr/bin/env python ] +BEGIN_SHELL [ /usr/bin/env python2 ] import os from properties import properties root = os.environ['QMCCHEM_PATH'] @@ -59,7 +59,7 @@ END_SHELL ! DIMENSIONS !==========================================================================! -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * make_dims() END_SHELL diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f index 1eb082e..d4d8b80 100644 --- a/src/PROPERTIES/properties_energy.irp.f +++ b/src/PROPERTIES/properties_energy.irp.f @@ -172,14 +172,6 @@ BEGIN_PROVIDER [ double precision, dmc_zv_weight_half ] dmc_zv_weight_half = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ integer, E_loc_traj_size ] - implicit none - BEGIN_DOC - ! Size of E_loc_traj - END_DOC - E_loc_traj_size = pdmc_projection(pdmc_n_diag) -END_PROVIDER - !==========================================================================! ! PROPERTIES ! @@ -272,13 +264,10 @@ BEGIN_PROVIDER [ double precision, E_loc ] END_PROVIDER -BEGIN_PROVIDER [ double precision, psi_no_cusp_weight ] - psi_no_cusp_weight = 1.d0 +BEGIN_PROVIDER [ double precision, psi_nocusp_weight ] + psi_nocusp_weight = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, psi_no_cusp_weight2 ] - psi_no_cusp_weight2 = 1.d0 -END_PROVIDER BEGIN_PROVIDER [ double precision, E_loc_nocusp ] implicit none @@ -290,142 +279,59 @@ BEGIN_PROVIDER [ double precision, E_loc_nocusp ] psi_cusp_inv = psi_value_inv do_nucl_fitcusp = .False. TOUCH do_nucl_fitcusp - psi_no_cusp_weight = psi_value * psi_cusp_inv - psi_no_cusp_weight2 = psi_no_cusp_weight * psi_no_cusp_weight - E_loc_nocusp = psi_no_cusp_weight * E_loc + psi_nocusp_weight = psi_value * psi_cusp_inv + pdmc_norm_vmc_zv = psi_nocusp_weight + E_loc_nocusp = E_loc do_nucl_fitcusp = .True. + E_vmc_zv = E_loc_nocusp TOUCH do_nucl_fitcusp else E_loc_nocusp = E_loc - psi_no_cusp_weight = 1.d0 + E_vmc_zv = E_loc + pdmc_norm_vmc_zv = psi_nocusp_weight + psi_nocusp_weight = 1.d0 endif END_PROVIDER -BEGIN_PROVIDER [ double precision, E_ref_zv ] + +BEGIN_PROVIDER [ double precision, E_vmc_zv ] implicit none BEGIN_DOC - ! E_loc - E_trial + ! PDMC weight, computed in PDMC END_DOC - E_ref_zv = E_loc_nocusp - - E_ref_zv_min = min(E_ref_zv,E_ref_zv_min) - E_ref_zv_max = max(E_ref_zv,E_ref_zv_max) - SOFT_TOUCH E_ref_zv_min E_ref_zv_max - + E_vmc_zv = 0.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv ] +BEGIN_PROVIDER [ double precision, pdmc_norm_vmc_zv ] implicit none BEGIN_DOC - ! Estimate of the VMC energy in PDMC + ! PDMC norm END_DOC - E_loc_zv = E_loc_nocusp * psi_no_cusp_weight * dmc_zv_weight - - E_loc_zv_min = min(E_loc_zv,E_loc_zv_min) - E_loc_zv_max = max(E_loc_zv,E_loc_zv_max) - SOFT_TOUCH E_loc_zv_min E_loc_zv_max - + pdmc_norm_vmc_zv = 1.d0 END_PROVIDER - -BEGIN_PROVIDER [ double precision, Weight_zv ] +BEGIN_PROVIDER [ double precision, pdmc_w ] implicit none BEGIN_DOC - ! Average weight to go from PDMC to VMC + ! PDMC weight, computed in PDMC END_DOC - Weight_zv = dmc_zv_weight - + pdmc_w = 1.d0 END_PROVIDER - -BEGIN_PROVIDER [ double precision, E_loc_traj, (E_loc_traj_size) ] + BEGIN_PROVIDER [ double precision, pdmc_norm ] implicit none BEGIN_DOC - ! E_loc_traj + ! PDMC norm END_DOC - integer :: k,l - - do k=1,E_loc_traj_size - l = pdmc_projection_step(pdmc_n_diag)+k - if (l > pdmc_projection(pdmc_n_diag)) then - l = l - pdmc_projection(pdmc_n_diag) - endif -! E_loc_traj(k) = E_loc !_trajectory(l,pdmc_n_diag) -! E_loc_traj(k) += ( E_trial - E_loc_trajectory(l,pdmc_n_diag) ) *dmc_zv_weight - E_loc_traj(k) = E_loc_trajectory(l,pdmc_n_diag) - enddo - + pdmc_norm = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_zv_traj ] +BEGIN_PROVIDER [ double precision, pdmc_norm_nocusp ] implicit none BEGIN_DOC - ! Zero-variance parameter on E_loc + ! PDMC norm END_DOC - integer :: k - E_loc_zv_traj = 0.d0 - do k=1,pdmc_projection(pdmc_n_diag) - E_loc_zv_traj += E_loc !_trajectory(k,pdmc_n_diag) - E_loc_zv_traj += (E_trial - E_loc_trajectory(k,pdmc_n_diag)) * dmc_zv_weight - enddo - E_loc_zv_traj *= 1.d0/dble(pdmc_projection(pdmc_n_diag)) - - E_loc_zv_traj_min = min(E_loc_zv_traj,E_loc_zv_traj_min) - E_loc_zv_traj_max = max(E_loc_zv_traj,E_loc_zv_traj_max) - SOFT_TOUCH E_loc_zv_traj_min E_loc_zv_traj_max + pdmc_norm_nocusp = 1.d0 END_PROVIDER -BEGIN_PROVIDER [ double precision, E_pdmc ] - implicit none - BEGIN_DOC - ! PDMC local energy - END_DOC - - E_pdmc = E_loc* w_pdmc - -! E_pdmc = E_loc* w_pdmc - (E_loc - E_trial) - - -! double precision :: delta, w -! integer, save :: istep=1 -! double precision,save :: E0=0, Ezv -! istep = istep-1 -! if (istep == 0) then -! istep = 1000 -! E0 = E_loc -! Ezv = E_loc - E_trial -! endif -! E_pdmc = E0* w_pdmc - Ezv - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, w_pdmc ] - implicit none - BEGIN_DOC - ! PDMC local energy - END_DOC - double precision :: delta, w - integer, save :: istep=1 - istep = istep-1 - if (istep == 0) then - istep = 1000 - w_pdmc = 1.d0 - endif - - delta = E_loc - E_ref - if (delta >= 0.d0) then - w = dexp(-time_step * delta) - else - w = 2.d0 - dexp(time_step * delta) - endif - - w_pdmc = w * w_pdmc - -END_PROVIDER - - - - - - diff --git a/src/PROPERTIES/properties_general.irp.f b/src/PROPERTIES/properties_general.irp.f index c5342e6..b6c0ae9 100644 --- a/src/PROPERTIES/properties_general.irp.f +++ b/src/PROPERTIES/properties_general.irp.f @@ -47,22 +47,6 @@ BEGIN_PROVIDER [ double precision, wf_extension ] SOFT_TOUCH wf_extension_min wf_extension_max END_PROVIDER -BEGIN_PROVIDER [ double precision, pop_weight ] - implicit none - BEGIN_DOC - ! Weight of the SRMC population - END_DOC - include '../types.F' - if (qmc_method == t_SRMC) then - pop_weight = srmc_pop_weight_mult - else if (qmc_method == t_PDMC) then - pop_weight = pdmc_pop_weight_mult(pdmc_n_diag) - endif - pop_weight_min = min(pop_weight,pop_weight_min) - pop_weight_max = max(pop_weight,pop_weight_max) - SOFT_TOUCH pop_weight_min pop_weight_max - -END_PROVIDER BEGIN_PROVIDER [ double precision, drift_mod, (size_drift_mod) ] implicit none diff --git a/src/SAMPLING/block.irp.f b/src/SAMPLING/block.irp.f index 2492e12..cc417b5 100644 --- a/src/SAMPLING/block.irp.f +++ b/src/SAMPLING/block.irp.f @@ -1,6 +1,6 @@ ! Providers of *_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ diff --git a/src/SAMPLING/dmc_step.irp.f b/src/SAMPLING/dmc_step.irp.f index 27397bc..259bd41 100644 --- a/src/SAMPLING/dmc_step.irp.f +++ b/src/SAMPLING/dmc_step.irp.f @@ -1,6 +1,6 @@ ! Providers of *_dmc_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ @@ -67,7 +67,7 @@ END_SHELL integer :: k, i_walk, i_step -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -150,7 +150,7 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -265,7 +265,7 @@ END_SHELL factor = 1.d0/block_weight SOFT_TOUCH block_weight -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then diff --git a/src/SAMPLING/fkmc_step.irp.f b/src/SAMPLING/fkmc_step.irp.f index c7e0e95..43efcf7 100644 --- a/src/SAMPLING/fkmc_step.irp.f +++ b/src/SAMPLING/fkmc_step.irp.f @@ -1,6 +1,6 @@ ! Providers of *_fkmc_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ @@ -77,7 +77,7 @@ END_SHELL integer :: k, i_walk, i_step -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -173,7 +173,7 @@ END_SHELL psi_value_save(i_walk) = psi_value E_loc_save(i_walk) = E_loc -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -339,7 +339,7 @@ END_SHELL factor = 1.d0/block_weight SOFT_TOUCH block_weight -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f index 4815cff..b42407c 100644 --- a/src/SAMPLING/pdmc_step.irp.f +++ b/src/SAMPLING/pdmc_step.irp.f @@ -1,6 +1,6 @@ ! Providers of *_pdmc_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ @@ -8,400 +8,369 @@ t = """ &BEGIN_PROVIDER [ $T, $X_pdmc_block_walk_kahan $D2 ] &BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk $D1 ] &BEGIN_PROVIDER [ $T, $X_2_pdmc_block_walk_kahan $D2 ] +!&BEGIN_PROVIDER [ $T, $X_vmc_pdmc_block_walk $D1 ] +!&BEGIN_PROVIDER [ $T, $X_vmc_pdmc_block_walk_kahan $D2 ] +!&BEGIN_PROVIDER [ $T, $X_2_vmc_pdmc_block_walk $D1 ] +!&BEGIN_PROVIDER [ $T, $X_2_vmc_pdmc_block_walk_kahan $D2 ] implicit none BEGIN_DOC -! pdMC averages of $X. Computed in E_loc_pdmc_block_walk + ! PDMC averages of $X. Computed in E_loc_pdmc_block_walk END_DOC $X_pdmc_block_walk = 0.d0 $X_pdmc_block_walk_kahan = 0.d0 $X_2_pdmc_block_walk = 0.d0 $X_2_pdmc_block_walk_kahan = 0.d0 +! $X_vmc_pdmc_block_walk = 0.d0 +! $X_vmc_pdmc_block_walk_kahan = 0.d0 +! $X_2_vmc_pdmc_block_walk = 0.d0 +! $X_2_vmc_pdmc_block_walk_kahan = 0.d0 END_PROVIDER + +BEGIN_PROVIDER [ $T, $X_pdmc_trajectory $D3 ] + implicit none + BEGIN_DOC + ! Value of $X along the PDMC trajectory + END_DOC + $X_pdmc_trajectory = 1.d0 +END_PROVIDER + """ for p in properties: if p[1] != 'e_loc': if p[2] == "": D1 = "" D2 = ", (3)" + D3 = ", (0:pdmc_n_projection_steps)" else: D1 = ", ("+p[2][1:-1]+")" D2 = ", ("+p[2][1:-1]+",3)" - print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2) + D3 = ", ("+p[2][1:-1]+",0:pdmc_n_projection_steps)" + print t.replace("$X",p[1]).replace("$T",p[0]).replace("$D1",D1).replace("$D2",D2).replace("$D3",D3) END_SHELL +BEGIN_PROVIDER [ double precision, E_loc_pdmc_trajectory, (0:pdmc_n_projection_steps) ] + E_loc_pdmc_trajectory = -huge(1.d0) +END_PROVIDER - BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk ] -&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk ] -&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan, (3) ] -&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ] - implicit none - include '../types.F' - BEGIN_DOC -! Properties averaged over the block using the PDMC method - END_DOC - - real, allocatable :: elec_coord_tmp(:,:,:) - integer :: mod_align - double precision :: E_loc_save(4,walk_num_dmc_max) - double precision :: psi_value_save(walk_num) - double precision :: psi_value_save_tmp(walk_num) - double precision :: pdmc_weight(walk_num) - double precision :: block_weight_vmc - double precision, allocatable :: psi_grad_psi_inv_save(:,:,:) - double precision, allocatable :: psi_grad_psi_inv_save_tmp(:,:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_grad_psi_inv_save_tmp + BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk ] +&BEGIN_PROVIDER [ double precision, E_loc_pdmc_block_walk_kahan , (3) ] +&BEGIN_PROVIDER [ double precision, E_loc_2_pdmc_block_walk_kahan, (3) ] +!&BEGIN_PROVIDER [ double precision, E_loc_vmc_pdmc_block_walk ] +!&BEGIN_PROVIDER [ double precision, E_loc_2_vmc_pdmc_block_walk ] +!&BEGIN_PROVIDER [ double precision, E_loc_vmc_pdmc_block_walk_kahan , (3) ] +!&BEGIN_PROVIDER [ double precision, E_loc_2_vmc_pdmc_block_walk_kahan, (3) ] + implicit none + include '../types.F' + BEGIN_DOC + ! Properties averaged over the block using the PDMC method + END_DOC + + integer :: i,j,l + double precision :: E_loc_save(4) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: E_loc_save - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: psi_value_save_tmp - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: pdmc_weight - allocate ( psi_grad_psi_inv_save(elec_num_8,3,walk_num) , & - psi_grad_psi_inv_save_tmp(elec_num_8,3,walk_num) , & - elec_coord_tmp(mod_align(elec_num+1),3,walk_num) ) - psi_value_save = 0.d0 - psi_value_save_tmp = 0.d0 - pdmc_weight = 1.d0 -! Initialization - if (vmc_algo /= t_Brownian) then - call abrt(irp_here,'PDMC should run with Brownian algorithm') - endif + calc_pdmc_norm = .True. + calc_pdmc_norm_nocusp = .True. - integer :: k, i_walk, i_step + PROVIDE time_step -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ - if (calc_$X) then - !DIR$ VECTOR ALIGNED - $X_pdmc_block_walk = 0.d0 - !DIR$ VECTOR ALIGNED - $X_pdmc_block_walk_kahan = 0.d0 - !DIR$ VECTOR ALIGNED - $X_2_pdmc_block_walk = 0.d0 - !DIR$ VECTOR ALIGNED - $X_2_pdmc_block_walk_kahan = 0.d0 - endif + if (calc_$X) then + $X_pdmc_block_walk = 0.d0 + $X_pdmc_block_walk_kahan = 0.d0 + $X_2_pdmc_block_walk = 0.d0 + $X_2_pdmc_block_walk_kahan = 0.d0 +! $X_vmc_pdmc_block_walk = 0.d0 +! $X_vmc_pdmc_block_walk_kahan = 0.d0 +! $X_2_vmc_pdmc_block_walk = 0.d0 +! $X_2_vmc_pdmc_block_walk_kahan = 0.d0 + $X_min = huge(1.) + $X_max =-huge(1.) + endif """ for p in properties: print t.replace("$X",p[1]) END_SHELL - logical :: loop - integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max + block_weight = 0.d0 - loop = .True. - call system_clock(cpu0, count_rate, count_max) - cpu2 = cpu0 + ! CPU time at the beginning of the block + ! -------------------------------------- - block_weight = 0.d0 - block_weight_vmc = 0.d0 + integer*8 :: cpu0, cpu1, cpu2, count_rate, count_max + call system_clock(cpu0, count_rate, count_max) + cpu2 = cpu0 - real, external :: accep_rate - double precision :: delta, thr - thr = 2.d0/time_step_sq - logical :: first_loop - first_loop = .True. - if (walk_num > 1) then - call abrt(irp_here,'walk_num > 1') - endif + integer :: istep, istep_t_2 + logical :: loop, first_loop + first_loop = .True. + loop = .True. - 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) -! 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 - do i=1,elec_num+1 - elec_coord(i,l) = elec_coord_full(i,l,i_walk) - enddo - do i=1,elec_num - psi_grad_psi_inv_x(i) = psi_grad_psi_inv_save(i,1,i_walk) - psi_grad_psi_inv_y(i) = psi_grad_psi_inv_save(i,2,i_walk) - psi_grad_psi_inv_z(i) = psi_grad_psi_inv_save(i,3,i_walk) - enddo - psi_value = psi_value_save(i_walk) - E_loc = E_loc_save(1,i_walk) - enddo - SOFT_TOUCH elec_coord psi_grad_psi_inv_x psi_grad_psi_inv_y psi_grad_psi_inv_z psi_value E_loc - else - do l=1,3 - do i=1,elec_num+1 - elec_coord(i,l) = elec_coord_full(i,l,i_walk) - enddo - enddo - TOUCH elec_coord - psi_value_save(i_walk) = psi_value - E_loc_save(:,i_walk) = E_loc - endif - - double precision :: p,q - real :: delta_x - logical :: accepted - call brownian_step(p,q,accepted,delta_x) - -! if ( psi_value * psi_value_save(i_walk) >= 0.d0 ) then - -!2 delta = (E_loc+E_loc_save(1,i_walk))*0.5d0 -!3 delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1,i_walk) - E_loc_save(2,i_walk))/12.d0 -! delta = (9.d0*E_loc+19.d0*E_loc_save(1,i_walk)-5.d0*E_loc_save(2,i_walk)+E_loc_save(3,i_walk))*(1.d0/24.d0) - delta = -((-251.d0*E_loc)-646.d0*E_loc_save(1,i_walk)+264.d0*E_loc_save(2,i_walk)-& - 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 - pdmc_weight(i_walk) = 2.d0-dexp(dtime_step*delta) - endif - - elec_coord(elec_num+1,1) += p*time_step - elec_coord(elec_num+1,2) = E_loc - elec_coord(elec_num+1,3) = pdmc_weight(i_walk) * pdmc_pop_weight_mult(pdmc_n_diag) - do l=1,3 - do i=1,elec_num+1 - elec_coord_full(i,l,i_walk) = elec_coord(i,l) - enddo - enddo - do i=1,elec_num - psi_grad_psi_inv_save(i,1,i_walk) = psi_grad_psi_inv_x(i) - psi_grad_psi_inv_save(i,2,i_walk) = psi_grad_psi_inv_y(i) - psi_grad_psi_inv_save(i,3,i_walk) = psi_grad_psi_inv_z(i) - enddo - - psi_value_save(i_walk) = psi_value - if (accepted) then - E_loc_save(4,i_walk) = E_loc_save(3,i_walk) - E_loc_save(3,i_walk) = E_loc_save(2,i_walk) - E_loc_save(2,i_walk) = E_loc_save(1,i_walk) - E_loc_save(1,i_walk) = E_loc - endif - if (calc_e_loc_zv) then - 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)) - else - dmc_zv_weight = 0.d0 - dmc_zv_weight_half = 0.d0 - endif - TOUCH dmc_zv_weight dmc_zv_weight_half E_loc_trajectory - endif - do k=1,pdmc_n_diag - E_loc_trajectory(pdmc_projection_step(k),k) = E_loc + do l=1,3 + do i=1,elec_num+1 + elec_coord(i,l) = elec_coord_full(i,l,1) enddo - -! do i=1,pdmc_n_diag+1 -! E_loc_zv(i) = E_loc * pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight + (E_trial-E_loc) * dmc_zv_weight -! E_loc_zv(i+pdmc_n_diag+1) = pdmc_pop_weight_mult(i-1) * pdmc_weight(i_walk) * dmc_zv_weight -! enddo + enddo + TOUCH elec_coord + E_loc_save = E_loc -BEGIN_SHELL [ /usr/bin/python ] + do while (loop) + +! istep_t_2 = pdmc_n_projection_steps/2 + istep_t_2 = -1 +! istep_t_2 = 0 + + do istep = 0, pdmc_n_projection_steps + + istep_t_2 = istep_t_2+1 + if (istep_t_2 > pdmc_n_projection_steps) then + istep_t_2 = 0 + endif + + ! Brownian step + ! ------------- + + double precision :: p,q + real :: delta_x + logical :: accepted + + E_loc_save(4) = E_loc_save(3) + E_loc_save(3) = E_loc_save(2) + E_loc_save(2) = E_loc_save(1) + E_loc_save(1) = E_loc + call brownian_step(p,q,accepted,delta_x) + elec_coord(elec_num+1,1) += p*time_step + elec_coord(elec_num+1,2) = E_loc + elec_coord(elec_num+1,3) = 1. + + psi_nocusp_weight_pdmc_trajectory(istep) = psi_nocusp_weight + + ! PDMC weight + ! ----------- + + double precision :: delta, w0, wt, w, w_prod + +! if ( psi_value * psi_value_save >= 0.d0 ) then + +!2 delta = (E_loc+E_loc_save(1))*0.5d0 +!3 delta = (5.d0 * E_loc + 8.d0 * E_loc_save(1) - E_loc_save(2))/12.d0 + delta = (9.d0*E_loc+19.d0*E_loc_save(1)-5.d0*E_loc_save(2)+E_loc_save(3))/24.d0 +!5 delta = -((-251.d0*E_loc)-646.d0*E_loc_save(1)+264.d0*E_loc_save(2)-& +!5 106.d0*E_loc_save(3)+19.d0*E_loc_save(4))/720.d0 + + delta = (delta - E_ref)*p + if (delta >= 0.d0) then + pdmc_w = dexp(-dtime_step*delta) + else + pdmc_w = 2.d0 - dexp(dtime_step*delta) + endif + pdmc_w_pdmc_trajectory(istep) = 1.d0 + w_prod = product(pdmc_w_pdmc_trajectory) + pdmc_w_pdmc_trajectory(istep) = pdmc_w + + ! Observables + ! ----------- + + double precision :: x, x2 + integer :: jstep + block_weight += 1.d0 + +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ - if (calc_$X) then - ! Kahan's summation algorithm to compute these sums reducing the rounding error: - ! $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_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) - endif + + if (calc_$X) then + + if ($nocusp) then + wt = psi_nocusp_weight_pdmc_trajectory(istep) + if (istep /= pdmc_n_projection_steps) then + jstep = istep + 1 + else + jstep = 0 + endif + w0 = psi_nocusp_weight_pdmc_trajectory(jstep) + else + w0 = 1.d0 + wt = 1.d0 + endif + + $X_pdmc_trajectory($D3 istep) = $X + + + ! Kahan's summation algorithm to compute these sums reducing the rounding error: + ! $X_pdmc_block_walk $D1 += x + ! $X_2_pdmc_block_walk $D1 += x_2 + ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm + + ! PDMC observables + + $X_pdmc_block_walk_kahan($D2 3) = $X_pdmc_trajectory($D3 istep_t_2) * w0 * wt * w_prod - $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_2_pdmc_block_walk_kahan($D2 3) = ($X_pdmc_trajectory($D3 istep_t_2) * w0 * wt * w_prod)**2 - $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) + + +! ! VMC observables +! +! $X_vmc_pdmc_block_walk_kahan($D2 3) = $X_pdmc_trajectory($D3 istep_t_2) * wt * wt - $X_vmc_pdmc_block_walk_kahan($D2 1) +! $X_vmc_pdmc_block_walk_kahan($D2 2) = $X_vmc_pdmc_block_walk $D1 + $X_vmc_pdmc_block_walk_kahan($D2 3) +! $X_vmc_pdmc_block_walk_kahan($D2 1) = ($X_vmc_pdmc_block_walk_kahan($D2 2) - $X_vmc_pdmc_block_walk $D1 ) & +! - $X_vmc_pdmc_block_walk_kahan($D2 3) +! $X_vmc_pdmc_block_walk $D1 = $X_vmc_pdmc_block_walk_kahan($D2 2) +! +! +! $X_2_vmc_pdmc_block_walk_kahan($D2 3) = $X_pdmc_trajectory($D3 istep_t_2)**2 * wt * wt - $X_2_vmc_pdmc_block_walk_kahan($D2 1) +! $X_2_vmc_pdmc_block_walk_kahan($D2 2) = $X_2_vmc_pdmc_block_walk $D1 + $X_2_vmc_pdmc_block_walk_kahan($D2 3) +! $X_2_vmc_pdmc_block_walk_kahan($D2 1) = ($X_2_vmc_pdmc_block_walk_kahan($D2 2) - $X_2_vmc_pdmc_block_walk $D1 ) & +! - $X_2_vmc_pdmc_block_walk_kahan($D2 3) +! $X_2_vmc_pdmc_block_walk $D1 = $X_2_vmc_pdmc_block_walk_kahan($D2 2) + endif """ for p in properties: - if p[2] == "": - D1 = "" - D2 = "" - else: - D1 = "("+":"*(p[2].count(',')+1)+")" - D2 = ":"*(p[2].count(',')+1)+"," - print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2) + if "vmc_zv" not in p[1]: + if p[2] == "": + D1 = "" + D2 = "" + D3 = "" + else: + D1 = "("+":"*(p[2].count(',')+1)+")" + D2 = ":"*(p[2].count(',')+1)+"," + D3 = ":"*(p[2].count(',')+1)+"," + if 'nocusp' in p[1]: + nocusp = '.True.' + else: + nocusp = '.False.' + print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2).replace("$D3",D3).replace("$nocusp",nocusp) END_SHELL - block_weight += pdmc_pop_weight_mult(pdmc_n_diag) * pdmc_weight(i_walk) - block_weight_vmc += 1.d0 +BEGIN_SHELL [ /usr/bin/env python2 ] +from properties import * +t = """ - pdmc_pop_weight_mult(0) = 1.d0/pdmc_weight(i_walk) -! do k=0,pdmc_n_diag/2 -! do l=0,pdmc_n_diag/2 -! H(k,l) += E_loc*pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk) -! S(k,l) += pdmc_pop_weight_mult(k+l) * pdmc_weight(i_walk) -! enddo -! enddo -! H = H + (E_trial - E_loc) + if (calc_$X) then -! else -! pdmc_weight(i_walk) = 1.d0 -! pdmc_pop_weight(:,:) = 1.d0 -! pdmc_pop_weight_mult(:) = 1.d0 -! endif - - + if ($nocusp) then + wt = psi_nocusp_weight_pdmc_trajectory(istep_t_2) + else + wt = 1.d0 + endif - do k=1,pdmc_n_diag - ! Move to the next projection step - if (pdmc_projection(pdmc_n_diag) > 0) then - pdmc_projection_step(k) = mod(pdmc_projection_step(k),pdmc_projection(k))+1 - else - pdmc_projection_step(k) = 1 - endif + $X_pdmc_trajectory($D3 istep) = $X + - ! Eventually, recompute the weight of the population - if (pdmc_projection_step(k) == 1) then - pdmc_pop_weight_mult(k) = 1.d0 - do l=1,pdmc_projection(k) - pdmc_pop_weight_mult(k) *= pdmc_pop_weight(l,k) - enddo - endif + ! Kahan's summation algorithm to compute these sums reducing the rounding error: + ! $X_pdmc_block_walk $D1 += x + ! $X_2_pdmc_block_walk $D1 += x_2 + ! see http://en.wikipedia.org/wiki/Kahan_summation_algorithm - ! 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) + ! PDMC observables - ! Update the running population weight - pdmc_pop_weight_mult(k) *= pdmc_pop_weight(pdmc_projection_step(k),k) + $X_pdmc_block_walk_kahan($D2 3) = $X_pdmc_trajectory($D3 istep_t_2) * wt * wt - $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_2_pdmc_block_walk_kahan($D2 3) = ($X_pdmc_trajectory($D3 istep_t_2) * wt * wt )**2 - $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) + + endif +""" +for p in properties: + if "vmc_zv" in p[1]: + if p[2] == "": + D1 = "" + D2 = "" + D3 = "" + else: + D1 = "("+":"*(p[2].count(',')+1)+")" + D2 = ":"*(p[2].count(',')+1)+"," + D3 = ":"*(p[2].count(',')+1)+"," + if 'nocusp' in p[1]: + nocusp = '.True.' + else: + nocusp = '.False.' + print t.replace("$X",p[1]).replace("$D1",D1).replace("$D2",D2).replace("$D3",D3).replace("$nocusp",nocusp) + +END_SHELL + + end do ! istep + + + call system_clock(cpu1, count_rate, count_max) + if (cpu1 < cpu0) then + cpu1 = cpu1+cpu0 + endif + loop = dble(cpu1-cpu0)/dble(count_rate) < block_time + if (cpu1-cpu2 > count_rate) then + integer :: do_run + call get_running(do_run) + loop = loop.and.(do_run == t_Running) + cpu2 = cpu1 + endif + + enddo ! while (loop) + + do l=1,3 + do i=1,elec_num+1 + elec_coord_full(i,l,1) = elec_coord(i,l) + enddo enddo + + double precision :: factor + factor = 1.d0/block_weight - call system_clock(cpu1, count_rate, count_max) - if (cpu1 < cpu0) then - cpu1 = cpu1+cpu0 - endif - loop = dble(cpu1-cpu0)/dble(count_rate) < block_time - if (cpu1-cpu2 > count_rate) then - integer :: do_run - call get_running(do_run) - loop = loop.and.(do_run == t_Running) - cpu2 = cpu1 - endif - - 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 - -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ - if (calc_$X) then - $X_pdmc_block_walk *= factor - $X_2_pdmc_block_walk *= factor - endif + if (calc_$X) then + $X_pdmc_block_walk *= factor + $X_2_pdmc_block_walk *= factor +! $X_vmc_pdmc_block_walk *= factor +! $X_2_vmc_pdmc_block_walk *= factor + endif """ for p in properties: - print t.replace("$X",p[1]) +# if p[1] in [ "pdmc_norm", "e_loc"]: + print t.replace("$X",p[1]) END_SHELL -!if (calc_e_loc_zv) then -! factor = block_weight/block_weight_vmc -! E_loc_zv_pdmc_block_walk *= factor -! E_loc_zv_pdmc_block_walk = E_trial + E_loc_pdmc_block_walk - E_loc_zv_pdmc_block_walk -! E_loc_zv_2_pdmc_block_walk *= factor -! E_loc_zv_2_pdmc_block_walk = E_trial*E_trial*block_weight_vmc/block_weight ! Approximation -!endif - - -! H(0,0) = H(3,3) -! H(1,0) = H(4,3) -! H(0,1) = H(3,4) -! H(1,1) = H(4,4) -! S(0,0) = S(3,3) -! S(1,0) = S(4,3) -! S(0,1) = S(3,4) -! S(1,1) = S(4,4) -! -! print *, H(0,0)/S(0,0) -! print *, H(1,1)/S(1,1) -! print *, '' -! -! call dsygv(1, 'N', 'U', pdmc_n_diag/2+1, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info) -! call dsygv(1, 'N', 'U', 2, H, pdmc_n_diag/2+1, S, pdmc_n_diag/2+1, w, work, 3*(pdmc_n_diag+1), info) -! E_loc_zv_diag_pdmc_block_walk = w(0) -! print *, w - - deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) + SOFT_TOUCH elec_coord_full block_weight END_PROVIDER -BEGIN_PROVIDER [ double precision, E_loc_trajectory, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ] - implicit none - BEGIN_DOC -! Local energies along the trajectory - END_DOC - E_loc_trajectory(:,:) = E_trial + + +BEGIN_PROVIDER [ integer, pdmc_n_projection_steps ] + implicit none + BEGIN_DOC + ! Number of projection steps for PDMC + END_DOC + real :: pdmc_projection_time + pdmc_projection_time = 1. + call get_simulation_srmc_projection_time(pdmc_projection_time) + pdmc_n_projection_steps = int( pdmc_projection_time/time_step) + 1 END_PROVIDER - - BEGIN_PROVIDER [ integer, pdmc_projection, (pdmc_n_diag) ] -&BEGIN_PROVIDER [ integer, pdmc_projection_step, (pdmc_n_diag) ] - implicit none - BEGIN_DOC -! Number of projection steps for PDMC - END_DOC - real :: pdmc_projection_time - pdmc_projection_time = 1. - call get_simulation_srmc_projection_time(pdmc_projection_time) - pdmc_projection(pdmc_n_diag) = int( pdmc_projection_time/time_step) - integer :: k - do k=1,pdmc_n_diag-1 - pdmc_projection(k) = k*pdmc_projection(pdmc_n_diag)/pdmc_n_diag - enddo - pdmc_projection_step(:) = 0 -END_PROVIDER - -BEGIN_PROVIDER [ double precision, pdmc_pop_weight, (0:pdmc_projection(pdmc_n_diag)+1,pdmc_n_diag) ] - implicit none - BEGIN_DOC -! Population weight of PDMC - END_DOC - pdmc_pop_weight(:,:) = 1.d0 -END_PROVIDER - -BEGIN_PROVIDER [ double precision, pdmc_pop_weight_mult, (0:pdmc_n_diag) ] - implicit none - BEGIN_DOC -! Population weight of PDMC - END_DOC - pdmc_pop_weight_mult(:) = 1.d0 -END_PROVIDER - - -BEGIN_PROVIDER [ integer, pdmc_n_diag ] - implicit none - BEGIN_DOC -! Size of the matrix to diagonalize - END_DOC - pdmc_n_diag = 1 -END_PROVIDER - - - diff --git a/src/SAMPLING/srmc_step.irp.f b/src/SAMPLING/srmc_step.irp.f index d4066c2..215ce0a 100644 --- a/src/SAMPLING/srmc_step.irp.f +++ b/src/SAMPLING/srmc_step.irp.f @@ -1,6 +1,6 @@ ! Providers of *_srmc_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ @@ -73,7 +73,7 @@ END_SHELL integer :: k, i_walk, i_step -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -210,17 +210,9 @@ END_SHELL E_loc_save(2,i_walk) = E_loc_save(1,i_walk) E_loc_save(1,i_walk) = E_loc endif - if (calc_e_loc_zv) then - if (dabs(srmc_weight(i_walk)*srmc_pop_weight_mult) > 1.d-15) then - dmc_zv_weight = 1.d0/(srmc_weight(i_walk)*srmc_pop_weight_mult) - else - dmc_zv_weight = 0.d0 - endif - TOUCH dmc_zv_weight - endif -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -353,7 +345,7 @@ END_SHELL factor = 1.d0/block_weight SOFT_TOUCH block_weight -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -365,13 +357,6 @@ for p in properties: print t.replace("$X",p[1]) END_SHELL - if (calc_e_loc_zv) then - factor = block_weight/block_weight_vmc - E_loc_zv_srmc_block_walk *= factor - E_loc_zv_srmc_block_walk = E_trial + E_loc_srmc_block_walk - E_loc_zv_srmc_block_walk - E_loc_zv_2_srmc_block_walk *= factor - E_loc_zv_2_srmc_block_walk = E_trial*E_trial*block_weight_vmc/block_weight ! Approximation - endif deallocate ( elec_coord_tmp, psi_grad_psi_inv_save, psi_grad_psi_inv_save_tmp ) END_PROVIDER diff --git a/src/SAMPLING/vmc_step.irp.f b/src/SAMPLING/vmc_step.irp.f index e09ee80..16ee78c 100644 --- a/src/SAMPLING/vmc_step.irp.f +++ b/src/SAMPLING/vmc_step.irp.f @@ -1,6 +1,6 @@ ! Providers of *_vmc_block_walk !============================== -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ @@ -39,7 +39,7 @@ END_SHELL PROVIDE time_step -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -93,7 +93,7 @@ END_SHELL block_weight += 1.d0 -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -152,8 +152,7 @@ END_SHELL double precision :: factor factor = 1.d0/block_weight - SOFT_TOUCH block_weight -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * t = """ if (calc_$X) then @@ -165,7 +164,7 @@ for p in properties: print t.replace("$X",p[1]) END_SHELL - SOFT_TOUCH elec_coord_full + SOFT_TOUCH elec_coord_full block_weight END_PROVIDER diff --git a/src/TOOLS/random.irp.f b/src/TOOLS/random.irp.f index 65bed3e..fe99635 100644 --- a/src/TOOLS/random.irp.f +++ b/src/TOOLS/random.irp.f @@ -54,7 +54,7 @@ BEGIN_PROVIDER [ integer*8, seed, (5) ] double precision :: r integer*8 :: pid8 read(current_PID,*) pid8 - pid8 = iand( ishft(pid8, 32), pid8) + pid8 = iand( ishft(pid8, 31), pid8) do i=1,12 clock(i) = i enddo diff --git a/src/ZMQ/qmc.irp.f b/src/ZMQ/qmc.irp.f index d2190cd..8e58cbb 100644 --- a/src/ZMQ/qmc.irp.f +++ b/src/ZMQ/qmc.irp.f @@ -101,7 +101,7 @@ subroutine run_qmc(cpu0) call zmq_send_real(msg,elec_coord_full,size(elec_coord_full)) endif -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] from properties import * derivlist = [] diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 515df24..950b66b 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -1,4 +1,4 @@ -BEGIN_SHELL [ /usr/bin/python ] +BEGIN_SHELL [ /usr/bin/env python2 ] data = [ \ ("electrons_elec_coord_pool_size" , "integer" , "" ), diff --git a/src/properties.py b/src/properties.py index ff180f2..0f02810 100755 --- a/src/properties.py +++ b/src/properties.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/env python2 import string import os