LCPQ/qmcchem
10
1
Fork 1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-07-26 20:57:26 +02:00
Code Releases Activity

Merge branch 'master' of gitlab.com:scemama/qmcchem

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-11 23:18:24 +01:00
commit 7aa1ee504d
27 changed files with 1768 additions and 373 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -8,4 +8,4 @@ make.config
src/ZMQ/f77_zmq.h
ezfio_config/properties.config
bin/
*.mod

6
.gitmodules vendored Normal file
View File

@ -0,0 +1,6 @@
[submodule "EZFIO"]
path = EZFIO
url = https://gitlab.com/scemama/EZFIO.git
[submodule "f77_zmq"]
path = f77_zmq
url = https://github.com/zeromq/f77_zmq.git

1
EZFIO Submodule

@ -0,0 +1 @@
Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c

11
Makefile
View File

@ -1,3 +1,14 @@
# QMC=Chem
FC = gfortran
LIBS =
CPPFLAGS =
LDFLAGS=
ZMQ_LIBS = -lzmq
ZMQ_LDFLAGS =
ZMQ_CPPFLAGS =
include make.config

76
Makefile.am Normal file
View File

@ -0,0 +1,76 @@
# QMC=Chem
FC = @FC@
LIBS = @LIBS@
CPPFLAGS = @CPPFLAGS@
LDFLAGS= @LDFLAGS@
ZMQ_LIBS = @ZMQ_LIBS@
ZMQ_LDFLAGS = @ZMQ_LDFLAGS@
ZMQ_CPPFLAGS = @ZMQ_CPPFLAGS@
default_target: all
.phony: clean always
clean:
./scripts/clean.sh
# put all files of PROPERTIES here
# --------------------------------
ezfio_config/properties.config ocaml/Property.ml: scripts/create_properties_python.py src/PROPERTIES/properties.irp.f src/PROPERTIES/properties_energy.irp.f src/PROPERTIES/properties_general.irp.f src/irpf90.make
bash -c "source qmcchemrc ; exec python2 ./scripts/create_properties_ezfio.py"
# EZFIO
# -----
build lib/libezfio.a lib/libezfio_irp.a EZFIO/lib/libezfio.a EZFIO/lib/libezfio_irp.a EZFIO/Ocaml/ezfio.ml EZFIO/Python/ezfio.py: ezfio_config/qmc.config ezfio_config/properties.config make.config scripts/create_properties_ezfio.py src/tags src/irpf90_entities src/irpf90.make src/irpf90.make
./scripts/compile_ezfio.sh
# Fortran executables
# -------------------
always: /dev/null
src/tags src/irpf90_entities src/irpf90.make: make.config always
./scripts/compile_irpf90.sh
src/MAIN/qmc src/MAIN/vmc_test src/MAIN/qmc_create_walkers src/MAIN/qmcchem_info: lib/libezfio.a lib/libezfio_irp.a src/tags src/irpf90_entities src/irpf90.make src/irpf90.make
./scripts/compile_src.sh
# OCaml
# -----
ocaml/qmcchem: EZFIO/Ocaml/ezfio.ml $(wildcard ocaml/*.ml)
./scripts/compile_ocaml.sh
# Archive
# -------
qmcchem.tar.gz: all
git archive --format tar.gz HEAD --prefix "QmcChem/" -7 -o qmcchem.tar.gz
# Binaries
# --------
bin/qmc: src/MAIN/qmc
cp src/MAIN/qmc bin/
bin/qmcchem_info: src/MAIN/qmcchem_info
cp src/MAIN/qmcchem_info bin/
bin/qmc_create_walkers: src/MAIN/qmc_create_walkers
cp src/MAIN/qmc_create_walkers bin/
bin/qmcchem: ocaml/qmcchem
cp ocaml/qmcchem bin/
all: bin/qmc bin/qmcchem_info bin/qmc_create_walkers bin/qmcchem

57
configure.ac Normal file
View File

@ -0,0 +1,57 @@
# -*- Autoconf -*-
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.69])
AC_INIT([QMC=Chem], [2.0.0], [https://gitlab.com/scemama/qmcchem/-/issues/new])
AM_INIT_AUTOMAKE([foreign subdir-objects silent-rules])
AC_CONFIG_SRCDIR([README.md])
AC_CONFIG_MACRO_DIR([m4])
AC_CONFIG_HEADERS([include/config.h])
# Checks for programs.
AC_PROG_AWK
AC_PROG_CC
AC_PROG_INSTALL
AC_PROG_LN_S
AC_PROG_MAKE_SET
AX_PROG_IRPF90
# Fortran compiler checks
AC_PROG_FC([ifort gfortran])
AC_FC_LINE_LENGTH([512])
# Checks for libraries.
AX_ZMQ([4.0], [], [ AC_MSG_ERROR([Please install ZeroMQ >= 4.0]) ])
AX_BLAS()
AX_LAPACK()
# Required by EZFIO
AC_CHECK_FUNCS([mkdir strerror])
AC_TYPE_SIZE_T
# Required by OCaml C bindings
AC_CHECK_FUNCS([inet_ntoa])
AC_CHECK_FUNCS([memset])
AC_CHECK_FUNCS([socket])
AC_CHECK_HEADERS([arpa/inet.h netinet/in.h sys/ioctl.h sys/socket.h])
# Required by ZeroMQ
AC_CHECK_HEADERS([stddef.h])
AC_TYPE_INT32_T
AC_TYPE_UINT16_T
AC_TYPE_UINT32_T
AC_TYPE_UINT8_T
# Required by QMCkl
#AC_CHECK_HEADER_STDBOOL
#AC_TYPE_INT64_T
#AC_TYPE_UINT64_T
AC_LANG([Fortran])
AC_CONFIG_FILES([Makefile])
AC_OUTPUT

1
configure.sh
View File

@ -41,6 +41,7 @@ export LD_LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LD_LIBRARY_PATH}"
export LIBRARY_PATH="\${QMCCHEM_PATH}/lib:\${LIBRARY_PATH}"
export QMCCHEM_MPIRUN="mpirun"
export QMCCHEM_MPIRUN_FLAGS=""
#export QMCCHEM_IO="B"
export C_INCLUDE_PATH="\${QMCCHEM_PATH}/include:\${C_INCLUDE_PATH}"
#export QMCCHEM_NIC=ib0
source \${QMCCHEM_PATH}/irpf90/bin/irpman

1
f77_zmq Submodule

@ -0,0 +1 @@
Subproject commit 934e063881553e42d81c5d1aaed35822988529f0

241
m4/ax_blas.m4 Normal file
View File

@ -0,0 +1,241 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_blas.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_BLAS([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]])
#
# DESCRIPTION
#
# This macro looks for a library that implements the BLAS linear-algebra
# interface (see http://www.netlib.org/blas/). On success, it sets the
# BLAS_LIBS output variable to hold the requisite library linkages.
#
# To link with BLAS, you should link with:
#
# $BLAS_LIBS $LIBS $FLIBS
#
# in that order. FLIBS is the output variable of the
# AC_F77_LIBRARY_LDFLAGS macro (called if necessary by AX_BLAS), and is
# sometimes necessary in order to link with F77 libraries. Users will also
# need to use AC_F77_DUMMY_MAIN (see the autoconf manual), for the same
# reason.
#
# Many libraries are searched for, from ATLAS to CXML to ESSL. The user
# may also use --with-blas=<lib> in order to use some specific BLAS
# library <lib>. In order to link successfully, however, be aware that you
# will probably need to use the same Fortran compiler (which can be set
# via the F77 env. var.) as was used to compile the BLAS library.
#
# ACTION-IF-FOUND is a list of shell commands to run if a BLAS library is
# found, and ACTION-IF-NOT-FOUND is a list of commands to run it if it is
# not found. If ACTION-IF-FOUND is not specified, the default action will
# define HAVE_BLAS.
#
# LICENSE
#
# Copyright (c) 2008 Steven G. Johnson <stevenj@alum.mit.edu>
# Copyright (c) 2019 Geoffrey M. Oxberry <goxberry@gmail.com>
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <https://www.gnu.org/licenses/>.
#
# As a special exception, the respective Autoconf Macro's copyright owner
# gives unlimited permission to copy, distribute and modify the configure
# scripts that are the output of Autoconf when processing the Macro. You
# need not follow the terms of the GNU General Public License when using
# or distributing such scripts, even though portions of the text of the
# Macro appear in them. The GNU General Public License (GPL) does govern
# all other use of the material that constitutes the Autoconf Macro.
#
# This special exception to the GPL applies to versions of the Autoconf
# Macro released by the Autoconf Archive. When you make and distribute a
# modified version of the Autoconf Macro, you may extend this special
# exception to the GPL to apply to your modified version as well.
#serial 17
AU_ALIAS([ACX_BLAS], [AX_BLAS])
AC_DEFUN([AX_BLAS], [
AC_PREREQ([2.55])
AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS])
AC_REQUIRE([AC_CANONICAL_HOST])
ax_blas_ok=no
AC_ARG_WITH(blas,
[AS_HELP_STRING([--with-blas=<lib>], [use BLAS library <lib>])])
case $with_blas in
yes | "") ;;
no) ax_blas_ok=disable ;;
-* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o)
BLAS_LIBS="$with_blas"
;;
*) BLAS_LIBS="-l$with_blas" ;;
esac
# Get fortran linker names of BLAS functions to check for.
AC_F77_FUNC(sgemm)
AC_F77_FUNC(dgemm)
ax_blas_save_LIBS="$LIBS"
LIBS="$LIBS $FLIBS"
# First, check BLAS_LIBS environment variable
if test $ax_blas_ok = no; then
if test "x$BLAS_LIBS" != x; then
save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS"
AC_MSG_CHECKING([for $sgemm in $BLAS_LIBS])
AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes], [BLAS_LIBS=""])
AC_MSG_RESULT($ax_blas_ok)
LIBS="$save_LIBS"
fi
fi
# BLAS linked to by default? (happens on some supercomputers)
if test $ax_blas_ok = no; then
save_LIBS="$LIBS"; LIBS="$LIBS"
AC_MSG_CHECKING([if $sgemm is being linked in already])
AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes])
AC_MSG_RESULT($ax_blas_ok)
LIBS="$save_LIBS"
fi
# BLAS in OpenBLAS library? (http://xianyi.github.com/OpenBLAS/)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(openblas, $sgemm, [ax_blas_ok=yes
BLAS_LIBS="-lopenblas"])
fi
# BLAS in ATLAS library? (http://math-atlas.sourceforge.net/)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(atlas, ATL_xerbla,
[AC_CHECK_LIB(f77blas, $sgemm,
[AC_CHECK_LIB(cblas, cblas_dgemm,
[ax_blas_ok=yes
BLAS_LIBS="-lcblas -lf77blas -latlas"],
[], [-lf77blas -latlas])],
[], [-latlas])])
fi
# BLAS in PhiPACK libraries? (requires generic BLAS lib, too)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(blas, $sgemm,
[AC_CHECK_LIB(dgemm, $dgemm,
[AC_CHECK_LIB(sgemm, $sgemm,
[ax_blas_ok=yes; BLAS_LIBS="-lsgemm -ldgemm -lblas"],
[], [-lblas])],
[], [-lblas])])
fi
# BLAS in Intel MKL library?
if test $ax_blas_ok = no; then
# MKL for gfortran
if test x"$ac_cv_fc_compiler_gnu" = xyes; then
# 64 bit
if test $host_cpu = x86_64; then
AC_CHECK_LIB(mkl_gf_lp64, $sgemm,
[ax_blas_ok=yes;BLAS_LIBS="-lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread"],,
[-lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread])
# 32 bit
elif test $host_cpu = i686; then
AC_CHECK_LIB(mkl_gf, $sgemm,
[ax_blas_ok=yes;BLAS_LIBS="-lmkl_gf -lmkl_sequential -lmkl_core -lpthread"],,
[-lmkl_gf -lmkl_sequential -lmkl_core -lpthread])
fi
# MKL for other compilers (Intel, PGI, ...?)
else
# 64-bit
if test $host_cpu = x86_64; then
AC_CHECK_LIB(mkl_intel_lp64, $sgemm,
[ax_blas_ok=yes;BLAS_LIBS="-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread"],,
[-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread])
# 32-bit
elif test $host_cpu = i686; then
AC_CHECK_LIB(mkl_intel, $sgemm,
[ax_blas_ok=yes;BLAS_LIBS="-lmkl_intel -lmkl_sequential -lmkl_core -lpthread"],,
[-lmkl_intel -lmkl_sequential -lmkl_core -lpthread])
fi
fi
fi
# Old versions of MKL
if test $ax_blas_ok = no; then
AC_CHECK_LIB(mkl, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lmkl -lguide -lpthread"],,[-lguide -lpthread])
fi
# BLAS in Apple vecLib library?
if test $ax_blas_ok = no; then
save_LIBS="$LIBS"; LIBS="-framework vecLib $LIBS"
AC_MSG_CHECKING([for $sgemm in -framework vecLib])
AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"])
AC_MSG_RESULT($ax_blas_ok)
LIBS="$save_LIBS"
fi
# BLAS in Alpha CXML library?
if test $ax_blas_ok = no; then
AC_CHECK_LIB(cxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-lcxml"])
fi
# BLAS in Alpha DXML library? (now called CXML, see above)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(dxml, $sgemm, [ax_blas_ok=yes;BLAS_LIBS="-ldxml"])
fi
# BLAS in Sun Performance library?
if test $ax_blas_ok = no; then
if test "x$GCC" != xyes; then # only works with Sun CC
AC_CHECK_LIB(sunmath, acosp,
[AC_CHECK_LIB(sunperf, $sgemm,
[BLAS_LIBS="-xlic_lib=sunperf -lsunmath"
ax_blas_ok=yes],[],[-lsunmath])])
fi
fi
# BLAS in SCSL library? (SGI/Cray Scientific Library)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(scs, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lscs"])
fi
# BLAS in SGIMATH library?
if test $ax_blas_ok = no; then
AC_CHECK_LIB(complib.sgimath, $sgemm,
[ax_blas_ok=yes; BLAS_LIBS="-lcomplib.sgimath"])
fi
# BLAS in IBM ESSL library? (requires generic BLAS lib, too)
if test $ax_blas_ok = no; then
AC_CHECK_LIB(blas, $sgemm,
[AC_CHECK_LIB(essl, $sgemm,
[ax_blas_ok=yes; BLAS_LIBS="-lessl -lblas"],
[], [-lblas $FLIBS])])
fi
# Generic BLAS library?
if test $ax_blas_ok = no; then
AC_CHECK_LIB(blas, $sgemm, [ax_blas_ok=yes; BLAS_LIBS="-lblas"])
fi
AC_SUBST(BLAS_LIBS)
LIBS="$ax_blas_save_LIBS"
# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
if test x"$ax_blas_ok" = xyes; then
ifelse([$1],,AC_DEFINE(HAVE_BLAS,1,[Define if you have a BLAS library.]),[$1])
:
else
ax_blas_ok=no
$2
fi
])dnl AX_BLAS

96
m4/ax_check_gnu_make.m4 Normal file
View File

@ -0,0 +1,96 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_check_gnu_make.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_CHECK_GNU_MAKE([run-if-true],[run-if-false])
#
# DESCRIPTION
#
# This macro searches for a GNU version of make. If a match is found:
#
# * The makefile variable `ifGNUmake' is set to the empty string, otherwise
# it is set to "#". This is useful for including a special features in a
# Makefile, which cannot be handled by other versions of make.
# * The makefile variable `ifnGNUmake' is set to #, otherwise
# it is set to the empty string. This is useful for including a special
# features in a Makefile, which can be handled
# by other versions of make or to specify else like clause.
# * The variable `_cv_gnu_make_command` is set to the command to invoke
# GNU make if it exists, the empty string otherwise.
# * The variable `ax_cv_gnu_make_command` is set to the command to invoke
# GNU make by copying `_cv_gnu_make_command`, otherwise it is unset.
# * If GNU Make is found, its version is extracted from the output of
# `make --version` as the last field of a record of space-separated
# columns and saved into the variable `ax_check_gnu_make_version`.
# * Additionally if GNU Make is found, run shell code run-if-true
# else run shell code run-if-false.
#
# Here is an example of its use:
#
# Makefile.in might contain:
#
# # A failsafe way of putting a dependency rule into a makefile
# $(DEPEND):
# $(CC) -MM $(srcdir)/*.c > $(DEPEND)
#
# @ifGNUmake@ ifeq ($(DEPEND),$(wildcard $(DEPEND)))
# @ifGNUmake@ include $(DEPEND)
# @ifGNUmake@ else
# fallback code
# @ifGNUmake@ endif
#
# Then configure.in would normally contain:
#
# AX_CHECK_GNU_MAKE()
# AC_OUTPUT(Makefile)
#
# Then perhaps to cause gnu make to override any other make, we could do
# something like this (note that GNU make always looks for GNUmakefile
# first):
#
# if ! test x$_cv_gnu_make_command = x ; then
# mv Makefile GNUmakefile
# echo .DEFAULT: > Makefile ;
# echo \ $_cv_gnu_make_command \$@ >> Makefile;
# fi
#
# Then, if any (well almost any) other make is called, and GNU make also
# exists, then the other make wraps the GNU make.
#
# LICENSE
#
# Copyright (c) 2008 John Darrington <j.darrington@elvis.murdoch.edu.au>
# Copyright (c) 2015 Enrico M. Crisostomo <enrico.m.crisostomo@gmail.com>
#
# Copying and distribution of this file, with or without modification, are
# permitted in any medium without royalty provided the copyright notice
# and this notice are preserved. This file is offered as-is, without any
# warranty.
#serial 12
AC_DEFUN([AX_CHECK_GNU_MAKE],dnl
[AC_PROG_AWK
AC_CACHE_CHECK([for GNU make],[_cv_gnu_make_command],[dnl
_cv_gnu_make_command="" ;
dnl Search all the common names for GNU make
for a in "$MAKE" make gmake gnumake ; do
if test -z "$a" ; then continue ; fi ;
if "$a" --version 2> /dev/null | grep GNU 2>&1 > /dev/null ; then
_cv_gnu_make_command=$a ;
AX_CHECK_GNU_MAKE_HEADLINE=$("$a" --version 2> /dev/null | grep "GNU Make")
ax_check_gnu_make_version=$(echo ${AX_CHECK_GNU_MAKE_HEADLINE} | ${AWK} -F " " '{ print $(NF); }')
break ;
fi
done ;])
dnl If there was a GNU version, then set @ifGNUmake@ to the empty string, '#' otherwise
AS_VAR_IF([_cv_gnu_make_command], [""], [AS_VAR_SET([ifGNUmake], ["#"])], [AS_VAR_SET([ifGNUmake], [""])])
AS_VAR_IF([_cv_gnu_make_command], [""], [AS_VAR_SET([ifnGNUmake], [""])], [AS_VAR_SET([ifnGNUmake], ["#"])])
AS_VAR_IF([_cv_gnu_make_command], [""], [AS_UNSET(ax_cv_gnu_make_command)], [AS_VAR_SET([ax_cv_gnu_make_command], [${_cv_gnu_make_command}])])
AS_VAR_IF([_cv_gnu_make_command], [""],[$2],[$1])
AC_SUBST([ifGNUmake])
AC_SUBST([ifnGNUmake])
])

73
m4/ax_irpf90.m4 Normal file
View File

@ -0,0 +1,73 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_irpf90.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_PROG_IRPF90
#
# DESCRIPTION
#
# Check for the program 'irpf90', let script continue if exists, pops up
# error message if not.
#
# Besides checking existence, this macro also set these environment
# variables upon completion:
#
# IRPF90 = which irpf90
#
# DEPENDENCIES
#
# AX_CHECK_GNU_MAKE
#
# LICENSE
#
# Copyright (c) 2021 Anthony Scemama
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <https://www.gnu.org/licenses/>.
#
# As a special exception, the respective Autoconf Macro's copyright owner
# gives unlimited permission to copy, distribute and modify the configure
# scripts that are the output of Autoconf when processing the Macro. You
# need not follow the terms of the GNU General Public License when using
# or distributing such scripts, even though portions of the text of the
# Macro appear in them. The GNU General Public License (GPL) does govern
# all other use of the material that constitutes the Autoconf Macro.
#
# This special exception to the GPL applies to versions of the Autoconf
# Macro released by the Autoconf Archive. When you make and distribute a
# modified version of the Autoconf Macro, you may extend this special
# exception to the GPL to apply to your modified version as well.
AU_ALIAS([AC_PROG_IRPF90], [AX_PROG_IRPF90])
AC_DEFUN([AX_PROG_IRPF90], [
# Requirements
AC_REQUIRE([AX_CHECK_GNU_MAKE])
AS_IF([test "x$ifGNUmake" = "x#"], [ AC_MSG_ERROR([GNU Make (gmake) is required with IRPF90]) ])
# IRPF90
AC_PATH_PROG([IRPF90], [irpf90], [nocommand])
AS_IF([test "x$IRPF90" = xnocommand], [
AC_MSG_ERROR([irpf90 not found in $PATH]) ])
AC_FC_FREEFORM
AC_FC_LINE_LENGTH
AC_FC_MODULE_EXTENSION
])

135
m4/ax_lapack.m4 Normal file
View File

@ -0,0 +1,135 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_lapack.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_LAPACK([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]])
#
# DESCRIPTION
#
# This macro looks for a library that implements the LAPACK linear-algebra
# interface (see http://www.netlib.org/lapack/). On success, it sets the
# LAPACK_LIBS output variable to hold the requisite library linkages.
#
# To link with LAPACK, you should link with:
#
# $LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS
#
# in that order. BLAS_LIBS is the output variable of the AX_BLAS macro,
# called automatically. FLIBS is the output variable of the
# AC_F77_LIBRARY_LDFLAGS macro (called if necessary by AX_BLAS), and is
# sometimes necessary in order to link with F77 libraries. Users will also
# need to use AC_F77_DUMMY_MAIN (see the autoconf manual), for the same
# reason.
#
# The user may also use --with-lapack=<lib> in order to use some specific
# LAPACK library <lib>. In order to link successfully, however, be aware
# that you will probably need to use the same Fortran compiler (which can
# be set via the F77 env. var.) as was used to compile the LAPACK and BLAS
# libraries.
#
# ACTION-IF-FOUND is a list of shell commands to run if a LAPACK library
# is found, and ACTION-IF-NOT-FOUND is a list of commands to run it if it
# is not found. If ACTION-IF-FOUND is not specified, the default action
# will define HAVE_LAPACK.
#
# LICENSE
#
# Copyright (c) 2009 Steven G. Johnson <stevenj@alum.mit.edu>
# Copyright (c) 2019 Geoffrey M. Oxberry <goxberry@gmail.com>
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <https://www.gnu.org/licenses/>.
#
# As a special exception, the respective Autoconf Macro's copyright owner
# gives unlimited permission to copy, distribute and modify the configure
# scripts that are the output of Autoconf when processing the Macro. You
# need not follow the terms of the GNU General Public License when using
# or distributing such scripts, even though portions of the text of the
# Macro appear in them. The GNU General Public License (GPL) does govern
# all other use of the material that constitutes the Autoconf Macro.
#
# This special exception to the GPL applies to versions of the Autoconf
# Macro released by the Autoconf Archive. When you make and distribute a
# modified version of the Autoconf Macro, you may extend this special
# exception to the GPL to apply to your modified version as well.
#serial 10
AU_ALIAS([ACX_LAPACK], [AX_LAPACK])
AC_DEFUN([AX_LAPACK], [
AC_REQUIRE([AX_BLAS])
ax_lapack_ok=no
AC_ARG_WITH(lapack,
[AS_HELP_STRING([--with-lapack=<lib>], [use LAPACK library <lib>])])
case $with_lapack in
yes | "") ;;
no) ax_lapack_ok=disable ;;
-* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o)
LAPACK_LIBS="$with_lapack"
;;
*) LAPACK_LIBS="-l$with_lapack" ;;
esac
# Get fortran linker name of LAPACK function to check for.
AC_F77_FUNC(cheev)
# We cannot use LAPACK if BLAS is not found
if test "x$ax_blas_ok" != xyes; then
ax_lapack_ok=noblas
LAPACK_LIBS=""
fi
# First, check LAPACK_LIBS environment variable
if test "x$LAPACK_LIBS" != x; then
save_LIBS="$LIBS"; LIBS="$LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS"
AC_MSG_CHECKING([for $cheev in $LAPACK_LIBS])
AC_LINK_IFELSE([AC_LANG_CALL([], [$cheev])], [ax_lapack_ok=yes], [LAPACK_LIBS=""])
AC_MSG_RESULT($ax_lapack_ok)
LIBS="$save_LIBS"
if test $ax_lapack_ok = no; then
LAPACK_LIBS=""
fi
fi
# LAPACK linked to by default? (is sometimes included in BLAS lib)
if test $ax_lapack_ok = no; then
save_LIBS="$LIBS"; LIBS="$LIBS $BLAS_LIBS $FLIBS"
AC_CHECK_FUNC($cheev, [ax_lapack_ok=yes])
LIBS="$save_LIBS"
fi
# Generic LAPACK library?
for lapack in lapack lapack_rs6k; do
if test $ax_lapack_ok = no; then
save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS"
AC_CHECK_LIB($lapack, $cheev,
[ax_lapack_ok=yes; LAPACK_LIBS="-l$lapack"], [], [$FLIBS])
LIBS="$save_LIBS"
fi
done
AC_SUBST(LAPACK_LIBS)
# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
if test x"$ax_lapack_ok" = xyes; then
ifelse([$1],,AC_DEFINE(HAVE_LAPACK,1,[Define if you have LAPACK library.]),[$1])
:
else
ax_lapack_ok=no
$2
fi
])dnl AX_LAPACK

91
m4/ax_zmq.m4 Normal file
View File

@ -0,0 +1,91 @@
# ===========================================================================
# https://www.gnu.org/software/autoconf-archive/ax_zmq.html
# ===========================================================================
#
# SYNOPSIS
#
# AX_ZMQ([MINIMUM-VERSION], [ACTION-IF-FOUND], [ACTION-IF-NOT-FOUND])
#
# DESCRIPTION
#
# Test for the ZMQ libraries of a particular version (or newer). The
# default version tested for is 4.0.0.
#
# The macro tests for ZMQ libraries in the library/include path, and, when
# provided, also in the path given by --with-zmq.
#
# This macro calls:
#
# AC_SUBST(ZMQ_CPPFLAGS) / AC_SUBST(ZMQ_LDFLAGS) / AC_SUBST(ZMQ_LIBS)
#
# And sets:
#
# HAVE_ZMQ
#
# LICENSE
#
# Copyright (c) 2016 Jeroen Meijer <jjgmeijer@gmail.com>
#
# Copying and distribution of this file, with or without modification, are
# permitted in any medium without royalty provided the copyright notice
# and this notice are preserved. This file is offered as-is, without any
# warranty.
#serial 3
AC_DEFUN([AX_ZMQ], [
AC_ARG_WITH([zmq], [AS_HELP_STRING([--with-zmq=<prefix>],[ZMQ prefix directory])], [
ZMQ_LDFLAGS="-L${with_zmq}/lib"
ZMQ_CPPFLAGS="-I${with_zmq}/include"
])
HAVE_ZMQ=0
if test "$with_zmq" != "no"; then
LD_FLAGS="$LDFLAGS $ZMQ_LDFLAGS"
CPPFLAGS="$CPPFLAGS $ZMQ_CPPFLAGS"
AC_LANG_PUSH([C])
AC_CHECK_HEADER(zmq.h, [zmq_h=yes], [zmq_h=no])
AC_LANG_POP([C])
if test "$zmq_h" = "yes"; then
version=ifelse([$1], ,4.0.0,$1)
AC_MSG_CHECKING([for ZMQ version >= $version])
version=$(echo $version | tr '.' ',')
AC_EGREP_CPP([version_ok], [
#include <zmq.h>
#if defined(ZMQ_VERSION) && ZMQ_VERSION >= ZMQ_MAKE_VERSION($version)
version_ok
#endif
],[
AC_MSG_RESULT(yes)
AC_CHECK_LIB([zmq],[zmq_send],[
HAVE_ZMQ=1
ZMQ_LIBS="-lzmq"
],[
AC_MSG_ERROR([libzmq.so not in LIBRARY_PATH])
])
AC_SUBST(ZMQ_LDFLAGS)
AC_SUBST(ZMQ_CPPFLAGS)
AC_SUBST(ZMQ_LIBS)
], AC_MSG_RESULT([no valid ZMQ version was found]))
else
AC_MSG_WARN([no valid ZMQ installation was found])
fi
if test $HAVE_ZMQ = 1; then
# execute ACTION-IF-FOUND (if present):
ifelse([$2], , :, [$2])
else
# execute ACTION-IF-NOT-FOUND (if present):
ifelse([$3], , :, [$3])
fi
else
AC_MSG_NOTICE([not checking for ZMQ])
fi
AC_DEFINE(HAVE_ZMQ,,[define if the ZMQ library is available])
])

218
ocaml/Block.ml
View File

@ -1,6 +1,6 @@
open Qptypes
type t =
type t =
{ property : Property.t ;
value : Sample.t ;
weight : Weight.t ;
@ -12,7 +12,7 @@ type t =
let re =
Str.regexp "[ |#|\n]+"
let of_string s =
let of_string s =
try
let lst =
@ -23,24 +23,24 @@ let of_string s =
| b :: pid :: c:: p :: w :: v :: [] -> Some
{ property = Property.of_string p ;
value = Sample.of_float (float_of_string v) ;
weight = Weight.of_float (float_of_string w) ;
weight = Weight.of_float (float_of_string w) ;
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b) ;
}
| b :: pid :: c:: p :: w :: v ->
let v =
| b :: pid :: c:: p :: w :: v ->
let v =
List.rev v
|> Array.of_list
|> Array.of_list
|> Array.map float_of_string
in
let dim =
let dim =
Array.length v
in
Some
{ property = Property.of_string p ;
value = Sample.of_float_array ~dim v ;
weight = Weight.of_float (float_of_string w) ;
weight = Weight.of_float (float_of_string w) ;
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b) ;
@ -50,20 +50,105 @@ let of_string s =
| _ -> None
let to_short_string b =
Printf.sprintf "%s # %s %d %d"
(Property.to_string b.property)
(Compute_node.to_string b.compute_node)
b.pid
(Block_id.to_int b.block_id)
let to_string b =
Printf.sprintf "%s %s # %s %s %s %d"
(Sample.to_string b.value )
(Weight.to_float b.weight |> string_of_float)
(Property.to_string b.property)
(Compute_node.to_string b.compute_node)
(string_of_int b.pid)
(Block_id.to_int b.block_id)
Printf.sprintf "%s %s # %s %s %s %d"
(Sample.to_string b.value )
(Weight.to_float b.weight |> string_of_float)
(Property.to_string b.property)
(Compute_node.to_string b.compute_node)
(string_of_int b.pid)
(Block_id.to_int b.block_id)
let zero =
bytes_of_int 0
let to_bytes b =
(* [ Length of b
[ Length of value ;
Value ;
Length of weight ;
Weight ;
... ] ] *)
let l =
[ Property.to_bytes b.property ;
Sample.to_bytes b.value ;
Weight.to_bytes b.weight ;
bytes_of_int b.pid ;
Block_id.to_bytes b.block_id ;
Compute_node.to_bytes b.compute_node ]
|> List.map (fun x -> [ bytes_of_int (Bytes.length x) ; x ] )
|> List.concat
in
let result =
Bytes.concat Bytes.empty (zero :: l)
in
Bytes.set_int64_le result 0 (Int64.of_int ((Bytes.length result) - 8));
result
let read_bytes b =
(* Reads m, the first 8 bytes as an int64 containing the number of bytes to read.
Then, read the next m bytes and return a tuple containing the decoded data and the rest.
*)
let l = Bytes.length b in
if l < 8 then
failwith "Zero-sized bytes"
else
let m =
Bytes.get_int64_le b 0
|> Int64.to_int
in
let nl = l-m-8 in
if nl > 0 then
(Bytes.sub b 8 m, Some (Bytes.sub b (8+m) nl))
else
(Bytes.sub b 8 m, None)
let of_bytes b =
let rec loop accu s =
match read_bytes s with
| data, None -> List.rev (data :: accu)
| data, (Some rest) -> loop (data :: accu) rest
in
let result =
match loop [] b with
| property :: value :: weight :: pid :: block_id :: compute_node :: [] ->
Some
{ property = Property.of_bytes property;
value = Sample.of_bytes value;
weight = Weight.of_bytes weight;
pid = int_of_bytes pid;
block_id = Block_id.of_bytes block_id;
compute_node = Compute_node.of_bytes compute_node;
}
| _ -> None
in
result
let of_string_or_bytes s =
if Qmcchem_config.binary_io then
Bytes.of_string s
|> of_bytes
else
of_string s
let dir_name = lazy(
let ezfio_filename =
let ezfio_filename =
Lazy.force Qputils.ezfio_filename
in
let md5 =
@ -72,8 +157,8 @@ let dir_name = lazy(
let d = Filename.concat ezfio_filename "blocks" in
if not ( Sys.file_exists d ) then
Unix.mkdir d 0o755;
List.fold_right Filename.concat
[ ezfio_filename ; "blocks" ; md5 ; Filename.dir_sep ] ""
List.fold_right Filename.concat
[ ezfio_filename ; "blocks" ; md5 ; Filename.dir_sep ] ""
)
@ -84,11 +169,11 @@ let _raw_data =
let update_raw_data ?(locked=true) () =
(* Create array of files to read *)
let dir_name =
let dir_name =
Lazy.force dir_name
in
let files =
let result =
let files =
let result =
if Sys.file_exists dir_name && Sys.is_directory dir_name then
begin
Sys.readdir dir_name
@ -102,7 +187,7 @@ let update_raw_data ?(locked=true) () =
else
List.filter (fun x ->
try
let _ =
let _ =
Str.search_backward (Str.regexp "locked") x ((String.length x) - 1)
in false
with
@ -110,7 +195,7 @@ let update_raw_data ?(locked=true) () =
) result
in
let rec transform new_list = function
let rec transform new_list = function
| [] -> new_list
| head :: tail ->
let head = String.trim head in
@ -120,33 +205,74 @@ let update_raw_data ?(locked=true) () =
| Some x -> transform (x::new_list) tail
in
let result =
let rec aux ic accu =
let l =
try
Some (input_line ic)
with
| End_of_file -> None
if Qmcchem_config.binary_io then
begin
let result =
let rec aux buf accu =
(* Read one block *)
let item, rest =
read_bytes buf
in
match of_bytes item with
| None -> []
| Some item ->
match rest with
| None -> List.rev (item::accu)
| Some rest -> (aux [@tailcall]) rest (item::accu)
in
List.concat_map (fun filename ->
let ic = open_in filename in
let length = in_channel_length ic in
let result =
if length > 0 then
let buf = Bytes.create length in
really_input ic buf 0 length;
aux buf []
else []
in
close_in ic;
result ) files
in
match l with
| None -> List.rev accu
| Some l -> (aux [@tailcall]) ic (l::accu)
in
List.concat_map (fun filename ->
let ic = open_in filename in
let result = aux ic [] in
close_in ic;
result ) files
|> transform []
in
result
result
end
else
begin
let result =
let rec aux ic accu =
let l =
try
Some (input_line ic)
with
| End_of_file -> None
in
match l with
| None -> List.rev accu
| Some l -> (aux [@tailcall]) ic (l::accu)
in
List.concat_map (fun filename ->
let ic = open_in filename in
let result = aux ic [] in
close_in ic;
result ) files
|> transform []
in
result
end
let raw_data ?(locked=true) () =
let to_string_or_bytes b =
if Qmcchem_config.binary_io then
to_bytes b
|> Bytes.to_string
else
to_string b
let raw_data ?(locked=true) () =
match !_raw_data with
| Some x -> x
| None ->
let result =
let result =
update_raw_data ~locked ()
in
_raw_data := Some result;
@ -156,7 +282,7 @@ let raw_data ?(locked=true) () =
let properties = lazy (
let h = Hashtbl.create 63 in
List.iter (fun x ->
List.iter (fun x ->
Hashtbl.replace h (Property.to_string x.property) x.property)
(raw_data ());
Hashtbl.fold (fun k v a -> v :: a) h []

34
ocaml/Message.ml
View File

@ -11,7 +11,7 @@ type t =
| Error of string
let create m =
let create m =
try
match m with
| [ "cpu" ; c ; pid ; b ; "1" ; v ] ->
@ -23,7 +23,7 @@ let create m =
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b);
}
}
| [ "accep" ; c ; pid ; b ; "1" ; v ] ->
let open Block in
Property
@ -33,8 +33,8 @@ let create m =
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b);
}
| [ prop ; c ; pid ; b ; w ; v ] ->
}
| [ prop ; c ; pid ; b ; w ; v ] ->
let open Block in
Property
{ property = Property.of_string prop;
@ -43,19 +43,19 @@ let create m =
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b);
}
}
| "elec_coord" :: c :: pid :: _ :: n ::walkers ->
begin
let elec_num =
Lazy.force Qputils.elec_num
and n =
and n =
int_of_string n
in
assert (n = List.length walkers);
let rec build_walker accu = function
| (0,tail) ->
let result =
List.rev accu
| (0,tail) ->
let result =
List.rev accu
|> List.rev_map float_of_string
|> List.rev
|> Array.of_list
@ -65,10 +65,10 @@ let create m =
build_walker (head::accu) (n-1, tail)
| _ -> failwith "Bad walkers"
in
let rec build accu = function
let rec build accu = function
| [] -> Array.of_list accu
| w ->
let (result, tail) =
| w ->
let (result, tail) =
build_walker [] (3*elec_num+3, w)
in
build (result::accu) tail
@ -80,13 +80,13 @@ let create m =
| [ "unregister" ; c ; pid ] -> Unregister (Compute_node.of_string c, int_of_string pid)
| [ "Test" ] -> Test
| [ "Ezfio" ; ezfio_msg ] -> Ezfio ezfio_msg
| prop :: c :: pid :: b :: d :: w :: l ->
| prop :: c :: pid :: b :: d :: w :: l ->
let property =
Property.of_string prop
in
begin
assert (not (Property.is_scalar property));
let a =
let a =
Array.of_list l
|> Array.map float_of_string
and dim =
@ -101,7 +101,7 @@ let create m =
compute_node = Compute_node.of_string c;
pid = int_of_string pid;
block_id = Block_id.of_int (int_of_string b);
}
}
end
| l -> Error (String.concat ":" l)
with
@ -114,9 +114,9 @@ let to_string = function
| Walkers (h,p,w) -> Printf.sprintf "Walkers : %s %d : %d walkers"
(Compute_node.to_string h) p (Array.length w)
| GetWalkers n -> Printf.sprintf "GetWalkers %d" (Strictly_positive_int.to_int n)
| Register (h,p) -> Printf.sprintf "Register : %s %d"
| Register (h,p) -> Printf.sprintf "Register : %s %d"
(Compute_node.to_string h) p
| Unregister (h,p) -> Printf.sprintf "Unregister : %s %d"
| Unregister (h,p) -> Printf.sprintf "Unregister : %s %d"
(Compute_node.to_string h) p
| Test -> "Test"
| Ezfio msg -> "Ezfio "^msg

19
ocaml/Qmcchem_config.ml
View File

@ -8,10 +8,10 @@ let root = lazy (
(* PATH environment variable as a list of strings *)
let path = lazy (
let p =
let p =
try Sys.getenv "PATH" with
| Not_found -> failwith "PATH environment variable is not set"
in
in
String.split_on_char ':' p
)
@ -23,8 +23,8 @@ let full_path exe =
| [] -> None
| head :: tail ->
begin
let fp =
Filename.concat head exe
let fp =
Filename.concat head exe
in
if Sys.file_exists fp then
Some fp
@ -82,8 +82,8 @@ let hostname = lazy (
)
external get_ipv4_address_for_interface : string -> string =
"get_ipv4_address_for_interface"
external get_ipv4_address_for_interface : string -> string =
"get_ipv4_address_for_interface"
let ip_address = lazy (
@ -106,7 +106,7 @@ let ip_address = lazy (
end
| Some interface ->
let result = get_ipv4_address_for_interface interface in
if String.sub result 0 5 = "error" then
if String.sub result 0 5 = "error" then
Printf.sprintf "Unable to use network interface %s" interface
|> failwith
else
@ -115,3 +115,8 @@ let ip_address = lazy (
let binary_io =
try
Sys.getenv "QMCCHEM_IO" = "B"
with Not_found -> false

264
ocaml/Qmcchem_dataserver.ml
View File

@ -17,9 +17,9 @@ let initialization_timeout = 600.
let bind_socket ~socket_type ~socket ~address =
try
Zmq.Socket.bind socket address
with
| Unix.Unix_error (_, message, f) ->
failwith @@ Printf.sprintf
with
| Unix.Unix_error (_, message, f) ->
failwith @@ Printf.sprintf
"\n%s\nUnable to bind the dataserver's %s socket :\n %s\n%s"
f socket_type address message
| other_exception -> raise other_exception
@ -30,7 +30,7 @@ let run ?(daemon=true) ezfio_filename =
Qputils.set_ezfio_filename ezfio_filename;
(** Check if walkers need to be created. *)
let () =
let () =
if ( not(Ezfio.has_electrons_elec_coord_pool ()) ) then
begin
Printf.printf "Generating initial walkers...\n%!";
@ -39,12 +39,12 @@ let run ?(daemon=true) ezfio_filename =
Unix.execvp
(Lazy.force Qmcchem_config.qmc_create_walkers)
[|"qmc_create_walkers" ; ezfio_filename|]
| pid ->
| pid ->
begin
ignore @@ Unix.waitpid [] pid;
Printf.printf "Initial walkers ready\n%!"
end
end
end
in
@ -74,7 +74,7 @@ let run ?(daemon=true) ezfio_filename =
in
(** Status variable (mutable) *)
let status =
let status =
ref Status.Queued
in
@ -91,7 +91,7 @@ let run ?(daemon=true) ezfio_filename =
let adress_prefix =
"tcp://*:"
in
let result =
let result =
List.fold_left (fun accu i ->
let address =
adress_prefix ^ (string_of_int (n+i))
@ -99,14 +99,14 @@ let run ?(daemon=true) ezfio_filename =
let socket =
Zmq.Socket.create zmq_context Zmq.Socket.rep
in
let result =
let result =
try
Zmq.Socket.bind socket address;
accu
with
| _ -> false
in
Zmq.Socket.close socket;
Zmq.Socket.close socket;
result
) true [0;1;2;3]
in
@ -120,7 +120,7 @@ let run ?(daemon=true) ezfio_filename =
(** Random port number between 49152 and 65535 *)
let port =
let newport =
let newport =
ref ( 1024 + (Random.int (49151-1024)))
in
while ((check_port !newport) = `Unavailable)
@ -131,7 +131,7 @@ let run ?(daemon=true) ezfio_filename =
in
let debug_socket =
let debug_socket =
Zmq.Socket.create zmq_context Zmq.Socket.xpub
and address =
Printf.sprintf "tcp://*:%d" (port+4)
@ -140,17 +140,17 @@ let run ?(daemon=true) ezfio_filename =
Zmq.Socket.set_linger_period debug_socket 100 ;
let close_debug_socket () =
Zmq.Socket.close debug_socket
Zmq.Socket.close debug_socket
in
(** Sends a log text to the debug socket *)
let send_log socket size t0 text =
let dt =
let dt =
delta_t t0
in
let message =
Printf.sprintf "%20s : %8d : %10s : %e"
socket size text dt
let message =
Printf.sprintf "%20s : %8d : %10s : %e"
socket size text dt
in
Zmq.Socket.send debug_socket message
in
@ -159,7 +159,7 @@ let run ?(daemon=true) ezfio_filename =
(** {2 Walkers} *)
(** Number of electrons *)
let elec_num =
let elec_num =
Lazy.force Qputils.elec_num
in
@ -181,14 +181,14 @@ let run ?(daemon=true) ezfio_filename =
in
(** Array of walkers. The size is [walk_num_tot]. *)
let walkers_array =
let t0 =
let walkers_array =
let t0 =
Unix.gettimeofday ()
in
let j =
3*elec_num + 3
in
let result =
let result =
let size =
Ezfio.get_electrons_elec_coord_pool_size ()
and ez =
@ -204,7 +204,7 @@ let run ?(daemon=true) ezfio_filename =
failwith "Walkers file is broken."
in
String.concat " " [ "Read" ; string_of_int (Array.length result) ;
"walkers"]
"walkers"]
|> send_log "status" 0 t0 ;
result
in
@ -217,7 +217,7 @@ let run ?(daemon=true) ezfio_filename =
(** Last time when the walkers were saved to disk. *)
let last_save_walkers =
let last_save_walkers =
ref 0.
in
@ -227,14 +227,12 @@ let run ?(daemon=true) ezfio_filename =
if (delta_t !last_save_walkers > 10. ) then
begin
Ezfio.set_electrons_elec_coord_pool_size walk_num_tot ;
let walkers_list =
Array.map Array.to_list walkers_array
|> Array.to_list
|> List.concat
|> List.rev_map float_of_string
|> List.rev
let walkers_list =
Array.to_list walkers_array
|> Array.concat
|> Array.map float_of_string
in
Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_list
Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_array
~rank:3 ~dim:[| elec_num+1 ; 3 ; walk_num_tot |] ~data:walkers_list);
let t0 =
Unix.gettimeofday ()
@ -263,7 +261,7 @@ let run ?(daemon=true) ezfio_filename =
(** {2 Set of workers} *)
(** A hash table is kept to track the running workers. The keys are the
(** A hash table is kept to track the running workers. The keys are the
built as string containing the couple ([Compute_node], [PID]), and
the values are the last communication time.
*)
@ -271,16 +269,16 @@ let run ?(daemon=true) ezfio_filename =
(** The hash table for workers *)
let workers_hash =
let workers_hash =
Hashtbl.create 63
in
(** Creates a key using the couple ([Compute_node], [PID]) *)
let key compute_node pid =
String.concat " " [
(Compute_node.to_string compute_node);
let key compute_node pid =
String.concat " " [
(Compute_node.to_string compute_node);
(string_of_int pid) ]
in
@ -294,7 +292,7 @@ let run ?(daemon=true) ezfio_filename =
in
match Hashtbl.find_opt workers_hash s with
| Some _ -> failwith (s^" already registered")
| None -> Hashtbl.add workers_hash s (Unix.gettimeofday ())
| None -> Hashtbl.add workers_hash s (Unix.gettimeofday ())
in
@ -307,14 +305,14 @@ let run ?(daemon=true) ezfio_filename =
in
match Hashtbl.find_opt workers_hash s with
| Some x -> Hashtbl.remove workers_hash s
| None -> failwith (s^" not registered")
| None -> failwith (s^" not registered")
in
(** Sets the last access of the worker to [Unix.gettimeofday ()] *)
let touch_worker w pid =
let s =
let touch_worker w pid =
let s =
key w pid
in
Hashtbl.replace workers_hash s (Unix.gettimeofday ())
@ -326,11 +324,11 @@ let run ?(daemon=true) ezfio_filename =
let delta =
initialization_timeout +. block_time *. 2.
in
Hashtbl.fold (fun k v accu ->
Hashtbl.fold (fun k v accu ->
if (now -. v) <= delta then
v :: accu
else
accu ) hash []
accu ) hash []
|> List.length
in
@ -344,14 +342,14 @@ let run ?(daemon=true) ezfio_filename =
(** Name of the blocks file written by the current process. *)
let block_channel_filename =
let dirname =
let block_channel_filename =
let dirname =
Lazy.force Block.dir_name
in
let () =
let () =
if not ( Sys.file_exists dirname ) then
Unix.mkdir dirname 0o755
in
in
Filename.concat dirname (
hostname ^ "." ^ (string_of_int dataserver_pid)
)
@ -387,7 +385,7 @@ let run ?(daemon=true) ezfio_filename =
the compute nodes. Happens when [max_file_size] is reached.
*)
let compress_block_file filename =
let t0 =
let t0 =
Unix.gettimeofday ()
in
close_out !block_channel;
@ -406,16 +404,16 @@ let run ?(daemon=true) ezfio_filename =
(** {3 Status thread} *)
let start_status_thread =
let start_status_thread =
let t0 =
Unix.gettimeofday ()
in
Thread.create (fun () ->
Thread.create (fun () ->
send_log "status" 0 t0 "Starting status thread";
let socket =
Zmq.Socket.create zmq_context Zmq.Socket.pub
and address =
and address =
Printf.sprintf "tcp://*:%d" (port+1)
in
bind_socket "PUB" socket address;
@ -423,15 +421,15 @@ let run ?(daemon=true) ezfio_filename =
and delay_read = 2.
in
let start_time =
let start_time =
Unix.gettimeofday ()
and stop_time =
and stop_time =
ref (Input.Stop_time.(read () |> to_float) )
in
let last_update =
ref start_time
in
in
while (!status <> Status.Stopped)
do
@ -439,46 +437,46 @@ let run ?(daemon=true) ezfio_filename =
let now =
Unix.gettimeofday ()
in
let status_string =
let status_string =
Status.to_string !status
in
Zmq.Socket.send socket status_string;
send_log "status" (String.length status_string) now status_string;
let test =
let test =
if (now -. !last_update > delay_read) then
let n_connect =
n_connected workers_hash now
let n_connect =
n_connected workers_hash now
in
`Update n_connect
else if (now -. start_time > !stop_time) then
`Terminate
`Terminate
else if (now -. start_time > initialization_timeout) then
`Timeout
`Timeout
else
`None
`None
in
match (daemon, !status, test) with
| (_ , _ , `None ) -> ()
| (_ , Status.Running , `Terminate ) -> change_status Status.Stopping
| (false, Status.Running , `Update 0 ) -> change_status Status.Stopped
| (true , Status.Running , `Update 0 ) -> change_status Status.Queued
| (_ , _ , `Update i ) ->
| (true , Status.Running , `Update 0 ) -> change_status Status.Queued
| (_ , _ , `Update i ) ->
begin
status := Status.read ();
last_update := now;
stop_time := Input.Stop_time.(read () |> to_float) ;
let n_tot =
Hashtbl.length workers_hash
Hashtbl.length workers_hash
in
if (i <> n_tot) then
if (i <> n_tot) then
begin
Printf.sprintf "Connected workers : %d / %d" i n_tot
|> send_log "status" 0 now
|> send_log "status" 0 now
end
end
| (false, Status.Queued , `Timeout ) -> change_status Status.Stopped
| (false, Status.Queued , `Timeout ) -> change_status Status.Stopped
| (_, _, _) -> ()
;
done;
@ -487,37 +485,37 @@ let run ?(daemon=true) ezfio_filename =
Zmq.Socket.close socket
)
in
(** {3 Log thread} *)
let start_log_thread =
let start_log_thread =
let t0 =
Unix.gettimeofday ()
in
Thread.create (fun () ->
Thread.create (fun () ->
send_log "status" 0 t0 "Starting log thread";
let socket =
Zmq.Socket.create zmq_context Zmq.Socket.xsub
and address =
and address =
Printf.sprintf "tcp://*:%d" (port+3)
in
bind_socket "XSUB" socket address;
let pollitem =
Zmq.Poll.mask_of
Zmq.Poll.mask_of
[| (socket , Zmq.Poll.In) ;
(debug_socket , Zmq.Poll.In)
(debug_socket , Zmq.Poll.In)
|]
in
while (!status <> Status.Stopped)
do
let polling =
Zmq.Poll.poll ~timeout:1000 pollitem
Zmq.Poll.poll ~timeout:1000 pollitem
in
if (polling.(0) = Some Zmq.Poll.In) then
begin
let message =
let message =
Zmq.Socket.recv_all ~block:false socket
|> String.concat " "
in
@ -530,7 +528,7 @@ let run ?(daemon=true) ezfio_filename =
begin
(* Forward subscription from XPUB to XSUB *)
Zmq.Socket.recv_all ~block:false debug_socket
|> Zmq.Socket.send_all socket
|> Zmq.Socket.send_all socket
end
done;
Zmq.Socket.set_linger_period socket 1000 ;
@ -539,17 +537,17 @@ let run ?(daemon=true) ezfio_filename =
in
(** {3 Main thread} *)
let random_walkers n_walks =
let random_walkers n_walks =
let rec walkers accu = function
| 0 -> accu
| n ->
| 0 -> accu
| n ->
let random_int =
Random.int (Strictly_positive_int.to_int n_walks)
in
let new_accu =
let new_accu =
walkers_array.(random_int) :: accu
in
walkers new_accu (n-1)
walkers new_accu (n-1)
in
walkers [] (Strictly_positive_int.to_int n_walks)
|> Array.concat
@ -560,21 +558,21 @@ let run ?(daemon=true) ezfio_filename =
let wall0 =
Unix.gettimeofday ()
in
let f () =
let f () =
change_status Status.Queued;
send_log "status" 0 wall0 "Starting main thread";
(** Reply socket *)
let rep_socket =
Zmq.Socket.create zmq_context Zmq.Socket.rep
and address =
and address =
Printf.sprintf "tcp://*:%d" port
in
bind_socket "REP" rep_socket address;
Zmq.Socket.set_receive_high_water_mark rep_socket 100_000;
Zmq.Socket.set_send_high_water_mark rep_socket 100_000;
Zmq.Socket.set_immediate rep_socket true;
Zmq.Socket.set_immediate rep_socket true;
Zmq.Socket.set_linger_period rep_socket 600_000 ;
(** EZFIO Cache *)
@ -595,14 +593,14 @@ let run ?(daemon=true) ezfio_filename =
in
List.iter (fun x ->
if handle_ezfio ("has_"^x) = "T" then
try ignore @@ handle_ezfio ("get_"^x)
try ignore @@ handle_ezfio ("get_"^x)
with Failure _ -> ())
Qptypes.all_ezfio_messages;
(** Pull socket for computed data *)
let pull_socket =
Zmq.Socket.create zmq_context Zmq.Socket.pull
and address =
and address =
Printf.sprintf "tcp://*:%d" (port+2)
in
bind_socket "PULL" pull_socket address;
@ -611,7 +609,7 @@ let run ?(daemon=true) ezfio_filename =
(** Address of the dataserver *)
let server_address =
let ip =
Lazy.force Qmcchem_config.ip_address
Lazy.force Qmcchem_config.ip_address
in
Printf.sprintf "tcp://%s:%d" ip port
in
@ -621,7 +619,7 @@ let run ?(daemon=true) ezfio_filename =
(** Polling item to poll REP and PULL sockets. *)
let pollitem =
Zmq.Poll.mask_of
Zmq.Poll.mask_of
[| ( rep_socket, Zmq.Poll.In) ;
( pull_socket, Zmq.Poll.In) ;
|]
@ -629,27 +627,27 @@ let run ?(daemon=true) ezfio_filename =
(** Handles messages coming into the REP socket. *)
let handle_rep () =
let handle_rep () =
let raw_msg =
Zmq.Socket.recv_all ~block:false rep_socket
in
let t0 =
let t0 =
Unix.gettimeofday ()
in
let msg =
List.rev_map String.trim raw_msg
|> List.rev
|> Message.create
|> Message.create
and msg_size =
List.fold_left (fun accu x -> accu + (String.length x)) 0 raw_msg
in
let handle = function
| Message.Error _ -> ()
| Message.Ezfio ezfio_msg ->
let result =
| Message.Ezfio ezfio_msg ->
let result =
handle_ezfio ezfio_msg
in
Zmq.Socket.send_all rep_socket
Zmq.Socket.send_all rep_socket
[ String.length result
|> Printf.sprintf "%d " ;
result ] ;
@ -657,53 +655,53 @@ let run ?(daemon=true) ezfio_filename =
| Message.GetWalkers n_walks ->
begin
send_log "req" msg_size t0 "get_walkers";
let result =
let result =
random_walkers n_walks
in
Zmq.Socket.send_all rep_socket result;
send_log "rep" walkers_size t0 "get_walkers"
end
| Message.Register (w,pid) ->
end
| Message.Register (w,pid) ->
begin
match !status with
| Status.Queued
| Status.Queued
| Status.Running ->
begin
String.concat " " [ "Register :" ;
Compute_node.to_string w ;
String.concat " " [ "Register :" ;
Compute_node.to_string w ;
string_of_int pid ]
|> send_log "req" msg_size t0;
add_worker w pid;
if (!status = Status.Queued) then
change_status Status.Running ;
Zmq.Socket.send rep_socket "OK";
send_log "rep" 2 t0 "Register : OK"
send_log "rep" 2 t0 "Register : OK"
end
| Status.Stopping
| Status.Stopping
| Status.Stopped ->
Zmq.Socket.send rep_socket "Failed";
end
| Message.Unregister (w,pid) ->
| Message.Unregister (w,pid) ->
begin
String.concat " " [ "Unregister :" ;
(Compute_node.to_string w) ;
(Compute_node.to_string w) ;
(string_of_int pid) ]
|> send_log "req" msg_size t0;
Zmq.Socket.send rep_socket "OK";
del_worker w pid;
String.concat " " [ "Unregister :";
String.concat " " [ "Unregister :";
(Hashtbl.length workers_hash) |> string_of_int ;
"remaining" ]
"remaining" ]
|> send_log "rep" 2 t0 ;
let n_connect =
let n_connect =
n_connected workers_hash t0
in
match (daemon,n_connect) with
match (daemon,n_connect) with
| (false,0) -> change_status Status.Stopped
| (true ,0) -> change_status Status.Queued
| (true ,0) -> change_status Status.Queued
| _ -> ()
end
| Message.Test ->
| Message.Test ->
begin
Zmq.Socket.send rep_socket "OK";
send_log "rep" 2 t0 "Test"
@ -719,17 +717,17 @@ let run ?(daemon=true) ezfio_filename =
let raw_msg =
Zmq.Socket.recv_all ~block:false pull_socket
in
let t0 =
let t0 =
Unix.gettimeofday ()
in
let msg =
List.rev_map String.trim raw_msg
|> List.rev
|> Message.create
|> Message.create
and msg_size =
List.fold_left (fun accu x -> accu + (String.length x)) 0 raw_msg
in
let recv_log =
let recv_log =
send_log "pull" msg_size t0
in
@ -739,7 +737,7 @@ let run ?(daemon=true) ezfio_filename =
begin
if (status = Status.Running) then
touch_worker h pid ;
let log_msg =
let log_msg =
Printf.sprintf "Walkers from %s : %d / %d / %d"
(key h pid) (Array.length w) (!last_walker) walk_num_tot
in
@ -754,13 +752,14 @@ let run ?(daemon=true) ezfio_filename =
(Unix.gettimeofday () -. wall0)
1. (Property.to_string Property.Wall)
hostname (string_of_int dataserver_pid) 1
|> Block.of_string
|> Block.of_string
in
match wall with
| Some wall ->
| Some wall ->
begin
output_string !block_channel (Block.to_string wall);
output_char !block_channel '\n';
output_string !block_channel (Block.to_string_or_bytes wall);
if not Qmcchem_config.binary_io then
output_char !block_channel '\n';
end
| _ -> ()
end
@ -768,16 +767,17 @@ let run ?(daemon=true) ezfio_filename =
begin
if (status = Status.Running) then
touch_worker b.Block.compute_node b.Block.pid ;
output_string !block_channel (Block.to_string b);
output_char !block_channel '\n';
recv_log (Block.to_string b)
output_string !block_channel (Block.to_string_or_bytes b);
if not Qmcchem_config.binary_io then
output_char !block_channel '\n';
recv_log (Block.to_short_string b)
end
| Message.Test
| Message.GetWalkers _
| Message.Ezfio _
| Message.Register (_, _)
| Message.Unregister (_, _)
-> failwith "Bad message"
-> failwith "Bad message"
in handle msg
in
@ -785,18 +785,18 @@ let run ?(daemon=true) ezfio_filename =
while (!status <> Status.Stopped)
do
let polling =
Zmq.Poll.poll ~timeout:1000 pollitem
Zmq.Poll.poll ~timeout:1000 pollitem
in
match polling.(1) with
| Some Zmq.Poll.In -> handle_pull !status
| _ ->
| _ ->
begin
match polling.(0) with
| Some Zmq.Poll.In -> handle_rep ()
| _ ->
| _ ->
begin
flush !block_channel ;
let file_size =
let file_size =
(Unix.stat block_channel_filename_locked).Unix.st_size
in
if (file_size > !max_file_size) then
@ -811,13 +811,13 @@ let run ?(daemon=true) ezfio_filename =
List.iter (fun socket ->
Zmq.Socket.set_linger_period socket 1000 ;
Zmq.Socket.close socket)
[ rep_socket ; pull_socket ]
[ rep_socket ; pull_socket ]
in
Thread.create f
in
(** {2 Finalization} *)
(** Cleans all the open files, sockets, etc.
@ -842,16 +842,16 @@ let run ?(daemon=true) ezfio_filename =
(** {3 Main function} *)
let t0 =
let t0 =
Unix.gettimeofday ()
in
(* Handle signals *)
let handler s =
let handler s =
Printf.printf "Dataserver received signal %d... killing\n%!" s;
Watchdog.kill ();
in
List.iter (fun s -> ignore @@ Sys.signal s (Sys.Signal_handle handler))
List.iter (fun s -> ignore @@ Sys.signal s (Sys.Signal_handle handler))
[
Sys.sigint ;
Sys.sigterm ;
@ -864,7 +864,7 @@ let run ?(daemon=true) ezfio_filename =
begin
try
(List.iter Thread.join
(List.iter Thread.join
[ start_status_thread () ;
start_log_thread () ;
start_main_thread () ;

265
ocaml/Random_variable.ml
View File

@ -1,13 +1,13 @@
open Qptypes
type t =
type t =
{ property : Property.t ;
data : Block.t list;
}
module Average = struct
include Sample
include Sample
end
module Error = struct
@ -19,42 +19,42 @@ module Variance = struct
end
module Skewness: sig
type t
type t
val to_float : t -> float
val of_float : float -> t
val to_string : t -> string
end = struct
type t = float
type t = float
let to_string = string_of_float
let to_float x = x
let of_float x = x
end
module Kurtosis: sig
type t
type t
val to_float : t -> float
val of_float : float -> t
val to_string : t -> string
end = struct
type t = float
type t = float
let to_string = string_of_float
let to_float x = x
let of_float x = x
end
module GaussianDist: sig
type t
type t
val create : mu:Average.t -> sigma2:Variance.t -> t
val eval : g:t -> x:float -> float
end = struct
type t = { mu: Average.t ; sigma2: Variance.t }
type t = { mu: Average.t ; sigma2: Variance.t }
let create ~mu ~sigma2 =
{ mu ; sigma2 }
let eval ~g ~x =
let { mu ; sigma2 } =
let { mu ; sigma2 } =
g
in
let mu =
let mu =
Average.to_float mu
and sigma2 =
Variance.to_float sigma2
@ -62,10 +62,10 @@ end = struct
let x2 =
(x -. mu) *. ( x -. mu) /. sigma2
in
let pi =
let pi =
acos (-1.)
in
let c =
let c =
1. /. (sqrt (sigma2 *. (pi +. pi)))
in
c *. exp ( -0.5 *. x2)
@ -78,7 +78,7 @@ let hashtbl_to_alist table =
let hashtbl_change table key f =
let elt =
try
try
Some (Hashtbl.find table key)
with
| Not_found -> None
@ -91,7 +91,7 @@ let hashtbl_change table key f =
(** Build from raw data. Range values are given in percent. *)
let of_raw_data ?(locked=true) ~range property =
let data =
let data =
Block.raw_data ~locked ()
|> List.filter (fun x -> x.Block.property = property)
|> List.sort (fun x y ->
@ -109,7 +109,7 @@ let of_raw_data ?(locked=true) ~range property =
(Weight.to_float x.Block.weight) +. accu
) 0. data
in
let wmin, wmax =
rmin *. total_weight *. 0.01,
rmax *. total_weight *. 0.01
@ -128,13 +128,13 @@ let of_raw_data ?(locked=true) ~range property =
(wsum_new, x::l)
else
(wsum_new, l)
end
end
) (0.,[]) data
in
List.rev new_data
in
let result =
let result =
match range with
| (0.,100.) -> { property ; data }
| (rmin,rmax) -> { property ; data=data_in_range rmin rmax }
@ -146,7 +146,7 @@ let of_raw_data ?(locked=true) ~range property =
(** Compute average *)
let average { property ; data } =
if Property.is_scalar property then
let (num,denom) =
let (num,denom) =
List.fold_left (fun (an, ad) x ->
let num =
(Weight.to_float x.Block.weight) *. (Sample.to_float x.Block.value)
@ -154,7 +154,7 @@ let average { property ; data } =
(Weight.to_float x.Block.weight)
in (an +. num, ad +. den)
) (0., 0.) data
in
in
num /. denom
|> Average.of_float
else
@ -163,15 +163,15 @@ let average { property ; data } =
| [] -> 1
| x :: tl -> Sample.dimension x.Block.value
in
let (num,denom) =
let (num,denom) =
List.fold_left (fun (an, ad) x ->
let num =
Array.map (fun y -> (Weight.to_float x.Block.weight) *. y)
(Sample.to_float_array x.Block.value)
(Sample.to_float_array x.Block.value)
and den = (Weight.to_float x.Block.weight)
in ( Array.mapi (fun i y -> y +. num.(i)) an , ad +. den)
in ( Array.mapi (fun i y -> y +. num.(i)) an , ad +. den)
) (Array.make dim 0. , 0.) data
in
in
let denom_inv =
1. /. denom
in
@ -180,22 +180,28 @@ let average { property ; data } =
(** Compute sum (for CPU/Wall time) *)
let sum { property ; data } =
List.fold_left (fun accu x ->
let num = (Weight.to_float x.Block.weight) *. (Sample.to_float x.Block.value)
in accu +. num
) 0. data
) 0. data
(** Calculation of the average and error bar *)
let ave_error { property ; data } =
let rec loop ~sum ~avsq ~ansum ~avsum ~n ?idx = function
| [] ->
(* sum: \sum_k x_k *. w_k
ansum: \sum_k w_k
avsum: \sum_k x_k *. w_k
avcu0: avsum / ansum
avsq: \sum_k (1. -. (w_k /. ansum_k)) *. (x_k -. avcu0)^2 *. w_k)
*)
let rec loop ~sum ~avsq ~ansum ~avsum ~n ?idx = function
| [] ->
begin
if (n > 0.) then
( Average.of_float (sum /. ansum),
@ -205,12 +211,8 @@ let ave_error { property ; data } =
end
| (x,w) :: tail ->
begin
let avcu0 =
avsum /. ansum
in
let xw =
x *. w
in
let avcu0 = avsum /. ansum in
let xw = x *. w in
let ansum, avsum, sum =
ansum +. w ,
avsum +. xw ,
@ -220,9 +222,9 @@ let ave_error { property ; data } =
~sum:sum
~avsq:(avsq +. (1. -. (w /. ansum)) *. (x -. avcu0)
*. (x -. avcu0) *. w)
~avsum:avsum
~avsum:avsum
~ansum:ansum
~n:(n +. 1.)
~n:(n +. 1.)
end
in
@ -242,7 +244,7 @@ let ave_error { property ; data } =
(Sample.to_float x.Block.value,
Weight.to_float x.Block.weight)
) data
|> ave_error_scalar
|> ave_error_scalar
else
match data with
| [] -> (Average.of_float 0., None)
@ -251,7 +253,7 @@ let ave_error { property ; data } =
head.Block.value
|> Sample.dimension
in
let result =
let result =
Array.init dim (fun idx ->
List.rev_map (fun x ->
(Sample.to_float ~idx x.Block.value,
@ -260,16 +262,16 @@ let ave_error { property ; data } =
|> ave_error_scalar
)
in
( Array.map (fun (x,_) -> Average.to_float x) result
|> Average.of_float_array ~dim ,
( Array.map (fun (x,_) -> Average.to_float x) result
|> Average.of_float_array ~dim ,
if (Array.length result < 2) then
None
else
Some (Array.map (function
Some (Array.map (function
| (_,Some y) -> Error.to_float y
| (_,None) -> 0.) result
| (_,None) -> 0.) result
|> Average.of_float_array ~dim)
)
)
@ -286,14 +288,14 @@ let fold_blocks ~f { property ; data } =
List.fold_left (fun accu block ->
let x = Sample.to_float block.Block.value
in f accu x
) init data
) init data
(** Convergence plot *)
let convergence { property ; data } =
let rec loop ~sum ~avsq ~ansum ~avsum ~n ~accu = function
let rec loop ~sum ~avsq ~ansum ~avsum ~n ~accu = function
| [] -> List.rev accu
| head :: tail ->
begin
@ -307,7 +309,7 @@ let convergence { property ; data } =
and avsum = avsum +. xw
and sum = sum +. xw
in
let accu =
let accu =
if (n > 0.) then
(sum /. ansum, sqrt ( abs_float ( avsq /.( ansum *. n))))::accu
else
@ -317,9 +319,9 @@ let convergence { property ; data } =
~sum:sum
~avsq:(avsq +. (1. -. (w /. ansum)) *. (x -. avcu0)
*. (x -. avcu0) *. w)
~avsum:avsum
~avsum:avsum
~ansum:ansum
~n:(n +. 1.)
~n:(n +. 1.)
~accu:accu
end
in
@ -365,7 +367,7 @@ let max_block =
(** Create a hash table for merging *)
let create_hash ~create_key ?(update_block_id=(fun x->x))
let create_hash ~create_key ?(update_block_id=(fun x->x))
?(update_value=(fun wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw) )
?(update_weight=(fun wc wb -> wc +. wb) ) t =
let table = Hashtbl.create 63
@ -374,7 +376,7 @@ let create_hash ~create_key ?(update_block_id=(fun x->x))
let key = create_key block
in
let open Block in
hashtbl_change table key (function
hashtbl_change table key (function
| Some current ->
let wc, wb =
Weight.to_float current.weight,
@ -384,7 +386,7 @@ let create_hash ~create_key ?(update_block_id=(fun x->x))
update_weight wc wb
in
if (Property.is_scalar current.property) then
let vc, vb =
let vc, vb =
Sample.to_float current.value,
Sample.to_float block.value
in Some
@ -396,15 +398,15 @@ let create_hash ~create_key ?(update_block_id=(fun x->x))
compute_node = block.compute_node;
}
else
let vc, vb =
let vc, vb =
Sample.to_float_array current.value,
Sample.to_float_array block.value
and dim =
and dim =
Sample.dimension current.value
in Some
{ property = current.property ;
weight = Weight.of_float sw ;
value =
value =
Array.init dim (fun i -> update_value wc vc.(i) wb vb.(i) sw)
|> Sample.of_float_array ~dim ;
block_id = update_block_id block.block_id;
@ -443,15 +445,15 @@ let merge ~create_key ?update_block_id ?update_value ?update_weight t =
(** Merge per block id *)
let merge_per_block_id =
merge
merge
~create_key:(fun block -> Block_id.to_string block.Block.block_id)
(** Merge per compute_node *)
let merge_per_compute_node =
merge
~create_key:(fun block ->
Printf.sprintf "%s"
~create_key:(fun block ->
Printf.sprintf "%s"
(Compute_node.to_string block.Block.compute_node) )
@ -459,8 +461,8 @@ let merge_per_compute_node =
(** Merge per Compute_node and PID *)
let merge_per_compute_node_and_pid =
merge
~create_key:(fun block ->
Printf.sprintf "%s %10.10d"
~create_key:(fun block ->
Printf.sprintf "%s %10.10d"
(Compute_node.to_string block.Block.compute_node)
(block.Block.pid) )
@ -469,8 +471,8 @@ let merge_per_compute_node_and_pid =
(** Merge per Compute_node and BlockId *)
let merge_per_compute_node_and_block_id =
merge
~create_key:(fun block ->
Printf.sprintf "%s %10.10d"
~create_key:(fun block ->
Printf.sprintf "%s %10.10d"
(Compute_node.to_string block.Block.compute_node)
(Block_id.to_int block.Block.block_id) )
@ -510,48 +512,48 @@ let error_x_over_y = function
(** Create float, variable operators *)
let one_variable_operator ~update_value p f =
{ p with
data = List.rev @@ List.rev_map (fun b -> { b with
let one_variable_operator ~update_value p f =
{ p with
data = List.rev @@ List.rev_map (fun b -> { b with
Block.value = Sample.of_float (update_value (Sample.to_float b.Block.value) ) }
) p.data }
let ( +@ ) p f = one_variable_operator p f
~update_value: (fun x -> x +. f )
let ( *@ ) p f = one_variable_operator p f
let ( *@ ) p f = one_variable_operator p f
~update_value: (fun x -> x *. f )
let ( -@ ) p f = one_variable_operator p f
let ( -@ ) p f = one_variable_operator p f
~update_value: (fun x -> x -. f )
let ( /@ ) p f = one_variable_operator p f
let ( /@ ) p f = one_variable_operator p f
~update_value: (fun x -> x /. f )
(** Create two variable operators *)
let two_variable_operator ~update_value p1 p2 =
let two_variable_operator ~update_value p1 p2 =
merge
~update_value
~create_key:(fun block ->
Printf.sprintf "%s %10.10d %10.10d"
Printf.sprintf "%s %10.10d %10.10d"
(Compute_node.to_string block.Block.compute_node)
(Block_id.to_int block.Block.block_id)
(Block_id.to_int block.Block.block_id)
(block.Block.pid) )
~update_weight:(fun wc wb -> wc )
{ property = p1.property ;
data = List.concat [ p1.data ; p2.data ] }
data = List.concat [ p1.data ; p2.data ] }
let ( +! ) = two_variable_operator
let ( +! ) = two_variable_operator
~update_value: (fun wc vc wb vb sw -> (vc +. vb) )
let ( *! ) = two_variable_operator
let ( *! ) = two_variable_operator
~update_value: (fun wc vc wb vb sw -> (vc *. vb) )
let ( -! ) = two_variable_operator
let ( -! ) = two_variable_operator
~update_value: (fun wc vc wb vb sw -> (vc -. vb) )
let ( /! ) = two_variable_operator
let ( /! ) = two_variable_operator
~update_value: (fun wc vc wb vb sw -> (vc /. vb) )
@ -560,11 +562,11 @@ let ( /! ) = two_variable_operator
(** Merge two consecutive blocks *)
let compress =
merge
~create_key:(fun block ->
~create_key:(fun block ->
Printf.sprintf "%s %10.10d %10.10d"
(Compute_node.to_string block.Block.compute_node) block.Block.pid
(((Block_id.to_int block.Block.block_id)+1)/2))
~update_block_id:(fun block_id ->
~update_block_id:(fun block_id ->
((Block_id.to_int block_id)+1)/2
|> Block_id.of_int )
@ -576,15 +578,15 @@ let max_value_per_compute_node t =
let table = Hashtbl.create 63
in
let create_key block =
Printf.sprintf "%s %10.10d"
Printf.sprintf "%s %10.10d"
(Compute_node.to_string block.Block.compute_node)
(block.Block.pid)
(block.Block.pid)
in
List.iter (fun block ->
let key = create_key block
in
let open Block in
hashtbl_change table key (function
hashtbl_change table key (function
| Some current ->
let vc = Sample.to_float current.value
and vb = Sample.to_float block.value
@ -610,36 +612,36 @@ let max_value_per_compute_node t =
(** String representation *)
let to_string p =
let to_string p =
match p.property with
| Property.Cpu -> Printf.sprintf "%s" (Time.string_of_sec (sum p))
| Property.Wall -> Printf.sprintf "%s" (Time.string_of_sec (sum (max_value_per_compute_node p)))
| Property.Accep -> Printf.sprintf "%16.10f" (average p |> Average.to_float)
| _ ->
| _ ->
begin
if Property.is_scalar p.property then
match ave_error p with
| (ave, Some error) ->
let (ave, error) =
Average.to_float ave,
| (ave, Some error) ->
let (ave, error) =
Average.to_float ave,
Error.to_float error
in
Printf.sprintf "%16.10f +/- %16.10f" ave error
| (ave, None) ->
let ave =
| (ave, None) ->
let ave =
Average.to_float ave
in
Printf.sprintf "%16.10f" ave
else
match ave_error p with
| (ave, Some error) ->
| (ave, Some error) ->
let idxmax =
Average.dimension ave
in
let rec f accu idx =
if (idx < idxmax) then
let (ave, error) =
Average.to_float ~idx ave,
let (ave, error) =
Average.to_float ~idx ave,
Error.to_float ~idx error
in
let s =
@ -650,9 +652,9 @@ let to_string p =
accu
in
(f "[ \n" 0) ^ " ]"
| (ave, None) ->
| (ave, None) ->
Average.to_float ave
|> Printf.sprintf "%16.10f"
|> Printf.sprintf "%16.10f"
end
@ -666,19 +668,19 @@ let compress_files () =
let properties =
Lazy.force Block.properties
in
(* Create temporary file *)
let dir_name =
Block.dir_name
in
let dir_name =
let dir_name =
Lazy.force dir_name
in
let files =
let files =
Sys.readdir dir_name
|> Array.to_list
|> List.filter (fun x ->
|> List.filter (fun x ->
try
Str.search_backward (Str.regexp "locked") x (String.length x) >= 0
with
@ -688,10 +690,10 @@ let compress_files () =
|> List.rev
in
let out_channel_dir =
let out_channel_dir =
let rand_num = Random.int 1000000 |> string_of_int in
let dirname =
Filename.concat !Ezfio.ezfio_filename "blocks"
let dirname =
Filename.concat !Ezfio.ezfio_filename "blocks"
in
if not ( Sys.file_exists dirname ) then
Unix.mkdir dirname 0o755;
@ -706,31 +708,31 @@ let compress_files () =
raise (Sys_error message)
in
let out_channel_name =
let out_channel_name =
let hostname =
Lazy.force Qmcchem_config.hostname
Lazy.force Qmcchem_config.hostname
and suffix =
Unix.getpid ()
|> string_of_int
|> string_of_int
in
String.concat "." [ hostname ; suffix ]
in
let block_channel =
let block_channel =
Filename.concat out_channel_dir out_channel_name
|> open_out
in
List.iter (fun p ->
let l =
let l =
match p with
| Property.Cpu
| Property.Cpu
| Property.Accep ->
of_raw_data ~locked:false ~range:(0.,100.) p
|> merge_per_compute_node
| Property.Wall ->
| Property.Wall ->
of_raw_data ~locked:false ~range:(0.,100.) p
|> max_value_per_compute_node
|> max_value_per_compute_node
| _ ->
of_raw_data ~locked:false ~range:(0.,100.) p
(*
@ -738,10 +740,11 @@ let compress_files () =
*)
in
List.iter (fun x ->
output_string block_channel (Block.to_string x);
output_char block_channel '\n';
) l.data
) properties ;
output_string block_channel (Block.to_string_or_bytes x);
if not Qmcchem_config.binary_io then
output_char block_channel '\n';
) l.data
) properties ;
close_out block_channel;
List.iter Unix.unlink files ;
@ -760,17 +763,17 @@ let autocovariance { property ; data } =
match (merge_per_block_id { property ; data })
with { property ; data } -> Array.of_list data
in
let x_t =
let x_t =
Array.map (fun x -> (Sample.to_float x.Block.value) -. ave) data
in
let f i =
let denom =
let f i =
let denom =
if (i > 1) then (float_of_int i) else 1.
in
let r =
let r =
Array.sub x_t 0 i
|> Array.fold_left (fun accu x ->
accu +. x *. x_t.(i)) 0.
accu +. x *. x_t.(i)) 0.
in
r /. denom
in
@ -786,14 +789,14 @@ let centered_cumulants { property ; data } =
|> Average.to_float
in
let centered_data =
List.rev_map (fun x ->
( (Weight.to_float x.Block.weight),
List.rev_map (fun x ->
( (Weight.to_float x.Block.weight),
(Sample.to_float x.Block.value) -. ave )
) data
|> List.rev
in
let var =
let (num, denom) =
let var =
let (num, denom) =
List.fold_left (fun (a2, ad) (w,x) ->
let x2 = x *. x
in
@ -802,18 +805,18 @@ let centered_cumulants { property ; data } =
in (a2 +. var, ad +. den)
) (0., 0.) centered_data
in num /. denom
in
in
let centered_data =
let sigma_inv =
1. /. (sqrt var)
in
List.rev_map (fun x ->
( (Weight.to_float x.Block.weight),
List.rev_map (fun x ->
( (Weight.to_float x.Block.weight),
( (Sample.to_float x.Block.value) -. ave ) *. sigma_inv )
) data
|> List.rev
in
let (cum3,cum4) =
let (cum3,cum4) =
let (cum3, cum4, denom) =
List.fold_left (fun (a3, a4, ad) (w,x) ->
let x2 = x *. x
@ -823,9 +826,9 @@ let centered_cumulants { property ; data } =
and den = w
in (a3 +. cum3, a4 +. cum4, ad +. den)
) (0., 0., 0.) centered_data
in
in
( cum3 /. denom, cum4 /. denom -. 3. )
in
in
[| ave ; var ; cum3 ; cum4 |]
@ -833,26 +836,26 @@ let centered_cumulants { property ; data } =
(** Computes a histogram *)
let histogram { property ; data } =
let min, max =
let min, max =
(min_block { property ; data }),
(max_block { property ; data })
in
let length =
let length =
max -. min
and n =
List.length data
|> float_of_int
|> sqrt
in
let delta_x =
let delta_x =
length /. (n-.1.)
and result =
Array.init (int_of_float n + 1) (fun _ -> 0.)
Array.init (int_of_float n + 1) (fun _ -> 0.)
in
List.iter (fun x ->
let w =
(Weight.to_float x.Block.weight)
and x =
and x =
(Sample.to_float x.Block.value)
in
let i =
@ -862,7 +865,7 @@ let histogram { property ; data } =
result.(i) <- result.(i) +. w
) data
;
let norm =
let norm =
1. /. ( delta_x *. (
Array.fold_left (fun accu x -> accu +. x) 0. result
) )

39
ocaml/Sample.ml
View File

@ -1,6 +1,6 @@
open Sexplib.Std
type t =
type t =
| One_dimensional of float
| Multidimensional of (float array * int)
[@@deriving sexp]
@ -9,27 +9,27 @@ let dimension = function
| One_dimensional _ -> 1
| Multidimensional (_,d) -> d
let to_float ?idx x =
let to_float ?idx x =
match (idx,x) with
| None , One_dimensional x
| Some 0, One_dimensional x -> x
| Some 0, One_dimensional x -> x
| Some i, One_dimensional x ->
failwith "Index should not be specified in One_dimensional"
| None , Multidimensional (x,_) -> x.(0)
| Some i, Multidimensional (x,s) when i < s -> x.(i)
| Some i, Multidimensional (x,s) ->
Printf.sprintf "Index out of bounds in Multidimensional
| Some i, Multidimensional (x,s) ->
Printf.sprintf "Index out of bounds in Multidimensional
%d not in [0,%d[ " i s
|> failwith
let to_float_array = function
let to_float_array = function
| One_dimensional _ -> failwith "Should be Multidimensional"
| Multidimensional (x,_) -> x
let of_float x =
let of_float x =
One_dimensional x
let of_float_array ~dim x =
let of_float_array ~dim x =
if (Array.length x) <> dim then
failwith "Inconsistent array size in of_float_array"
else
@ -39,9 +39,26 @@ let of_float_array ~dim x =
let to_string = function
| One_dimensional x -> string_of_float x
| Multidimensional (x,_) ->
| Multidimensional (x,_) ->
Array.map string_of_float x
|> Array.to_list
|> String.concat " "
|> Printf.sprintf "%s"
|> String.concat " "
let to_bytes = function
| One_dimensional x -> Qptypes.bytes_of_float x
| Multidimensional (x,_) ->
let b = Bytes.create (8 * Array.length x) in
Array.iteri (fun i x ->
Int64.bits_of_float x
|> Bytes.set_int64_le b (i*8) ) x;
b
let of_bytes b =
match Bytes.length b with
| 8 -> let x = Qptypes.float_of_bytes b in
One_dimensional x
| l -> let len = l/8 in
Multidimensional ( Array.init len (fun i ->
Bytes.get_int64_le b (i*8)
|> Int64.float_of_bits ),
len )

4
ocaml/Sample.mli
View File

@ -2,7 +2,9 @@ type t [@@deriving sexp]
val to_float : ?idx:int -> t -> float
val to_float_array : t -> float array
val of_float : float -> t
val of_float_array : dim:int -> float array -> t
val of_float_array : dim:int -> float array -> t
val to_string : t -> string
val to_bytes : t -> bytes
val of_bytes : bytes -> t
val dimension : t -> int

88
ocaml/qptypes_generator.ml
View File

@ -4,7 +4,15 @@ let global_replace x =
|> Str.global_replace (Str.regexp "Float.of_string") "float_of_string"
|> Str.global_replace (Str.regexp "Int.to_string") "string_of_int"
|> Str.global_replace (Str.regexp "Int.of_string") "int_of_string"
|> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") ""
|> Str.global_replace (Str.regexp "Int.to_bytes") "bytes_of_int"
|> Str.global_replace (Str.regexp "Int64.to_bytes") "bytes_of_int64"
|> Str.global_replace (Str.regexp "Float.to_bytes") "bytes_of_float"
|> Str.global_replace (Str.regexp "Float.of_bytes") "float_of_bytes"
|> Str.global_replace (Str.regexp "Int.of_bytes") "int_of_bytes"
|> Str.global_replace (Str.regexp "Int64.of_bytes") "int64_of_bytes"
|> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") ""
|> Str.global_replace (Str.regexp "String.to_bytes") "Bytes.of_string"
|> Str.global_replace (Str.regexp "String.of_bytes") "Bytes.to_string"
let input_data = "
* Positive_float : float
@ -38,7 +46,7 @@ let input_data = "
* Negative_int : int
if not (x <= 0) then
raise (Invalid_argument (Printf.sprintf \"Negative_int : (x <= 0) : x=%d\" x));
assert (x <= 0) ;
assert (x <= 0) ;
* Det_coef : float
if (x < -1.) || (x > 1.) then
@ -167,6 +175,34 @@ let input_ezfio = "
let untouched = "
let bytes_of_int64 i =
let result = Bytes.create 8 in
Bytes.set_int64_le result 0 i;
result
let bytes_of_int i =
Int64.of_int i
|> bytes_of_int64
let int64_of_bytes b =
Bytes.get_int64_le b 0
let int_of_bytes b =
int64_of_bytes b
|> Int64.to_int
let float_of_bytes b =
int64_of_bytes b
|> Int64.float_of_bits
let bytes_of_float f =
Int64.bits_of_float f
|> bytes_of_int64
"
let template = format_of_string "
@ -175,11 +211,15 @@ module %s : sig
val to_%s : t -> %s
val of_%s : %s %s -> t
val to_string : t -> string
val to_bytes : t -> bytes
val of_bytes : bytes -> t
end = struct
type t = %s [@@deriving sexp]
let to_%s x = x
let of_%s %s x = ( %s x )
let to_string x = %s.to_string x
let to_bytes x = %s.to_bytes x
let of_bytes b = %s.of_bytes b
end
"
@ -203,7 +243,7 @@ let parse_input input=
and name = String_ext.strip name in
let typ_cap = String.capitalize_ascii typ in
let newstring = Printf.sprintf template name typ typ typ params_val typ typ
typ typ params ( String_ext.strip text ) typ_cap
typ typ params ( String_ext.strip text ) typ_cap typ_cap typ_cap
in
List.rev (parse (newstring::result) tail )
in
@ -223,9 +263,11 @@ module %s : sig
val get_max : unit -> %s
val of_%s : ?min:%s -> ?max:%s -> %s -> t
val to_string : t -> string
val to_bytes : t -> bytes
end = struct
type t = %s [@@deriving sexp]
let to_string x = %s.to_string x
let to_bytes x = %s.to_bytes x
let get_max () =
if (Ezfio.has_%s ()) then
Ezfio.get_%s ()
@ -251,6 +293,10 @@ end = struct
end
"
(*
val of_bytes : bytes -> t
let of_bytes x = %s.of_bytes x
*)
let parse_input_ezfio input=
let parse s =
@ -271,9 +317,10 @@ let parse_input_ezfio input=
| [ name ; typ ; ezfio_func ; min ; max ; msg ] -> (name, typ, ezfio_func, min, max, msg)
| _ -> assert false
in
let typ_cap = String.capitalize_ascii typ in
Printf.sprintf ezfio_template
name typ typ typ typ typ typ typ typ (String.capitalize_ascii typ)
ezfio_func ezfio_func max min typ typ max msg min name (String.capitalize_ascii typ)
name typ typ typ typ typ typ typ typ typ_cap typ_cap
ezfio_func ezfio_func max min typ typ max msg min name typ_cap
end
| _ -> failwith "Error in input_ezfio"
in
@ -294,29 +341,30 @@ let input_lines filename =
let create_ezfio_handler () =
let lines =
let lines =
input_lines "ezfio.ml"
|> List.mapi (fun i l -> if i > 417 then Some l else None)
(* /!\ Change when ezfio.ml changes *)
|> List.mapi (fun i l -> if i > 442 then Some l else None)
|> List.filter (fun x -> x <> None)
|> List.map (fun x ->
match x with
| Some x -> x
| None -> assert false)
in
let functions =
let functions =
List.map (fun x ->
match String.split_on_char ' ' x with
| _ :: x :: "()" :: "=" :: f :: dir :: item :: _-> (x, f, dir, item)
| _ :: x :: "=" :: f :: dir :: item :: _-> (x, f, dir, item)
| _ -> ("","","","")
) lines
) lines
in
let has_functions =
List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "has_") functions
and get_functions =
let has_functions =
List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "has_") functions
and get_functions =
List.filter (fun (x,_,_,_) -> String.sub x 0 4 = "get_") functions
in
let chop s =
let chop s =
match (Str.split_delim (Str.regexp ";;") s) with
| x :: _ -> x
| _ -> assert false
@ -329,17 +377,17 @@ match msg with " ] @
List.map (fun (x,f,d,i) ->
let i = chop i in
if (String.sub f ((String.length f)-6) 6 = "_array") then
Printf.sprintf " | \"%s\" ->
Printf.sprintf " | \"%s\" ->
Ezfio.read_string_array %s %s
|> Ezfio.flattened_ezfio
|> Ezfio.flattened_ezfio
|> Array.to_list
|> String.concat \" \"" x d i
else
Printf.sprintf " | \"%s\" -> Ezfio.read_string %s %s" x d i
) get_functions
) get_functions
) @ (
List.map (fun (x,_,_,_) ->
Printf.sprintf " | \"%s\" -> if (Ezfio.%s ()) then \"T\" else \"F\"" x x
Printf.sprintf " | \"%s\" -> if (Ezfio.%s ()) then \"T\" else \"F\"" x x
) has_functions
)
@ [" | x -> failwith (x^\" : Unknown EZFIO function\")\n;;" ;
@ -347,17 +395,17 @@ match msg with " ] @
(
List.rev_map (fun (x,_,_,_) ->
Printf.sprintf " \"%s\" ; " (String.sub x 4 ((String.length x)-4))
) has_functions
) has_functions
) @ ["]"]
in
String.concat "\n" result
String.concat "\n" result
|> print_endline
(** Main *)
let () =
print_endline untouched;
parse_input input_data ;
parse_input_ezfio input_ezfio;
print_endline untouched;
create_ezfio_handler ()

6
scripts/compile_src.sh
View File

@ -18,12 +18,6 @@ OBJ="${OBJ} IRPF90_temp/ZMQ/f77_zmq_module.o"
INCLUDES="${INCLUDES} -I AO -I SAMPLING -I TOOLS -I JASTROW -I MAIN -I PROPERTIES -I ZMQ"
IRPF90_FLAGS="${IRPF90_FLAGS} ${INCLUDES}"
# Check IRPF90 version
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
fi
export IRPF90 IRPF90_FLAGS INCLUDES LIB SRC OBJ

8
scripts/create_properties_ezfio.py
View File

@ -145,6 +145,14 @@ let to_string = function
for p in properties_qmcvar:
print >>file, """| %(P)s -> "%(P)s" """%{'P':p[1].capitalize(), 'p':p[1]}
print >>file, """;;
let of_bytes x =
Bytes.to_string x
|> of_string
let to_bytes x =
to_string x
|> Bytes.of_string
"""
# is_scalar

277
src/MAIN/admc.org Normal file
View File

@ -0,0 +1,277 @@
#+TITLE: Asynchronous DMC
#+AUTHOR: Anthony Scemama
#+EMAIL: scemama@irsamc.ups-tlse.fr
#+PROPERTY: header-args :tangle no :noweb yes
* Main program
** Declarations
#+NAME: declarations
#+begin_src f90
include '../types.F'
integer :: iter
integer :: k_full ! Index of walkers in elec_coord_full
integer :: k_new ! Index of walkers in elec_coord_new
integer :: iw ! Number of copies in branching
integer :: l
real, allocatable :: elec_coord_new(:,:,:)
double precision :: w
double precision, allocatable :: E_out(:), w_sum(:)
double precision, external :: qmc_ranf
allocate(elec_coord_new(elec_num+1,3,walk_num))
allocate(E_out(walk_num), w_sum(walk_num))
#+end_src
** Main flow
- Fetch ~walk_num~ electron coordinates in ~elec_coord_full~
- For each set of coordinates,
- Make a PDMC trajectory, and output the weight ~w~
- Perform branching depending on the value of the weight
- Store the new sets of coordinates in ~elec_coord_new~
- When ~elec_coord_new~ is full, send it to the server
#+begin_src f90 :tangle "admc.irp.f"
program admc
call run
call ezfio_finish
end program admc
subroutine run
implicit none
<<declarations>>
! Initialization
if (vmc_algo /= t_Brownian) then
call abrt(irp_here,'DMC should run with Brownian algorithm')
endif
do iter=1,1000
! call read_coords()
k_new = 1
do k_full=1,walk_num
call pdmc_trajectory(k_full, w, E_out(k_full), w_sum(k_full))
<<branching>>
if (k_new >= walk_num) then
w_sum(k_full+1:) = 0.d0
exit
end if
end do
k_new = k_new-1
elec_coord_full(1:elec_num+1,1:3,1:k_new) = &
elec_coord_new(1:elec_num+1,1:3,1:k_new)
! call write_coords(k_new)
call write_energy(walk_num, E_out, w_sum)
end do
end subroutine run
<<read_coords>>
<<write_coords>>
<<write_energy>>
<<pdmc_trajectory>>
#+end_src
** Branching
#+NAME: branching
#+begin_src f90
! Find number of copies
iw = int(w)
w = w - int(w)
if (qmc_ranf() < w) then
iw = iw+1
end if
! Duplicate walker
do l=1,iw
elec_coord_new(1:elec_num+1,1:3,k_new) = &
elec_coord(1:elec_num+1,1:3)
k_new = k_new+1
if (k_new >= walk_num) exit
end do
#+end_src
* Read/write
** Read coordinates
Fetch a new set of coordinates for ~walk_num~ walkers from the pool of coordinates.
#+NAME: read_coords
#+begin_src f90
subroutine read_coords()
implicit none
integer :: i, k
do k=1,walk_num
do i=1,elec_num
read(*,*) elec_coord_full(i,1:3,k)
end do
end do
SOFT_TOUCH elec_coord_full
end subroutine read_coords
#+end_src
** Write coordinates
Send the current set of coordinates for ~walk_num~ walkers to the pool of coordinates.
#+NAME: write_coords
#+begin_src f90
subroutine write_coords()
implicit none
integer :: i, k
do k=1,walk_num
do i=1,elec_num
write(*,*) 'C', elec_coord_full(i,1:3,k)
end do
end do
end subroutine write_coords
#+end_src
** Write energy
Compute the weighted average over the computed energies.
\[
E = \frac{\sum_i w_i E_i}{\sum_i w_i}
\]
#+NAME: write_energy
#+begin_src f90
subroutine write_energy(walk_num_, E_out, w_sum)
implicit none
integer, intent(in) :: walk_num_
double precision, intent(in) :: E_out(walk_num_)
double precision, intent(in) :: w_sum(walk_num_)
integer :: i, k
double precision :: E, S
E = 0.d0
S = 0.d0
do k=1,walk_num
S = S + w_sum(k)
E = E + w_sum(k) * E_out(k)
end do
write(*,*) 'E', E/S, S
end subroutine write_energy
#+end_src
* PDMC trajectory
Computes a PDMC trajectory until the weight ~w~ is $1/2 < w < 3/2$.
The energy of the trajectory is computed as
\[
E = \frac{\sum_i w_i E(R_i)}{\sum_i w_i}
\]
The function returns:
- ~w~: the last of all $w_i$
- ~E_out~: The average energy $E$ of the trajectory
- ~w_sum~: The sum of the weights
#+NAME: declarations_pdmc
#+begin_src f90
integer :: i,j,l
double precision :: delta
! If true, continue to make more steps
logical :: loop
! Max number of steps
integer :: imax
integer, parameter :: nmax=10000
! Brownian step variables
double precision :: p,q
real :: delta_x
logical :: accepted
! Local energies from the past
double precision :: E_loc_save(4)
double precision :: w
#+end_src
#+NAME: pdmc_trajectory
#+begin_src f90
subroutine pdmc_trajectory(k_full, pdmc_weight, E_out, w_sum)
implicit none
integer, intent(in) :: k_full
double precision, intent(out) :: pdmc_weight, E_out, w_sum
<<declarations_pdmc>>
elec_coord(1:elec_num+1,1:3) = elec_coord_full(1:elec_num+1,1:3,k_full)
TOUCH elec_coord
E_out = 0.d0
w_sum = 0.d0
E_loc_save(1:4) = E_loc
pdmc_weight = 1.d0
loop = .True.
do imax = 1, nmax
call brownian_step(p,q,accepted,delta_x)
! delta = (9.d0*E_loc+19.d0*E_loc_save(1)-5.d0*E_loc_save(2)+E_loc_save(3))/24.d0
delta = E_loc
delta = (delta - E_ref)*p
if (delta >= 0.d0) then
w = dexp(-dtime_step*delta)
else
w = 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
if (accepted) then
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
endif
w_sum = w_sum + pdmc_weight
E_out = E_out + pdmc_weight * E_loc
pdmc_weight = pdmc_weight * w
loop = pdmc_weight > 0.5d0 .and. pdmc_weight < 2.0d0
if (.not.loop) exit
end do
E_out = E_out / w_sum
end subroutine pdmc_trajectory
#+end_src

125
src/MAIN/admc.py Executable file
View File

@ -0,0 +1,125 @@
#!/usr/bin/env python3
from mpi4py import MPI
import sys
import gzip
import random
import math
import subprocess
admc_exec = "/home/scemama/qmcchem/src/MAIN/admc"
n_walk_per_proc = 10
def start():
return subprocess.Popen(
[ admc_exec, sys.argv[1] ],
stdin=subprocess.PIPE,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE)
def read(process,len_walk):
line = process.stdout.readline().decode("utf-8").strip()
walk_num = int(line)
walkers = []
print(walk_num)
for k in range(walk_num):
w = []
for i in range(len_walk):
line = process.stdout.readline().decode("utf-8").strip()
w.append( line )
w = '\n'.join(w)
walkers.append(w)
_, E, W = process.stdout.readline().decode("utf-8").split()
return walkers, float(E), float(W)
def write(process, message):
process.stdin.write(f"{message}\n".encode("utf-8"))
process.stdin.flush()
def terminate(process):
process.stdin.close()
process.terminate()
process.wait(timeout=0.2)
def print_energy(EnergyWeight, Energy2Weight, Weight, N):
e = EnergyWeight / Weight
e2 = Energy2Weight / Weight
err = math.sqrt(abs(e*e - e2) / max(1,(N-1)) )
print("%f +/- %f"%(e, err))
return err
def main():
try:
input_dir = sys.argv[1]
except:
print("syntax: argv[0] [FILE]")
sys.exit(-1)
# Pool of electron coordinates
with gzip.open(input_dir+"/electrons/elec_coord_pool.gz","r") as f:
data = f.read().decode("utf-8").split()
len_walk = int(data[1])*int(data[2])
icount = 0
buffer = []
walkers = []
for d in data[4:]:
buffer.append(d)
icount += 1
if (icount == len_walk):
walkers.append(buffer)
buffer = []
icount = 0
walkers = [ '\n'.join(x) for x in walkers ]
do_loop = True
EnergyWeight = 0.
Energy2Weight = 0.
Weight = 0.
NSamples = 0.
# Start processes
proc = start()
while do_loop:
# Once every 1000, shuffle the list of walkers
if random.random() < 0.01:
print("SHUFFLE")
random.shuffle(walkers)
# Pick new walkers
new_coords = walkers[:n_walk_per_proc]
walkers = walkers[n_walk_per_proc:]
# Send new walkers to the process
write(proc, '\n'.join(new_coords))
# Fetch new walkers from the process
new_coords, e_new, w_new = read(proc, len_walk)
walkers += new_coords
# Print energy
ew = e_new * w_new
EnergyWeight += ew
Energy2Weight += e_new * ew
Weight += w_new
NSamples += 1.
print (len(walkers))
err = print_energy(EnergyWeight, Energy2Weight, Weight, NSamples)
if err < 1.e-3:
do_loop = False
terminate(proc)
return
if __name__ == "__main__":
main()

2
src/MAIN/vmc_test.irp.f
View File

@ -1,4 +1,4 @@
program vmc_test
program vmc_test
real :: t1,t2
print *, 'Ndet=',det_num
print *, 'Ndet alpha beta =',det_alpha_num, det_beta_num

1
src/Makefile
View File

@ -1,3 +1,4 @@
IRPF90=irpf90
IRPF90+= $(IRPF90_FLAGS)
include irpf90.make
export