diff --git a/INSTALL.rst b/INSTALL.rst index a9cee40b..979c8edd 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -143,7 +143,7 @@ IRPF90 to Parameters (IRP) method. * Download the latest version of IRPF90 - here : ``_ and move + here : ``_ and move the downloaded archive in the :file:`${QP_ROOT}/external` directory * Extract the archive and go into the :file:`irpf90-*` directory to run diff --git a/REPLACE b/REPLACE index 10a53958..d0b6adce 100755 --- a/REPLACE +++ b/REPLACE @@ -834,3 +834,4 @@ qp_name potential_sr_c_alpha_ao_pbe --rename=potential_c_alpha_ao_sr_pbe qp_name potential_sr_c_beta_ao_pbe --rename=potential_c_beta_ao_sr_pbe qp_name potential_sr_xc_alpha_ao_pbe --rename=potential_xc_alpha_ao_sr_pbe qp_name potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe +qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion diff --git a/bin/qp_tunnel b/bin/qp_tunnel deleted file mode 100755 index 554c28ae..00000000 --- a/bin/qp_tunnel +++ /dev/null @@ -1,71 +0,0 @@ -#!/usr/bin/env python2 - - -""" -Creates an ssh tunnel for using slaves on another network. -Launch a server on the front-end node of the cluster on which the master -process runs. Then start a client ont the front-end node of the distant -cluster. - -Usage: - qp_tunnel server EZFIO_DIR - qp_tunnel client
EZFIO_DIR - -Options: - -h --help - -""" - -import os -import sys -import zmq - -try: - import qp_path -except ImportError: - print "source .quantum_package.rc" - raise - -from docopt import docopt -from ezfio import ezfio - - -def get_address(filename): - with open(os.path.join(filename,'work','qp_run_address'),'r') as f: - a = f.readlines()[0].strip() - return a - - -def set_address(filename,address): - with open(os.path.join(filename,'work','qp_run_address'),'r') as f: - backup = f.readlines() - - with open(os.path.join(filename,'work','qp_run_address'),'w') as f: - f.write('\n'.join([address]+backup)) - - -def main_server(arguments,filename): - destination = get_address(filename) - print destination - - -def main_client(arguments,filename): - destination = arguments["
"] - print destination - - -def main(arguments): - """Main function""" - - print arguments - filename = arguments["EZFIO_DIR"] - - if arguments["server"]: - return main_server(arguments, filename) - if arguments["client"]: - return main_client(arguments, filename) - - -if __name__ == '__main__': - ARGUMENTS = docopt(__doc__) - main(ARGUMENTS) diff --git a/configure b/configure index f8160aef..a940ade2 100755 --- a/configure +++ b/configure @@ -290,6 +290,9 @@ EOF # Special commands for Travis CI chmod +x "${QP_ROOT}"/external/opam_installer.sh rm --force ${QP_ROOT}/bin/opam + if [[ -n ${NO_CACHE} ]] ; then + rm -rf ${HOME}/.opam + fi export OPAMROOT=${HOME}/.opam cat << EOF | bash ${QP_ROOT}/external/opam_installer.sh --no-backup ${QP_ROOT}/bin @@ -310,21 +313,21 @@ EOF else # Conventional commands execute << EOF - chmod +x "\${QP_ROOT}"/external/opam_installer.sh - "\${QP_ROOT}"/external/opam_installer.sh --no-backup + chmod +x "${QP_ROOT}"/external/opam_installer.sh + "${QP_ROOT}"/external/opam_installer.sh --no-backup EOF execute << EOF - rm --force \${QP_ROOT}/bin/opam - export OPAMROOT=\${OPAMROOT:-\${QP_ROOT}/external/opam} - echo \${QP_ROOT}/bin \ - | sh \${QP_ROOT}/external/opam_installer.sh + rm --force ${QP_ROOT}/bin/opam + export OPAMROOT=${OPAMROOT:-${QP_ROOT}/external/opam} + echo ${QP_ROOT}/bin \ + | sh ${QP_ROOT}/external/opam_installer.sh EOF rm ${QP_ROOT}/external/opam_installer.sh # source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true # opam switch create ocaml-base-compiler.4.07.1 || exit 1 + opam init --verbose --yes --compiler=4.07.1 --disable-sandboxing eval $(opam env) -EOF execute << EOF opam install -y \${OCAML_PACKAGES} || exit 1 EOF @@ -423,18 +426,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then fail fi -BWRAP=$(find_exe bwrap) -if [[ ${BWRAP} = $(not_found) ]] ; then - error "Bubblewrap (bwrap) is not installed." - fail -fi - LIBCAP=$(find_lib -lcap) if [[ ${LIBCAP} = $(not_found) ]] ; then error "Libcap (libcap) is not installed." fail fi +BWRAP=$(find_exe bwrap) +if [[ ${BWRAP} = $(not_found) ]] ; then + error "Bubblewrap (bwrap) is not installed." + fail +fi + OPAM=$(find_exe opam) if [[ ${OPAM} = $(not_found) ]] ; then error "OPAM (ocaml) package manager is not installed." diff --git a/data/basis/6-31++g b/data/basis/6-31++g index 7e36d0f9..e1994202 100644 --- a/data/basis/6-31++g +++ b/data/basis/6-31++g @@ -8,6 +8,16 @@ S 1 S 1 1 0.0360000 1.0000000 +HELIUM +S 3 +1 0.3842163400E+02 0.4013973935E-01 +2 0.5778030000E+01 0.2612460970E+00 +3 0.1241774000E+01 0.7931846246E+00 +S 1 +1 0.2979640000E+00 1.0000000 +S 1 +1 0.8600000000E-01 0.1000000000E+01 + LITHIUM S 6 1 642.4189200 0.0021426 diff --git a/data/basis/6-31++g_star_star b/data/basis/6-31++g_star_star index ae124833..4f4ead25 100644 --- a/data/basis/6-31++g_star_star +++ b/data/basis/6-31++g_star_star @@ -1,14 +1,27 @@ + HYDROGEN S 3 - 1 18.7311370 0.03349460 - 2 2.8253937 0.23472695 - 3 0.6401217 0.81375733 +1 0.1873113696E+02 0.3349460434E-01 +2 0.2825394365E+01 0.2347269535E+00 +3 0.6401216923E+00 0.8137573261E+00 S 1 - 1 0.1612778 1.0000000 +1 0.1612777588E+00 1.0000000 S 1 - 1 0.0360000 1.0000000 +1 0.3600000000E-01 0.1000000000E+01 P 1 - 1 1.1000000 1.0000000 +1 0.1100000000E+01 1.0000000 + +HELIUM +S 3 +1 0.3842163400E+02 0.4013973935E-01 +2 0.5778030000E+01 0.2612460970E+00 +3 0.1241774000E+01 0.7931846246E+00 +S 1 +1 0.2979640000E+00 1.0000000 +S 1 +1 0.8600000000E-01 0.1000000000E+01 +P 1 +1 0.1100000000E+01 1.0000000 LITHIUM S 6 diff --git a/data/basis/6-31+g_star b/data/basis/6-31+g_star index 663f7af8..bc112c64 100644 --- a/data/basis/6-31+g_star +++ b/data/basis/6-31+g_star @@ -1,3 +1,19 @@ +HYDROGEN +S 3 +1 0.1873113696E+02 0.3349460434E-01 +2 0.2825394365E+01 0.2347269535E+00 +3 0.6401216923E+00 0.8137573261E+00 +S 1 +1 0.1612777588E+00 1.0000000 + +HELIUM +S 3 +1 0.3842163400E+02 0.4013973935E-01 +2 0.5778030000E+01 0.2612460970E+00 +3 0.1241774000E+01 0.7931846246E+00 +S 1 +1 0.2979640000E+00 1.0000000 + LITHIUM S 6 1 642.4189200 0.0021426 diff --git a/data/basis/6-311++g_2d_2p b/data/basis/6-311++g_2d_2p index 5cf54355..de8f54b4 100644 --- a/data/basis/6-311++g_2d_2p +++ b/data/basis/6-311++g_2d_2p @@ -14,6 +14,18 @@ P 1 P 1 1 0.3750000 1.0000000 +HELIUM +S 3 +1 98.12430 0.0287452 +2 14.76890 0.208061 +3 3.318830 0.837635 +S 1 +1 0.874047 1.000000 +S 1 +1 0.244564 1.000000 +P 1 +1 0.750 1.000000 + LITHIUM S 6 1 900.4600000 0.00228704 diff --git a/data/basis/6-311+g_star b/data/basis/6-311+g_star index 9a6bd86f..14ec3a44 100644 --- a/data/basis/6-311+g_star +++ b/data/basis/6-311+g_star @@ -1,3 +1,23 @@ +HYDROGEN +S 3 +1 33.86500 0.0254938 +2 5.094790 0.190373 +3 1.158790 0.852161 +S 1 +1 0.325840 1.000000 +S 1 +1 0.102741 1.000000 + +HELIUM +S 3 +1 98.12430 0.0287452 +2 14.76890 0.208061 +3 3.318830 0.837635 +S 1 +1 0.874047 1.000000 +S 1 +1 0.244564 1.000000 + LITHIUM S 6 1 900.4600000 0.00228704 diff --git a/data/basis/6-311G_star b/data/basis/6-311G_star index 21621a45..c9cb224c 100644 --- a/data/basis/6-311G_star +++ b/data/basis/6-311G_star @@ -1,3 +1,23 @@ +HYDROGEN +S 3 +1 33.86500 0.0254938 +2 5.094790 0.190373 +3 1.158790 0.852161 +S 1 +1 0.325840 1.000000 +S 1 +1 0.102741 1.000000 + +HELIUM +S 3 +1 98.12430 0.0287452 +2 14.76890 0.208061 +3 3.318830 0.837635 +S 1 +1 0.874047 1.000000 +S 1 +1 0.244564 1.000000 + LITHIUM S 6 1 900.4600000 0.00228704 diff --git a/data/basis/6-31g_star b/data/basis/6-31g_star index 2e79dff2..12291b80 100644 --- a/data/basis/6-31g_star +++ b/data/basis/6-31g_star @@ -1,3 +1,19 @@ +HYDROGEN +S 3 + 1 18.7311370 0.03349460 + 2 2.8253937 0.23472695 + 3 0.6401217 0.81375733 +S 1 + 1 0.1612778 1.0000000 + +HELIUM +S 3 + 1 38.4216340 0.0237660 + 2 5.7780300 0.1546790 + 3 1.2417740 0.4696300 +S 1 + 1 0.2979640 1.0000000 + LITHIUM S 6 1 642.4189200 0.0021426 diff --git a/data/qp2.png b/data/qp2.png index 55dac420..a81e611e 100644 Binary files a/data/qp2.png and b/data/qp2.png differ diff --git a/data/qp2_hd.png b/data/qp2_hd.png new file mode 100644 index 00000000..55dac420 Binary files /dev/null and b/data/qp2_hd.png differ diff --git a/docs/source/research.bib b/docs/source/research.bib index 124e1539..a57bf58a 100644 --- a/docs/source/research.bib +++ b/docs/source/research.bib @@ -1,6 +1,16 @@ %%% ARXIV TO BE UPDATED %%% +@article{Hollett2019Aug, + author = {Hollett, Joshua W. and Loos, Pierre-Fran{\c{c}}ois}, + title = {{Capturing static and dynamic correlation with $\Delta \text{NO}$-MP2 and $\Delta \text{NO}$-CCSD}}, + journal = {arXiv}, + year = {2019}, + month = {Aug}, + eprint = {1908.09914}, + url = {https://arxiv.org/abs/1908.09914} +} + @article{Giner2019Jul, - author = {Giner, Emmanuel and Scemama, Anthony and Toulouse, Julien and Loos, Pierre-Fran{\ifmmode\mbox{\c{c}}\else\c{c}\fi}ois}, + author = {Giner, Emmanuel and Scemama, Anthony and Toulouse, Julien and Loos, Pierre-Fran{\c{c}}ois}, title = {{Chemically Accurate Excitation Energies With Small Basis Sets}}, journal = {arXiv}, year = {2019}, @@ -9,29 +19,66 @@ url = {https://arxiv.org/abs/1907.01245} } -@article{Dash2019May, - author = {Dash, Monika and Feldt, Jonas and Moroni, Saverio and Scemama, Anthony and Filippi, Claudia}, - title = {{Excited states with selected CI-QMC: chemically accurate excitation energies and geometries}}, - journal = {arXiv}, - year = {2019}, - month = {May}, - eprint = {1905.06737}, - url = {https://arxiv.org/abs/1905.06737} -} - -@article{Burton2019May, - author = {Burton, Hugh G. A. and Thom, Alex J. W.}, - title = {{A General Approach for Multireference Ground and Excited States using Non-Orthogonal Configuration Interaction}}, - journal = {arXiv}, - year = {2019}, - month = {May}, - eprint = {1905.02626}, - url = {https://arxiv.org/abs/1905.02626} -} %%%% PUBLISHED PAPERS +@article{Burton2019Sep, + author = {Burton, Hugh G. A. and Thom, Alex J. W.}, + title = {{General Approach for Multireference Ground and Excited States Using Nonorthogonal Configuration Interaction}}, + journal = {J. Chem. Theory Comput.}, + volume = {15}, + number = {9}, + pages = {4851--4861}, + year = {2019}, + month = {Sep}, + issn = {1549-9618}, + publisher = {American Chemical Society}, + doi = {10.1021/acs.jctc.9b00441} +} + +@article{Dash_2019, + author = {Dash, Monika and Feldt, Jonas and Moroni, Saverio and Scemama, Anthony and Filippi, Claudia}, + title = {{Excited States with Selected Configuration Interaction-Quantum Monte Carlo: Chemically Accurate Excitation Energies and Geometries}}, + journal = {J. Chem. Theory Comput.}, + volume = {15}, + number = {9}, + pages = {4896--4906}, + year = {2019}, + month = {Sep}, + issn = {1549-9618}, + publisher = {American Chemical Society}, + doi = {10.1021/acs.jctc.9b00476} +} + +@article{Ferte_2019, + doi = {10.1063/1.5082638}, + url = {https://doi.org/10.1063%2F1.5082638}, + year = 2019, + month = {feb}, + publisher = {{AIP} Publishing}, + volume = {150}, + number = {8}, + pages = {084103}, + author = {Anthony Fert{\'{e}} and Emmanuel Giner and Julien Toulouse}, + title = {Range-separated multideterminant density-functional theory with a short-range correlation functional of the on-top pair density}, + journal = {The Journal of Chemical Physics} +} + +@article{Caffarel_2019, + doi = {10.1063/1.5114703}, + url = {https://doi.org/10.1063%2F1.5114703}, + year = 2019, + month = {aug}, + publisher = {{AIP} Publishing}, + volume = {151}, + number = {6}, + pages = {064101}, + author = {Michel Caffarel}, + title = {Evaluating two-electron-repulsion integrals over arbitrary orbitals using zero variance Monte Carlo: Application to full configuration interaction calculations with Slater-type orbitals}, + journal = {The Journal of Chemical Physics} +} + @article{Loos_2019, doi = {10.1021/acs.jpclett.9b01176}, url = {https://doi.org/10.1021%2Facs.jpclett.9b01176}, diff --git a/docs/source/users_guide/qp_tunnel.rst b/docs/source/users_guide/qp_tunnel.rst new file mode 100644 index 00000000..254216fa --- /dev/null +++ b/docs/source/users_guide/qp_tunnel.rst @@ -0,0 +1,89 @@ +.. _qp_tunnel: + +========= +qp_tunnel +========= + +.. TODO + +.. program:: qp_tunnel + +Establishes a tunnel to allow communications between machines within +different networks, for example multiple MPI slave jobs running on +different clusters. + + +Usage +----- + +.. code:: bash + + qp_tunnel [-g] (ADDRESS|EZFIO_DIR) + +``EZFIO_DIR`` is the name of the |EZFIO| directory containing the data, +and ``ADDRESS`` is the address of another tunnel. + + +.. option:: -h, --help + + Displays the help message + + +.. option:: -g, --get-input + + Download the EZFIO directory from the remote instance of qp_tunnel. + + +Example +------- + +.. code:: text + + +-------------------+ +------------------+ + | | | | + | N1_1 N1_2 N1_3 | | N2_1 N2_2 N2_3 | + | | | | | | | | | | + | +----+----+ | | +----+----+ | + | | | | | | + | C1 F1 | | F2 C2 | + | +---------=----=--------+ | + | | | | + +-------------------+ +------------------+ + + +Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH +on a front-end named respectively F1 and F2. Groups of nodes N1 and N2 have +been reserved by the batch scheduling system on both clusters. +Each node in N1 is on the same network as the other nodes of N1, but they +can't access the network on which the nodes of N2 are. + +1) Start a parallel simulation on the cluster C1, running on nodes N1. + We assume that there is a shared file system, such that F1 can access + the EZFIO directory. We also assume that F1 can communicate with the + nodes of N1. + +2) Run a tunnel on the front-end F1 and keep it running: + +.. code:: bash + + me@f1 $ qp_tunnel my_directory.ezfio + Connect to: + tcp://31.122.230.47:42379 + Ready + +3) On the front-end F2, run another instance connecting to the other one, + which will fetch the |EZFIO| directory: + +.. code:: bash + + me@f2 $ qp_tunnel --get-input tcp://31.122.230.47:42379 + Connect to: + tcp://31.122.209.139:42379 + Communication [ OK ] + Getting input... my_directory.ezfio ...done + Ready + +4) Keep the tunnel running, and you can now run a slave simulation within the + nodes N2. + + diff --git a/docs/source/users_guide/quickstart.rst b/docs/source/users_guide/quickstart.rst index 2883ed96..f0620c5a 100644 --- a/docs/source/users_guide/quickstart.rst +++ b/docs/source/users_guide/quickstart.rst @@ -115,7 +115,7 @@ create an |EZFIO| database with the 6-31G basis set: .. code:: bash - qp create_ezfio -b "6-31G" hcn.xyz -o hcn + qp create_ezfio -b "6-31g" hcn.xyz -o hcn The EZFIO database now contains data relative to the nuclear coordinates and the atomic basis set: diff --git a/man/qp_tunnel.1 b/man/qp_tunnel.1 new file mode 100644 index 00000000..42fa9288 --- /dev/null +++ b/man/qp_tunnel.1 @@ -0,0 +1,139 @@ +.\" Man page generated from reStructuredText. +. +.TH "QP_TUNNEL" "1" "Jun 15, 2019" "2.0" "Quantum Package" +.SH NAME +qp_tunnel \- | Quantum Package > +. +.nr rst2man-indent-level 0 +. +.de1 rstReportMargin +\\$1 \\n[an-margin] +level \\n[rst2man-indent-level] +level margin: \\n[rst2man-indent\\n[rst2man-indent-level]] +- +\\n[rst2man-indent0] +\\n[rst2man-indent1] +\\n[rst2man-indent2] +.. +.de1 INDENT +.\" .rstReportMargin pre: +. RS \\$1 +. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin] +. nr rst2man-indent-level +1 +.\" .rstReportMargin post: +.. +.de UNINDENT +. RE +.\" indent \\n[an-margin] +.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]] +.nr rst2man-indent-level -1 +.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]] +.in \\n[rst2man-indent\\n[rst2man-indent-level]]u +.. +.sp +Establishes a tunnel to allow communications between machines within +different networks, for example multiple MPI slave jobs running on +different clusters. +.SH USAGE +.INDENT 0.0 +.INDENT 3.5 +.sp +.nf +.ft C +qp_tunnel [\-g] (ADDRESS|EZFIO_DIR) +.ft P +.fi +.UNINDENT +.UNINDENT +.sp +\fBEZFIO_DIR\fP is the name of the \fI\%EZFIO\fP directory containing the data, +and \fBADDRESS\fP is the address of another tunnel. +.INDENT 0.0 +.TP +.B \-h, \-\-help +Displays the help message +.UNINDENT +.INDENT 0.0 +.TP +.B \-g, \-\-get\-input +Download the EZFIO directory from the remote instance of qp_tunnel. +.UNINDENT +.SH EXAMPLE +.INDENT 0.0 +.INDENT 3.5 +.sp +.nf +.ft C ++\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +| | | | +| N1_1 N1_2 N1_3 | | N2_1 N2_2 N2_3 | +| | | | | | | | | | +| +\-\-\-\-+\-\-\-\-+ | | +\-\-\-\-+\-\-\-\-+ | +| | | | | | +| C1 F1 | | F2 C2 | +| +\-\-\-\-\-\-\-\-\-=\-\-\-\-=\-\-\-\-\-\-\-\-+ | +| | | | ++\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-+ +.ft P +.fi +.UNINDENT +.UNINDENT +.sp +Imagine you have two clusters, C1 and C2. Each cluster is accessible via SSH +on a front\-end named respectively F1 and F2. Groups of nodes N1 and N2 have +been reserved by the batch scheduling system on both clusters. +Each node in N1 is on the same network as the other nodes of N1, but they +can’t access the network on which the nodes of N2 are. +.INDENT 0.0 +.IP 1. 3 +Start a parallel simulation on the cluster C1, running on nodes N1. +We assume that there is a shared file system, such that F1 can access +the EZFIO directory. We also assume that F1 can communicate with the +nodes of N1. +.IP 2. 3 +Run a tunnel on the front\-end F1 and keep it running: +.UNINDENT +.INDENT 0.0 +.INDENT 3.5 +.sp +.nf +.ft C +me@f1 $ qp_tunnel my_directory.ezfio +Connect to: +tcp://31.122.230.47:42379 +Ready +.ft P +.fi +.UNINDENT +.UNINDENT +.INDENT 0.0 +.IP 3. 3 +On the front\-end F2, run another instance connecting to the other one, +which will fetch the \fI\%EZFIO\fP directory: +.UNINDENT +.INDENT 0.0 +.INDENT 3.5 +.sp +.nf +.ft C +me@f2 $ qp_tunnel \-\-get\-input tcp://31.122.230.47:42379 +Connect to: +tcp://31.122.209.139:42379 +Communication [ OK ] +Getting input... my_directory.ezfio ...done +Ready +.ft P +.fi +.UNINDENT +.UNINDENT +.INDENT 0.0 +.IP 4. 3 +Keep the tunnel running, and you can now run a slave simulation within the +nodes N2. +.UNINDENT +.SH AUTHOR +A. Scemama, E. Giner +.SH COPYRIGHT +2019, A. Scemama, E. Giner +.\" Generated by docutils manpage writer. +. diff --git a/man/rotate_mos.1 b/man/rotate_mos.1 index 54acbc0d..27dcaedb 100644 --- a/man/rotate_mos.1 +++ b/man/rotate_mos.1 @@ -34,7 +34,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]] .INDENT 3.5 Rotates molecular orbitals i and j by combining them as $1/sqrt{2} ( phi_i + phi_j )$ and -$1/sqrt{2} ( phi_i \- phi_j )$. +$1/sqrt{2} ( phi_i - phi_j )$. .sp Needs: .INDENT 0.0 diff --git a/man/test.1 b/man/test.1 new file mode 100644 index 00000000..33d2dc79 --- /dev/null +++ b/man/test.1 @@ -0,0 +1,55 @@ +.\" Man page generated from reStructuredText. +. +.TH "TEST" "1" "Jun 15, 2019" "2.0" "Quantum Package" +.SH NAME +test \- | Quantum Package > +. +.nr rst2man-indent-level 0 +. +.de1 rstReportMargin +\\$1 \\n[an-margin] +level \\n[rst2man-indent-level] +level margin: \\n[rst2man-indent\\n[rst2man-indent-level]] +- +\\n[rst2man-indent0] +\\n[rst2man-indent1] +\\n[rst2man-indent2] +.. +.de1 INDENT +.\" .rstReportMargin pre: +. RS \\$1 +. nr rst2man-indent\\n[rst2man-indent-level] \\n[an-margin] +. nr rst2man-indent-level +1 +.\" .rstReportMargin post: +.. +.de UNINDENT +. RE +.\" indent \\n[an-margin] +.\" old: \\n[rst2man-indent\\n[rst2man-indent-level]] +.nr rst2man-indent-level -1 +.\" new: \\n[rst2man-indent\\n[rst2man-indent-level]] +.in \\n[rst2man-indent\\n[rst2man-indent-level]]u +.. +.INDENT 0.0 +.INDENT 3.5 +Calls: +.INDENT 0.0 +.INDENT 2.0 +.IP \(bu 2 +\fBtwo_e_integrals_index()\fP +.UNINDENT +.INDENT 2.0 +.IP \(bu 2 +\fBtwo_e_integrals_index_reverse()\fP +.UNINDENT +.INDENT 2.0 +.UNINDENT +.UNINDENT +.UNINDENT +.UNINDENT +.SH AUTHOR +A. Scemama, E. Giner +.SH COPYRIGHT +2019, A. Scemama, E. Giner +.\" Generated by docutils manpage writer. +. diff --git a/ocaml/.merlin b/ocaml/.merlin index 3683fed6..84403fc0 100644 --- a/ocaml/.merlin +++ b/ocaml/.merlin @@ -1,4 +1,4 @@ -PKG core ZMQ cryptokit +PKG core zmq cryptokit B _build/ diff --git a/ocaml/Bitlist.ml b/ocaml/Bitlist.ml index 88f9b4dd..b6792bb7 100644 --- a/ocaml/Bitlist.ml +++ b/ocaml/Bitlist.ml @@ -7,82 +7,61 @@ Type for bits strings list of Bits *) -type t = Bit.t list +type t = int64 array + + +let n_int = Array.length + +(* Create a zero bit list *) +let zero n_int = + Array.make (N_int_number.to_int n_int) 0L (* String representation *) let to_string b = - let rec do_work accu = function - | [] -> accu - | head :: tail -> - let new_accu = (Bit.to_string head) ^ accu - in do_work new_accu tail + let int64_to_string x = + String.init 64 (fun i -> + if Int64.logand x @@ Int64.shift_left 1L i <> 0L then + '+' + else + '-') in - do_work "" b + Array.map int64_to_string b + |> Array.to_list + |> String.concat "" let of_string ?(zero='0') ?(one='1') s = - List.init (String.length s) (String.get s) - |> List.rev_map ( fun c -> - if (c = zero) then Bit.Zero - else if (c = one) then Bit.One - else (failwith ("Error in bitstring ") ) ) + let n_int = ( (String.length s - 1) lsr 6 ) + 1 in + let result = Array.make n_int 0L in + String.iteri (fun i c -> + if c = one then + begin + let iint = i lsr 6 in (* i / 64 *) + let k = i - (iint lsl 6) in + result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k; + end) s; + result -let of_string_mp s = - List.init (String.length s) (String.get s) - |> List.rev_map (function - | '-' -> Bit.Zero - | '+' -> Bit.One - | _ -> failwith ("Error in bitstring ") ) +let of_string_mp = of_string ~zero:'-' ~one:'+' (* Create a bit list from an int64 *) -let of_int64 i = - - let rec do_work accu = function - | 0L -> Bit.Zero :: accu |> List.rev - | 1L -> Bit.One :: accu |> List.rev - | i -> - let b = - match (Int64.logand i 1L ) with - | 0L -> Bit.Zero - | 1L -> Bit.One - | _ -> raise (Failure "i land 1 not in (0,1)") - in - do_work (b :: accu) (Int64.shift_right_logical i 1) - in - - let adjust_length result = - let rec do_work accu = function - | 64 -> List.rev accu - | i when i>64 -> raise (Failure "Error in of_int64 > 64") - | i when i<0 -> raise (Failure "Error in of_int64 < 0") - | i -> do_work (Bit.Zero :: accu) (i+1) - in - do_work (List.rev result) (List.length result) - in - adjust_length (do_work [] i) - +let of_int64 i = [| i |] (* Create an int64 from a bit list *) -let to_int64 l = - assert ( (List.length l) <= 64) ; - let rec do_work accu = function - | [] -> accu - | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail - | Bit.One::tail -> do_work Int64.(logor one (shift_left accu 1)) tail - in do_work Int64.zero (List.rev l) +let to_int64 = function + | [| i |] -> i + | _ -> failwith "N_int > 1" -(* Create a bit list from a list of int64 *) -let of_int64_list l = - List.map of_int64 l - |> List.concat - (* Create a bit list from an array of int64 *) -let of_int64_array l = - Array.map of_int64 l - |> Array.to_list - |> List.concat +external of_int64_array : int64 array -> t = "%identity" +external to_int64_array : t -> int64 array = "%identity" + + +(* Create a bit list from a list of int64 *) +let of_int64_list l = + Array.of_list l |> of_int64_array (* Compute n_int *) @@ -91,101 +70,64 @@ let n_int_of_mo_num mo_num = N_int_number.of_int ( (mo_num-1)/bit_kind_size + 1 ) -(* Create a zero bit list *) -let zero n_int = - let n_int = N_int_number.to_int n_int in - let a = Array.init n_int (fun i-> 0L) in - of_int64_list ( Array.to_list a ) - (* Create an int64 list from a bit list *) let to_int64_list l = - let rec do_work accu buf counter = function - | [] -> - begin - match buf with - | [] -> accu - | _ -> (List.rev buf)::accu - end - | i::tail -> - if (counter < 64) then - do_work accu (i::buf) (counter+1) tail - else - do_work ( (List.rev (i::buf))::accu) [] 1 tail - in - let l = do_work [] [] 1 l - in - List.rev_map to_int64 l + to_int64_array l |> Array.to_list -(* Create an array of int64 from a bit list *) -let to_int64_array l = - to_int64_list l - |> Array.of_list (* Create a bit list from a list of MO indices *) let of_mo_number_list n_int l = - let n_int = N_int_number.to_int n_int in - let length = n_int*64 in - let a = Array.make length (Bit.Zero) in - List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l; - Array.to_list a + let result = zero n_int in + List.iter (fun j -> + let i = (MO_number.to_int j) - 1 in + let iint = i lsr 6 in (* i / 64 *) + let k = i - (iint lsl 6) in + result.(iint) <- Int64.logor result.(iint) @@ Int64.shift_left 1L k; + ) l; + result let to_mo_number_list l = - let a = Array.of_list l in - let mo_num = MO_number.get_max () in - let rec do_work accu = function - | 0 -> accu - | i -> - begin - let new_accu = - match a.(i-1) with - | Bit.One -> (MO_number.of_int ~max:mo_num i)::accu - | Bit.Zero -> accu - in - do_work new_accu (i-1) - end + let rec aux_one x shift accu = function + | -1 -> accu + | i -> if Int64.logand x (Int64.shift_left 1L i) <> 0L then + aux_one x shift ( (i+shift) ::accu) (i-1) + else + aux_one x shift accu (i-1) in - do_work [] (List.length l) - + Array.mapi (fun i x -> + let shift = (i lsr 6) lsl 6 + 1 in + aux_one x shift [] 63 + ) l + |> Array.to_list + |> List.concat + |> List.map MO_number.of_int (* logical operations on bit_list *) -let logical_operator2 op a b = - let rec do_work_binary result a b = - match a, b with - | [], [] -> result - | [], _ | _ , [] -> raise (Failure "Lists should have same length") - | (ha::ta), (hb::tb) -> - let newbit = op ha hb - in do_work_binary (newbit::result) ta tb +let and_operator a b = Array.map2 Int64.logand a b +let xor_operator a b = Array.map2 Int64.logxor a b +let or_operator a b = Array.map2 Int64.logor a b +let not_operator b = Array.map Int64.lognot b + + + +let pop_sign = + let mask = + (Int64.pred (Int64.shift_left 1L 63)) in - List.rev (do_work_binary [] a b) - - -let logical_operator1 op b = - let rec do_work_unary result b = - match b with - | [] -> result - | (hb::tb) -> - let newbit = op hb - in do_work_unary (newbit::result) tb - in - List.rev (do_work_unary [] b) - - -let and_operator a b = logical_operator2 Bit.and_operator a b -let xor_operator a b = logical_operator2 Bit.xor_operator a b -let or_operator a b = logical_operator2 Bit.or_operator a b -let not_operator b = logical_operator1 Bit.not_operator b + fun x -> Int64.logand mask x let popcnt b = - List.fold_left (fun accu -> function - | Bit.One -> accu+1 - | Bit.Zero -> accu - ) 0 b + Array.fold_left (fun accu x -> + if x >= 0L then + accu + (Z.popcount @@ Z.of_int64 x) + else + accu + 1 + (Z.popcount @@ Z.of_int64 (pop_sign x)) + ) 0 b diff --git a/ocaml/Bitlist.mli b/ocaml/Bitlist.mli index 1f230f08..8d421c7f 100644 --- a/ocaml/Bitlist.mli +++ b/ocaml/Bitlist.mli @@ -1,4 +1,4 @@ -type t = Bit.t list +type t (** The zero bit list *) val zero : Qptypes.N_int_number.t -> t diff --git a/ocaml/Determinant.ml b/ocaml/Determinant.ml index c27aefc4..0a18a4c6 100644 --- a/ocaml/Determinant.ml +++ b/ocaml/Determinant.ml @@ -3,7 +3,8 @@ open Sexplib.Std type t = int64 array [@@deriving sexp] -let to_int64_array (x:t) = (x:int64 array) +external to_int64_array : t -> int64 array = "%identity" +external of_int64_array_no_check : int64 array -> t = "%identity" let to_alpha_beta x = @@ -24,19 +25,6 @@ let to_bitlist_couple x = in (xa,xb) -let bitlist_to_string ~mo_num x = - let len = - MO_number.to_int mo_num - in - let s = - List.map (function - | Bit.Zero -> "-" - | Bit.One -> "+" - ) x - |> String.concat "" - in - String.sub s 0 len - let of_int64_array ~n_int ~alpha ~beta x = @@ -47,38 +35,29 @@ let of_int64_array ~n_int ~alpha ~beta x = in if ( (Bitlist.popcnt a) <> alpha) then begin - let mo_num = MO_number.get_max () in - let mo_num = MO_number.of_int mo_num ~max:mo_num in failwith (Printf.sprintf "Expected %d electrons in alpha determinant -%s" alpha (bitlist_to_string ~mo_num:mo_num a) ) +%s" alpha (Bitlist.to_string a) ) end; if ( (Bitlist.popcnt b) <> beta ) then begin - let mo_num = MO_number.get_max () in - let mo_num = MO_number.of_int mo_num ~max:mo_num in failwith (Printf.sprintf "Expected %d electrons in beta determinant -%s" beta (bitlist_to_string ~mo_num:mo_num b) ) +%s" beta (Bitlist.to_string b) ) end; x -let of_int64_array_no_check x = x -let of_bitlist_couple ?n_int ~alpha ~beta (xa,xb) = +let of_bitlist_couple ~n_int ~alpha ~beta (xa,xb) = let ba, bb = Bitlist.to_int64_array xa , Bitlist.to_int64_array xb - and n_int = - match n_int with - | Some x -> x - | None -> Bitlist.n_int_of_mo_num (List.length xa) in of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb]) let to_string ~mo_num x = let (xa,xb) = to_bitlist_couple x in - [ " " ; bitlist_to_string ~mo_num xa ; "\n" ; - " " ; bitlist_to_string ~mo_num xb ] + [ " " ; Bitlist.to_string xa ; "\n" ; + " " ; Bitlist.to_string xb ] |> String.concat "" diff --git a/ocaml/Determinant.mli b/ocaml/Determinant.mli index 49ba1057..73f91dc7 100644 --- a/ocaml/Determinant.mli +++ b/ocaml/Determinant.mli @@ -24,7 +24,7 @@ val to_alpha_beta : t -> (int64 array)*(int64 array) val to_bitlist_couple : t -> Bitlist.t * Bitlist.t (** Create from a bit list *) -val of_bitlist_couple : ?n_int:Qptypes.N_int_number.t -> +val of_bitlist_couple : n_int:Qptypes.N_int_number.t -> alpha:Qptypes.Elec_alpha_number.t -> beta:Qptypes.Elec_beta_number.t -> Bitlist.t * Bitlist.t -> t diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index e4c6ff2a..6c449c1b 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -7,14 +7,14 @@ module Determinants_by_hand : sig { n_int : N_int_number.t; bit_kind : Bit_kind.t; n_det : Det_number.t; + n_det_qp_edit : Det_number.t; n_states : States_number.t; expected_s2 : Positive_float.t; psi_coef : Det_coef.t array; psi_det : Determinant.t array; state_average_weight : Positive_float.t array; } [@@deriving sexp] - val read : unit -> t - val read_maybe : unit -> t option + val read : ?full:bool -> unit -> t option val write : t -> unit val to_string : t -> string val to_rst : t -> Rst_string.t @@ -28,6 +28,7 @@ end = struct { n_int : N_int_number.t; bit_kind : Bit_kind.t; n_det : Det_number.t; + n_det_qp_edit : Det_number.t; n_states : States_number.t; expected_s2 : Positive_float.t; psi_coef : Det_coef.t array; @@ -38,8 +39,6 @@ end = struct let get_default = Qpackage.get_ezfio_default "determinants";; - let n_det_read_max = 10_000 ;; - let read_n_int () = if not (Ezfio.has_determinants_n_int()) then Ezfio.get_mo_basis_mo_num () @@ -80,11 +79,27 @@ end = struct |> Det_number.of_int ;; + let read_n_det_qp_edit () = + if not (Ezfio.has_determinants_n_det_qp_edit ()) then + begin + let n_det = read_n_det () |> Det_number.to_int in + Ezfio.set_determinants_n_det_qp_edit n_det + end; + Ezfio.get_determinants_n_det_qp_edit () + |> Det_number.of_int + ;; + let write_n_det n = Det_number.to_int n |> Ezfio.set_determinants_n_det ;; + let write_n_det_qp_edit n = + let n_det = read_n_det () |> Det_number.to_int in + min n_det (Det_number.to_int n) + |> Ezfio.set_determinants_n_det_qp_edit + ;; + let read_n_states () = if not (Ezfio.has_determinants_n_states ()) then Ezfio.set_determinants_n_states 1 @@ -178,7 +193,7 @@ end = struct |> Ezfio.set_determinants_expected_s2 ;; - let read_psi_coef () = + let read_psi_coef ~read_only () = if not (Ezfio.has_determinants_psi_coef ()) then begin let n_states = @@ -189,24 +204,33 @@ end = struct ~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. )) |> Ezfio.set_determinants_psi_coef end; - Ezfio.get_determinants_psi_coef () + begin + if read_only then + Ezfio.get_determinants_psi_coef_qp_edit () + else + Ezfio.get_determinants_psi_coef () + end |> Ezfio.flattened_ezfio |> Array.map Det_coef.of_float ;; let write_psi_coef ~n_det ~n_states c = let n_det = Det_number.to_int n_det - and c = Array.to_list c - |> List.map Det_coef.to_float + and c = + Array.map Det_coef.to_float c + |> Array.to_list and n_states = States_number.to_int n_states in - Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c - |> Ezfio.set_determinants_psi_coef + let r = + Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c + in + Ezfio.set_determinants_psi_coef r; + Ezfio.set_determinants_psi_coef_qp_edit r ;; - let read_psi_det () = + let read_psi_det ~read_only () = let n_int = read_n_int () and n_alpha = Ezfio.get_electrons_elec_alpha_num () |> Elec_alpha_number.of_int @@ -232,19 +256,26 @@ end = struct |> Ezfio.set_determinants_psi_det ; end ; let n_int = N_int_number.to_int n_int in - let psi_det_array = Ezfio.get_determinants_psi_det () in + let psi_det_array = + if read_only then + Ezfio.get_determinants_psi_det_qp_edit () + else + Ezfio.get_determinants_psi_det () + in let dim = psi_det_array.Ezfio.dim and data = Ezfio.flattened_ezfio psi_det_array in assert (n_int = dim.(0)); assert (dim.(1) = 2); - assert (dim.(2) = (Det_number.to_int (read_n_det ()))); - List.init dim.(2) (fun i -> + if read_only then + assert (dim.(2) = (Det_number.to_int (read_n_det_qp_edit ()))) + else + assert (dim.(2) = (Det_number.to_int (read_n_det ()))); + Array.init dim.(2) (fun i -> Array.sub data (2*n_int*i) (2*n_int) ) - |> List.map (Determinant.of_int64_array + |> Array.map (Determinant.of_int64_array ~n_int:(N_int_number.of_int n_int) ~alpha:n_alpha ~beta:n_beta ) - |> Array.of_list ;; let write_psi_det ~n_int ~n_det d = @@ -252,40 +283,45 @@ end = struct |> Array.concat |> Array.to_list in - Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data - |> Ezfio.set_determinants_psi_det + let r = + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data + in + Ezfio.set_determinants_psi_det r; + Ezfio.set_determinants_psi_det_qp_edit r ;; - let read () = + let read ?(full=true) () = + + let n_det_qp_edit = read_n_det_qp_edit () in + let n_det = read_n_det () in + let read_only = + if full then false else n_det_qp_edit <> n_det + in + if (Ezfio.has_mo_basis_mo_num ()) then + try + Some { n_int = read_n_int () ; bit_kind = read_bit_kind () ; n_det = read_n_det () ; + n_det_qp_edit = read_n_det_qp_edit () ; expected_s2 = read_expected_s2 () ; - psi_coef = read_psi_coef () ; - psi_det = read_psi_det () ; + psi_coef = read_psi_coef ~read_only () ; + psi_det = read_psi_det ~read_only () ; n_states = read_n_states () ; state_average_weight = read_state_average_weight () ; } + with _ -> None else - failwith "No molecular orbitals, so no determinants" - ;; - - let read_maybe () = - let n_det = - read_n_det () - in - if ( (Det_number.to_int n_det) < n_det_read_max ) then - try Some (read ()) with - | Failure _ -> None - else + (* No molecular orbitals, so no determinants *) None ;; let write { n_int ; bit_kind ; n_det ; + n_det_qp_edit ; expected_s2 ; psi_coef ; psi_det ; @@ -297,9 +333,13 @@ end = struct write_n_det n_det; write_n_states n_states; write_expected_s2 expected_s2; - write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ; - write_psi_det ~n_int:n_int ~n_det:n_det psi_det; - write_state_average_weight state_average_weight; + if n_det <= n_det_qp_edit then + begin + write_n_det_qp_edit n_det; + write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ; + write_psi_det ~n_int:n_int ~n_det:n_det psi_det + end; + write_state_average_weight state_average_weight ;; @@ -316,11 +356,13 @@ end = struct |> States_number.to_int and ndet = Det_number.to_int b.n_det + and ndet_qp_edit = + Det_number.to_int b.n_det_qp_edit in let coefs_string i = Array.init nstates (fun j -> let ishift = - j*ndet + j*ndet_qp_edit in if (ishift < Array.length b.psi_coef) then b.psi_coef.(i+ishift) @@ -331,7 +373,7 @@ end = struct ) |> Array.to_list |> String.concat "\t" in - Array.init ndet (fun i -> + Array.init ndet_qp_edit (fun i -> Printf.sprintf " %s\n%s\n" (coefs_string i) (Determinant.to_string ~mo_num:mo_num b.psi_det.(i) @@ -363,7 +405,7 @@ Determinants :: " (b.expected_s2 |> Positive_float.to_string) (b.n_det |> Det_number.to_string) - (b.state_average_weight |> Array.to_list |> List.map Positive_float.to_string |> String.concat "\t") + (b.state_average_weight |> Array.map Positive_float.to_string |> Array.to_list |> String.concat "\t") det_text |> Rst_string.of_string ;; @@ -387,10 +429,10 @@ psi_det = %s (b.n_states |> States_number.to_string) (b.expected_s2 |> Positive_float.to_string) (b.state_average_weight |> Array.to_list |> List.map Positive_float.to_string |> String.concat ",") - (b.psi_coef |> Array.to_list |> List.map Det_coef.to_string + (b.psi_coef |> Array.map Det_coef.to_string |> Array.to_list |> String.concat ", ") - (b.psi_det |> Array.to_list |> List.map (Determinant.to_string - ~mo_num) |> String.concat "\n\n") + (b.psi_det |> Array.map (Determinant.to_string ~mo_num) |> Array.to_list + |> String.concat "\n\n") ;; let of_rst r = @@ -472,6 +514,7 @@ psi_det = %s (* Handle determinants *) let psi_det = + let n_int = N_int_number.of_int @@ (MO_number.get_max () - 1) / 64 + 1 in let n_alpha = Ezfio.get_electrons_elec_alpha_num () |> Elec_alpha_number.of_int and n_beta = Ezfio.get_electrons_elec_beta_num () @@ -483,8 +526,8 @@ psi_det = %s begin let newdet = (Bitlist.of_string ~zero:'-' ~one:'+' alpha , - Bitlist.of_string ~zero:'-' ~one:'+' beta) - |> Determinant.of_bitlist_couple ~alpha:n_alpha ~beta:n_beta + Bitlist.of_string ~zero:'-' ~one:'+' beta) + |> Determinant.of_bitlist_couple ~n_int ~alpha:n_alpha ~beta:n_beta |> Determinant.sexp_of_t |> Sexplib.Sexp.to_string in @@ -492,9 +535,6 @@ psi_det = %s end | _::tail -> read_dets accu tail in - let dets = - List.map String_ext.rev dets - in let a = read_dets [] dets |> String.concat "" @@ -510,9 +550,11 @@ psi_det = %s Printf.sprintf "(n_int %d)" (N_int_number.get_max ()) and n_states = Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_n_states ()) + and n_det_qp_edit = + Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ()) in let s = - String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det] + String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ] in @@ -527,7 +569,9 @@ psi_det = %s Det_number.to_int n_det_new in let det = - read () + match read () with + | Some x -> x + | None -> failwith "No determinants in file" in let n_det_old, n_states = Det_number.to_int det.n_det, @@ -558,7 +602,9 @@ psi_det = %s let extract_state istate = Printf.printf "Extracting state %d\n" (States_number.to_int istate); let det = - read () + match read () with + | Some x -> x + | None -> failwith "No determinants in file" in let n_det, n_states = Det_number.to_int det.n_det, @@ -588,7 +634,9 @@ psi_det = %s let extract_states range = Printf.printf "Extracting states %s\n" (Range.to_string range); let det = - read () + match read () with + | Some x -> x + | None -> failwith "No determinants in file" in let n_det, n_states = Det_number.to_int det.n_det, @@ -614,7 +662,8 @@ psi_det = %s j*n_det in for i=0 to (n_det-1) do - det.psi_coef.(!state_shift+i) <- det.psi_coef.(i+ishift) + det.psi_coef.(!state_shift+i) <- + det.psi_coef.(i+ishift) done end; state_shift := !state_shift + n_det diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index 1402845f..12654aad 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -65,8 +65,15 @@ end = struct let read_mo_num () = - Ezfio.get_mo_basis_mo_num () - |> MO_number.of_int + let elec_alpha_num = + Ezfio.get_electrons_elec_alpha_num () + in + let result = + Ezfio.get_mo_basis_mo_num () + in + if result < elec_alpha_num then + failwith "More alpha electrons than MOs"; + MO_number.of_int result let read_mo_class () = diff --git a/ocaml/Makefile b/ocaml/Makefile index aaf8c2cc..6ff91273 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m $(QP_ROOT)/ocaml/element_create_db.byte external_libs: - opam install cryptokit core + opam install cryptokit sexplib qpackage.odocl: $(MLIFILES) ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl diff --git a/ocaml/_tags b/ocaml/_tags index 8c354c8b..55b1c681 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,4 +1,4 @@ -true: package(cryptokit,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt) +true: package(cryptokit,zarith,zmq,str,sexplib,ppx_sexp_conv,ppx_deriving,getopt) true: thread false: profile <*byte> : linkdep(c_bindings.o), custom diff --git a/ocaml/qp_print_basis.ml b/ocaml/qp_print_basis.ml index 74c36761..d34894ad 100644 --- a/ocaml/qp_print_basis.ml +++ b/ocaml/qp_print_basis.ml @@ -44,9 +44,12 @@ let psi_det () = let psi_det = Input.Determinants_by_hand.read () in - Input.Determinants_by_hand.to_rst psi_det - |> Rst_string.to_string - |> print_endline + match psi_det with + | Some psi_det -> + Input.Determinants_by_hand.to_rst psi_det + |> Rst_string.to_string + |> print_endline + | None -> () diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index e806082c..942e2cc2 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/ocaml/qp_set_mo_class.ml @@ -112,9 +112,9 @@ let set ~core ~inact ~act ~virt ~del = and av = Excitation.create_single act virt in let single_excitations = [ ia ; aa ; av ] - |> List.map (fun x -> + |> List.map (fun z -> let open Excitation in - match x with + match z with | Single (x,y) -> ( MO_class.to_bitlist n_int (Hole.to_mo_class x), MO_class.to_bitlist n_int (Particle.to_mo_class y) ) @@ -187,9 +187,10 @@ let set ~core ~inact ~act ~virt ~del = match aa with | Double _ -> assert false | Single (x,y) -> - ( MO_class.to_bitlist n_int (Hole.to_mo_class x) ) @ - ( MO_class.to_bitlist n_int (Particle.to_mo_class y) ) - |> Bitlist.to_int64_list + Bitlist.to_int64_list + ( MO_class.to_bitlist n_int ( Hole.to_mo_class x) ) @ + Bitlist.to_int64_list + ( MO_class.to_bitlist n_int (Particle.to_mo_class y) ) in Ezfio.set_bitmasks_n_mask_cas 1; Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result diff --git a/ocaml/qp_tunnel.ml b/ocaml/qp_tunnel.ml new file mode 100644 index 00000000..75112e66 --- /dev/null +++ b/ocaml/qp_tunnel.ml @@ -0,0 +1,470 @@ +open Qputils +open Qptypes + +type ezfio_or_address = EZFIO of string | ADDRESS of string +type req_or_sub = REQ | SUB + +let localport = 42379 + + +let in_time_sum = ref 1.e-9 +and in_size_sum = ref 0. + + +let () = + let open Command_line in + begin + "Creates an ssh tunnel for using slaves on another network. Launch a server on the front-end node of the cluster on which the master process runs. Then start a client ont the front-end node of the distant cluster." + + |> set_footer_doc ; + + [ { short='g' ; long="get-input" ; opt=Optional ; + doc="Downloads the EZFIO directory." ; + arg=Without_arg; } ; + + { short='v' ; long="verbose" ; opt=Optional ; + doc="Prints the transfer speed." ; + arg=Without_arg; } ; + + anonymous + "(EZFIO_DIR|ADDRESS)" + Mandatory + "EZFIO directory or address."; + ] |> set_specs + end; + + let arg = + let x = + match Command_line.anon_args () with + | [x] -> x + | _ -> begin + Command_line.help () ; + failwith "EZFIO_FILE or ADDRESS is missing" + end + in + if Sys.file_exists x && Sys.is_directory x then + EZFIO x + else + ADDRESS x + in + + let verbose = + Command_line.get_bool "verbose" + in + + + + let localhost = + Lazy.force TaskServer.ip_address + in + + + let long_address = + match arg with + | ADDRESS x -> x + | EZFIO x -> + let ic = + Filename.concat (Qpackage.ezfio_work x) "qp_run_address" + |> open_in + in + let result = + input_line ic + |> String.trim + in + close_in ic; + result + in + + let protocol, address, port = + match String.split_on_char ':' long_address with + | t :: a :: p :: [] -> t, a, int_of_string p + | _ -> failwith @@ + Printf.sprintf "%s : Malformed address" long_address + in + + + let zmq_context = + Zmq.Context.create () + in + + + (** Check availability of the ports *) + let localport = + let dummy_socket = + Zmq.Socket.create zmq_context Zmq.Socket.rep + in + let rec try_new_port port_number = + try + List.iter (fun i -> + let address = + Printf.sprintf "tcp://%s:%d" localhost (port_number+i) + in + Zmq.Socket.bind dummy_socket address; + Zmq.Socket.unbind dummy_socket address + ) [ 0;1;2;3;4;5;6;7;8;9 ] ; + port_number + with + | Unix.Unix_error _ -> try_new_port (port_number+100) + in + let result = + try_new_port localport + in + Zmq.Socket.close dummy_socket; + result + in + + + let create_socket sock_type bind_or_connect addr = + let socket = + Zmq.Socket.create zmq_context sock_type + in + let () = + try + bind_or_connect socket addr + with + | _ -> failwith @@ + Printf.sprintf "Unable to establish connection to %s." addr + in + socket + in + + + (* Handle termination *) + let run_status = ref true in + let handler = + Sys.Signal_handle (fun signum -> + run_status := false; + Sys.set_signal signum Sys.Signal_default + ) + in + Sys.set_signal Sys.sigusr1 handler; + Sys.set_signal Sys.sigint handler; + + + let new_thread req_or_sub addr_in addr_out = + let socket_in, socket_out = + match req_or_sub with + | REQ -> + create_socket Zmq.Socket.router Zmq.Socket.bind addr_in, + create_socket Zmq.Socket.dealer Zmq.Socket.connect addr_out + | SUB -> + create_socket Zmq.Socket.sub Zmq.Socket.connect addr_in, + create_socket Zmq.Socket.pub Zmq.Socket.bind addr_out + in + + if req_or_sub = SUB then + Zmq.Socket.subscribe socket_in ""; + + + (* + let action = + if verbose then + begin + match req_or_sub with + | REQ -> (fun () -> + let msg = + Zmq.Socket.recv_all socket_in + in + let t0 = Unix.gettimeofday () in + Zmq.Socket.send_all socket_out msg; + let in_size = + float_of_int ( List.fold_left (fun accu x -> accu + String.length x) 0 msg ) + /. 8192. /. 1024. + in + let msg = + Zmq.Socket.recv_all socket_out + in + let t1 = Unix.gettimeofday () in + Zmq.Socket.send_all socket_in msg; + let in_time = t1 -. t0 in + in_time_sum := !in_time_sum +. in_time; + in_size_sum := !in_size_sum +. in_size; + Printf.printf " %16.2f MiB/s -- %16.2f MiB/s\n%!" (in_size /. in_time) (!in_size_sum /. !in_time_sum); + ) + | SUB -> (fun () -> + Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + end + else + begin + match req_or_sub with + | REQ -> (fun () -> + Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out; + Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in ) + | SUB -> (fun () -> + Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + end + in + *) + + let action_in = + match req_or_sub with + | REQ -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + | SUB -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + in + + let action_out = + match req_or_sub with + | REQ -> (fun () -> Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in ) + | SUB -> (fun () -> () ) + in + + let pollitem = + Zmq.Poll.mask_of + [| (socket_in, Zmq.Poll.In) ; (socket_out, Zmq.Poll.In) |] + in + + + while !run_status do + + let polling = + Zmq.Poll.poll ~timeout:1000 pollitem + in + + match polling with + | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |] -> ( action_out () ; action_in () ) + | [| _ ; Some Zmq.Poll.In |] -> action_out () + | [| Some Zmq.Poll.In ; _ |] -> action_in () + | _ -> () + done; + + Zmq.Socket.close socket_in; + Zmq.Socket.close socket_out; + in + + + + let ocaml_thread = + let addr_out = + Printf.sprintf "tcp:%s:%d" address port + in + + let addr_in = + Printf.sprintf "tcp://*:%d" localport + in + + let f () = + new_thread REQ addr_in addr_out + in + + (Thread.create f) () + in + Printf.printf "Connect to:\ntcp://%s:%d\n%!" localhost localport; + + + let fortran_thread = + let addr_out = + Printf.sprintf "tcp:%s:%d" address (port+2) + in + + let addr_in = + Printf.sprintf "tcp://*:%d" (localport+2) + in + + let f () = + new_thread REQ addr_in addr_out + in + (Thread.create f) () + in + + + let pub_thread = + let addr_in = + Printf.sprintf "tcp:%s:%d" address (port+1) + in + + let addr_out = + Printf.sprintf "tcp://*:%d" (localport+1) + in + + let f () = + new_thread SUB addr_in addr_out + in + (Thread.create f) () + in + + + + let input_thread = + let f () = + let addr_out = + match arg with + | EZFIO _ -> None + | ADDRESS _ -> Some ( + Printf.sprintf "tcp:%s:%d" address (port+9) ) + in + + let addr_in = + Printf.sprintf "tcp://*:%d" (localport+9) + in + + let socket_in = + create_socket Zmq.Socket.rep Zmq.Socket.bind addr_in + in + + let socket_out = + match addr_out with + | Some addr_out -> Some ( + create_socket Zmq.Socket.req Zmq.Socket.connect addr_out) + | None -> None + in + + let temp_file = + Filename.temp_file "qp_tunnel" ".tar.gz" + in + + let get_ezfio_filename () = + match arg with + | EZFIO x -> x + | ADDRESS _ -> + begin + match socket_out with + | None -> assert false + | Some socket_out -> ( + Zmq.Socket.send socket_out "get_ezfio_filename" ; + Zmq.Socket.recv socket_out +) + end + in + + let get_input () = + match arg with + | EZFIO x -> + begin + Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x + |> Sys.command |> ignore; + let fd = + Unix.openfile temp_file [Unix.O_RDONLY] 0o640 + in + let len = + Unix.lseek fd 0 Unix.SEEK_END + in + ignore @@ Unix.lseek fd 0 Unix.SEEK_SET ; + let bstr = + Unix.map_file fd Bigarray.char + Bigarray.c_layout false [| len |] + |> Bigarray.array1_of_genarray + in + let result = + String.init len (fun i -> bstr.{i}) ; + in + Unix.close fd; + Sys.remove temp_file; + result + end + | ADDRESS _ -> + begin + match socket_out with + | None -> assert false + | Some socket_out -> ( + Zmq.Socket.send socket_out "get_input" ; + Zmq.Socket.recv socket_out +) + end + in + + let () = + match socket_out with + | None -> () + | Some socket_out -> + Zmq.Socket.send socket_out "test"; + Printf.printf "Communication [ %s ]\n%!" (Zmq.Socket.recv socket_out); + in + + (* Download input if asked *) + if Command_line.get_bool "get-input" then + begin + match arg with + | EZFIO _ -> () + | ADDRESS _ -> + begin + Printf.printf "Getting input... %!"; + let ezfio_filename = + get_ezfio_filename () + in + Printf.printf "%s%!" ezfio_filename; + let oc = + open_out temp_file + in + get_input () + |> output_string oc; + close_out oc; + Printf.sprintf "tar -zxf %s" temp_file + |> Sys.command |> ignore ; + let oc = + Filename.concat (Qpackage.ezfio_work ezfio_filename) "qp_run_address" + |> open_out + in + Printf.fprintf oc "tcp://%s:%d\n" localhost localport; + close_out oc; + Printf.printf " ...done\n%!" + end + end; + + (* Main loop *) + let pollitem = + Zmq.Poll.mask_of [| (socket_in, Zmq.Poll.In) |] + in + + let action () = + match Zmq.Socket.recv socket_in with + | "get_input" -> get_input () + |> Zmq.Socket.send socket_in + | "get_ezfio_filename" -> get_ezfio_filename () + |> Zmq.Socket.send socket_in + | "test" -> Zmq.Socket.send socket_in "OK" + | x -> Printf.sprintf "Message '%s' not understood" x + |> Zmq.Socket.send socket_in + in + + Printf.printf " +On remote hosts, create ssh tunnel using: + ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s & +Or from this host connect to clients using: + ssh -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d & +%!" + (port ) localhost (localport ) + (port+1) localhost (localport+1) + (port+2) localhost (localport+2) + (port+9) localhost (localport+9) + (Unix.gethostname ()) + (port ) (localport ) + (port+1) (localport+1) + (port+2) (localport+2) + (port+9) (localport+9); + Printf.printf "Ready\n%!"; + while !run_status do + + let polling = + Zmq.Poll.poll ~timeout:1000 pollitem + in + + match polling.(0) with + | Some Zmq.Poll.In -> action () + | None -> () + | Some Zmq.Poll.In_out + | Some Zmq.Poll.Out -> () + + done; + + let () = + match socket_out with + | Some socket_out -> Zmq.Socket.close socket_out + | None -> () + in + Zmq.Socket.close socket_in + in + + (Thread.create f) () + in + + (* Termination *) + Thread.join input_thread; + Thread.join fortran_thread; + Thread.join pub_thread; + Thread.join ocaml_thread; + Zmq.Context.terminate zmq_context; + Printf.printf "qp_tunnel exited properly.\n" + + + + diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index 353c91c4..a63a19cc 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -58,7 +58,7 @@ let input_data = " * Det_number_max : int assert (x > 0) ; - if (x > 50_00_000_000) then + if (x > 50_000_000_000) then warning \"More than 50 billion determinants\"; * States_number : int diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index b5196294..d7c3fd32 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -79,7 +79,7 @@ let get s = | Ao_basis -> f Ao_basis.(read, to_rst) | Determinants_by_hand -> - f Determinants_by_hand.(read_maybe, to_rst) + f Determinants_by_hand.(read ~full:false, to_rst) {section_to_rst} end with diff --git a/src/ao_one_e_ints/README.rst b/src/ao_one_e_ints/README.rst index bc8832c3..a87bdb17 100644 --- a/src/ao_one_e_ints/README.rst +++ b/src/ao_one_e_ints/README.rst @@ -6,9 +6,8 @@ All the one-electron integrals in the |AO| basis are here. The most important providers for usual quantum-chemistry calculation are: -* `ao_kinetic_integral` which are the kinetic operator integrals on the |AO| basis (see :file:`kin_ao_ints.irp.f`) -* `ao_nucl_elec_integral` which are the nuclear-elctron operator integrals on the |AO| basis (see :file:`pot_ao_ints.irp.f`) -* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis (see :file:`ao_mono_ints.irp.f`) +* `ao_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis +* `ao_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis +* `ao_one_e_integrals` which are the the h_core operator integrals on the |AO| basis -Note that you can find other interesting integrals related to the position operator in :file:`spread_dipole_ao.irp.f`. diff --git a/src/bitmask/bitmask_cas_routines.irp.f b/src/bitmask/bitmask_cas_routines.irp.f index 568005b3..c0c8cd11 100644 --- a/src/bitmask/bitmask_cas_routines.irp.f +++ b/src/bitmask/bitmask_cas_routines.irp.f @@ -2,6 +2,29 @@ use bitmasks integer function number_of_holes(key_in) BEGIN_DOC ! Function that returns the number of holes in the inact space + ! +! popcnt( +! xor( +! iand( +! reunion_of_core_inact_bitmask(1,1), +! xor( +! key_in(1,1), +! iand( +! key_in(1,1), +! cas_bitmask(1,1,1)) +! ) +! ), +! reunion_of_core_inact_bitmask(1,1)) ) +! +! (key_in && cas_bitmask) +! +---------------------+ +! electrons in cas xor key_in +! +---------------------------------+ +! electrons outside of cas && reunion_of_core_inact_bitmask +! +------------------------------------------------------------------+ +! electrons in the core/inact space xor reunion_of_core_inact_bitmask +! +---------------------------------------------------------------------------------+ +! holes END_DOC implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) @@ -78,22 +101,6 @@ integer function number_of_holes(key_in) + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )& + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,1), xor(key_in(7,1),iand(key_in(7,1),cas_bitmask(7,1,1)))), reunion_of_core_inact_bitmask(7,1)) )& + popcnt( xor( iand(reunion_of_core_inact_bitmask(7,2), xor(key_in(7,2),iand(key_in(7,2),cas_bitmask(7,2,1)))), reunion_of_core_inact_bitmask(7,2)) ) - else if(N_int == 8)then - number_of_holes = number_of_holes & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,1), xor(key_in(1,1),iand(key_in(1,1),cas_bitmask(1,1,1)))), reunion_of_core_inact_bitmask(1,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(1,2), xor(key_in(1,2),iand(key_in(1,2),cas_bitmask(1,2,1)))), reunion_of_core_inact_bitmask(1,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,1), xor(key_in(2,1),iand(key_in(2,1),cas_bitmask(2,1,1)))), reunion_of_core_inact_bitmask(2,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(2,2), xor(key_in(2,2),iand(key_in(2,2),cas_bitmask(2,2,1)))), reunion_of_core_inact_bitmask(2,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,1), xor(key_in(3,1),iand(key_in(3,1),cas_bitmask(3,1,1)))), reunion_of_core_inact_bitmask(3,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(3,2), xor(key_in(3,2),iand(key_in(3,2),cas_bitmask(3,2,1)))), reunion_of_core_inact_bitmask(3,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,1), xor(key_in(4,1),iand(key_in(4,1),cas_bitmask(4,1,1)))), reunion_of_core_inact_bitmask(4,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(4,2), xor(key_in(4,2),iand(key_in(4,2),cas_bitmask(4,2,1)))), reunion_of_core_inact_bitmask(4,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,1), xor(key_in(5,1),iand(key_in(5,1),cas_bitmask(5,1,1)))), reunion_of_core_inact_bitmask(5,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(5,2), xor(key_in(5,2),iand(key_in(5,2),cas_bitmask(5,2,1)))), reunion_of_core_inact_bitmask(5,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,1), xor(key_in(6,1),iand(key_in(6,1),cas_bitmask(6,1,1)))), reunion_of_core_inact_bitmask(6,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(6,2), xor(key_in(6,2),iand(key_in(6,2),cas_bitmask(6,2,1)))), reunion_of_core_inact_bitmask(6,2)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,1), xor(key_in(8,1),iand(key_in(8,1),cas_bitmask(8,1,1)))), reunion_of_core_inact_bitmask(8,1)) )& - + popcnt( xor( iand(reunion_of_core_inact_bitmask(8,2), xor(key_in(8,2),iand(key_in(8,2),cas_bitmask(8,2,1)))), reunion_of_core_inact_bitmask(8,2)) ) else do i = 1, N_int number_of_holes = number_of_holes & @@ -108,7 +115,7 @@ integer function number_of_holes(key_in) ) & ) & ), reunion_of_core_inact_bitmask(i,1)) ) & - + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,1), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,1,1)))), reunion_of_core_inact_bitmask(i,1)) ) + + popcnt( xor( iand(reunion_of_core_inact_bitmask(i,2), xor(key_in(i,2),iand(key_in(i,2),cas_bitmask(i,2,1)))), reunion_of_core_inact_bitmask(i,2)) ) enddo endif end diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index cb1bc83f..b1adfa20 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -141,8 +141,8 @@ subroutine run_pt2_slave_small(thread,iproc,energy) b%cur=0 ! ! Try to adjust n_tasks around nproc/2 seconds per job -! n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0))) - n_tasks = 1 + n_tasks = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0))) +! n_tasks = 1 end do integer, external :: disconnect_from_taskserver diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index ac889794..70ad543f 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -85,7 +85,7 @@ subroutine run_selection_slave(thread,iproc,energy) if(ctask > 0) then call sort_selection_buffer(buf) ! call merge_selection_buffers(buf,buf2) -print *, task_id(1), pt2(1), buf%cur, ctask +!print *, task_id(1), pt2(1), buf%cur, ctask call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask) ! buf%mini = buf2%mini pt2(:) = 0d0 diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 95128969..93a91933 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -66,6 +66,12 @@ interface: ezfio doc: Number of determinants in the current wave function type: integer +[n_det_qp_edit] +interface: ezfio +doc: Number of determinants to print in qp_edit +type: integer +interface: ezfio + [psi_coef] interface: ezfio doc: Coefficients of the wave function @@ -78,6 +84,18 @@ doc: Determinants of the variational space type: integer*8 size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det) +[psi_coef_qp_edit] +interface: ezfio +doc: Coefficients of the wave function +type: double precision +size: (determinants.n_det_qp_edit,determinants.n_states) + +[psi_det_qp_edit] +interface: ezfio +doc: Determinants of the variational space +type: integer*8 +size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det_qp_edit) + [expected_s2] interface: ezfio doc: Expected value of |S^2| diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f index ddb9ae0f..cec87901 100644 --- a/src/determinants/create_excitations.irp.f +++ b/src/determinants/create_excitations.irp.f @@ -12,6 +12,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok) integer(bit_kind), intent(inout) :: key_in(N_int,2) integer, intent(out) :: i_ok integer :: k,j,i + integer(bit_kind) :: mask use bitmasks ASSERT (i_hole > 0 ) ASSERT (i_particle <= mo_num) @@ -19,31 +20,66 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok) ! hole k = shiftr(i_hole-1,bit_kind_shift)+1 j = i_hole-shiftl(k-1,bit_kind_shift)-1 + mask = ibset(0_bit_kind,j) ! check whether position j is occupied - if (ibits(key_in(k,ispin),j,1).eq.1) then + if (iand(key_in(k,ispin),mask) /= 0_bit_kind) then key_in(k,ispin) = ibclr(key_in(k,ispin),j) else i_ok= -1 + return end if ! particle k = shiftr(i_particle-1,bit_kind_shift)+1 j = i_particle-shiftl(k-1,bit_kind_shift)-1 - key_in(k,ispin) = ibset(key_in(k,ispin),j) + mask = ibset(0_bit_kind,j) + if (iand(key_in(k,ispin),mask) == 0_bit_kind) then + key_in(k,ispin) = ibset(key_in(k,ispin),j) + else + i_ok= -1 + return + end if - integer :: n_elec_tmp - n_elec_tmp = 0 - do i = 1, N_int - n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2)) - enddo - if(n_elec_tmp .ne. elec_num)then - !print*, n_elec_tmp,elec_num - !call debug_det(key_in,N_int) - i_ok = -1 - endif +! integer :: n_elec_tmp +! n_elec_tmp = 0 +! do i = 1, N_int +! n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2)) +! enddo +! if(n_elec_tmp .ne. elec_num)then +! print*, n_elec_tmp,elec_num +! call debug_det(key_in,N_int) +! stop -1 +! endif end +subroutine build_singly_excited_wavefunction(i_hole,i_particle,ispin,det_out,coef_out) + implicit none + BEGIN_DOC + ! Applies the single excitation operator : a^{dager}_(i_particle) a_(i_hole) of + ! spin = ispin to the current wave function (psi_det, psi_coef) + END_DOC + integer, intent(in) :: i_hole,i_particle,ispin + integer(bit_kind), intent(out) :: det_out(N_int,2,N_det) + double precision, intent(out) :: coef_out(N_det,N_states) + + integer :: k + integer :: i_ok + double precision :: phase + do k=1,N_det + coef_out(k,:) = psi_coef(k,:) + det_out(:,:,k) = psi_det(:,:,k) + call do_single_excitation(det_out(1,1,k),i_hole,i_particle,ispin,i_ok) + if (i_ok == 1) then + call get_phase(psi_det(1,1,k), det_out(1,1,k),phase,N_int) + coef_out(k,:) = phase * coef_out(k,:) + else + coef_out(k,:) = 0.d0 + det_out(:,:,k) = psi_det(:,:,k) + endif + enddo +end + logical function is_spin_flip_possible(key_in,i_flip,ispin) implicit none BEGIN_DOC diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index e4f76bca..a930d70b 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -269,7 +269,6 @@ subroutine save_natural_mos END_DOC call set_natural_mos call save_mos - end diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index fbb3a813..71ee3d89 100644 --- a/src/determinants/determinants.irp.f +++ b/src/determinants/determinants.irp.f @@ -44,6 +44,16 @@ BEGIN_PROVIDER [ integer, N_det ] ASSERT (N_det > 0) END_PROVIDER + +BEGIN_PROVIDER [ integer, N_det_qp_edit ] + implicit none + BEGIN_DOC +! Number of determinants to print in qp_edit + END_DOC + + N_det_qp_edit = min(N_det,10000) +END_PROVIDER + BEGIN_PROVIDER [integer, max_degree_exc] implicit none integer :: i,degree @@ -476,7 +486,7 @@ subroutine save_wavefunction endif if (mpi_master) then call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) - endif + endif end @@ -504,12 +514,16 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) integer*8, allocatable :: psi_det_save(:,:,:) double precision, allocatable :: psi_coef_save(:,:) - integer :: i,j,k + double precision :: accu_norm + integer :: i,j,k, ndet_qp_edit if (mpi_master) then + ndet_qp_edit = min(ndet,N_det_qp_edit) + call ezfio_set_determinants_N_int(N_int) call ezfio_set_determinants_bit_kind(bit_kind) call ezfio_set_determinants_N_det(ndet) + call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit) call ezfio_set_determinants_n_states(nstates) call ezfio_set_determinants_mo_label(mo_label) @@ -522,10 +536,10 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) enddo enddo call ezfio_set_determinants_psi_det(psi_det_save) + call ezfio_set_determinants_psi_det_qp_edit(psi_det_save) deallocate (psi_det_save) allocate (psi_coef_save(ndet,nstates)) - double precision :: accu_norm do k=1,nstates do i=1,ndet psi_coef_save(i,k) = psicoef(i,k) @@ -535,6 +549,18 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) call ezfio_set_determinants_psi_coef(psi_coef_save) deallocate (psi_coef_save) + + allocate (psi_coef_save(ndet_qp_edit,nstates)) + do k=1,nstates + do i=1,ndet_qp_edit + psi_coef_save(i,k) = psicoef(i,k) + enddo + call normalize(psi_coef_save(1,k),ndet_qp_edit) + enddo + + call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save) + deallocate (psi_coef_save) + call write_int(6,ndet,'Saved determinants') endif end @@ -559,54 +585,80 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde integer :: N_int2 equivalence (det_8, det_bk) - integer :: i,k + integer :: i,j,k, ndet_qp_edit - call ezfio_set_determinants_N_int(N_int) - call ezfio_set_determinants_bit_kind(bit_kind) - call ezfio_set_determinants_N_det(ndetsave) - call ezfio_set_determinants_n_states(nstates) - call ezfio_set_determinants_mo_label(mo_label) + if (mpi_master) then + ndet_qp_edit = min(ndetsave,N_det_qp_edit) + call ezfio_set_determinants_N_int(N_int) + call ezfio_set_determinants_bit_kind(bit_kind) + call ezfio_set_determinants_N_det(ndetsave) + call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit) + call ezfio_set_determinants_n_states(nstates) + call ezfio_set_determinants_mo_label(mo_label) - N_int2 = (N_int*bit_kind)/8 - allocate (psi_det_save(N_int2,2,ndetsave)) - do i=1,ndetsave - do k=1,N_int - det_bk(k) = psidet(k,1,index_det_save(i)) - enddo - do k=1,N_int2 - psi_det_save(k,1,i) = det_8(k) - enddo - do k=1,N_int - det_bk(k) = psidet(k,2,index_det_save(i)) - enddo - do k=1,N_int2 - psi_det_save(k,2,i) = det_8(k) - enddo - enddo - call ezfio_set_determinants_psi_det(psi_det_save) - deallocate (psi_det_save) - - allocate (psi_coef_save(ndetsave,nstates)) - double precision :: accu_norm(nstates) - accu_norm = 0.d0 - do k=1,nstates + N_int2 = (N_int*bit_kind)/8 + allocate (psi_det_save(N_int2,2,ndetsave)) do i=1,ndetsave - accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k) - psi_coef_save(i,k) = psicoef(index_det_save(i),k) + do k=1,N_int + det_bk(k) = psidet(k,1,index_det_save(i)) + enddo + do k=1,N_int2 + psi_det_save(k,1,i) = det_8(k) + enddo + do k=1,N_int + det_bk(k) = psidet(k,2,index_det_save(i)) + enddo + do k=1,N_int2 + psi_det_save(k,2,i) = det_8(k) + enddo enddo - enddo - do k = 1, nstates - accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) - enddo - do k=1,nstates - do i=1,ndetsave - psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) - enddo - enddo + call ezfio_set_determinants_psi_det(psi_det_save) + call ezfio_set_determinants_psi_det_qp_edit(psi_det_save) + deallocate (psi_det_save) - call ezfio_set_determinants_psi_coef(psi_coef_save) - call write_int(6,ndet,'Saved determinants') - deallocate (psi_coef_save) + allocate (psi_coef_save(ndetsave,nstates)) + double precision :: accu_norm(nstates) + accu_norm = 0.d0 + do k=1,nstates + do i=1,ndetsave + accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k) + psi_coef_save(i,k) = psicoef(index_det_save(i),k) + enddo + enddo + do k = 1, nstates + accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) + enddo + do k=1,nstates + do i=1,ndetsave + psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) + enddo + enddo + + call ezfio_set_determinants_psi_coef(psi_coef_save) + deallocate (psi_coef_save) + + allocate (psi_coef_save(ndet_qp_edit,nstates)) + accu_norm = 0.d0 + do k=1,nstates + do i=1,ndet_qp_edit + accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k) + psi_coef_save(i,k) = psicoef(index_det_save(i),k) + enddo + enddo + do k = 1, nstates + accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) + enddo + do k=1,nstates + do i=1,ndet_qp_edit + psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) + enddo + enddo + + call ezfio_set_determinants_psi_coef(psi_coef_save) + deallocate (psi_coef_save) + + call write_int(6,ndet,'Saved determinants') + endif end diff --git a/src/determinants/occ_pattern.irp.f b/src/determinants/occ_pattern.irp.f index d5f458a0..5f37b289 100644 --- a/src/determinants/occ_pattern.irp.f +++ b/src/determinants/occ_pattern.irp.f @@ -53,7 +53,17 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint) use bitmasks implicit none BEGIN_DOC - ! Generate all possible determinants for a give occ_pattern + ! Generate all possible determinants for a given occ_pattern + ! + ! Input : + ! o : occupation pattern : (doubly occupied, singly occupied) + ! sze : Number of produced determinants, computed by `occ_pattern_to_dets_size` + ! n_alpha : Number of $\alpha$ electrons + ! Nint : N_int + ! + ! Output: + ! d : determinants + ! END_DOC integer ,intent(in) :: Nint integer ,intent(in) :: n_alpha ! Number of alpha electrons diff --git a/src/dft_utils_in_r/kin_dens.irp.f b/src/dft_utils_in_r/kin_dens.irp.f new file mode 100644 index 00000000..008d6b2d --- /dev/null +++ b/src/dft_utils_in_r/kin_dens.irp.f @@ -0,0 +1,15 @@ +BEGIN_PROVIDER [double precision, kinetic_density_generalized, (n_points_final_grid)] + implicit none + integer :: i,j,m,i_point + kinetic_density_generalized = 0.d0 + do i_point = 1, n_points_final_grid + do i = 1, mo_num + do j = 1, mo_num + do m = 1, 3 + kinetic_density_generalized(i_point) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i_point) * mos_grad_in_r_array_tranp(m,i,i_point) * one_e_dm_mo_for_dft(j,i,1) + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 9c92481c..60cd59f2 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -22,6 +22,7 @@ BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_num,n_points_final_grid,3)] +&BEGIN_PROVIDER[double precision, mos_grad_in_r_array_tranp,(3,mo_num,n_points_final_grid)] implicit none BEGIN_DOC ! mos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith mo on the jth grid point @@ -35,8 +36,37 @@ do m=1,3 call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_num) enddo + integer :: i,j + do i = 1, n_points_final_grid + do j = 1, mo_num + do m = 1, 3 + mos_grad_in_r_array_tranp(m,j,i) = mos_grad_in_r_array(j,i,m) + enddo + enddo + enddo END_PROVIDER + BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)] +&BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)] + implicit none + integer :: i,m,j + alpha_dens_kin_in_r = 0.d0 + beta_dens_kin_in_r = 0.d0 + do i = 1, n_points_final_grid + do j = 1, elec_alpha_num + do m = 1, 3 + alpha_dens_kin_in_r(i) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0 + enddo + enddo + do j = 1, elec_beta_num + do m = 1, 3 + beta_dens_kin_in_r(i) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i)**2.d0 + enddo + enddo + enddo + + END_PROVIDER + BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_num,n_points_final_grid,3)] implicit none BEGIN_DOC diff --git a/src/dft_utils_one_e/ec_lyp.irp.f b/src/dft_utils_one_e/ec_lyp.irp.f index 0e3daa65..22d15a9c 100644 --- a/src/dft_utils_one_e/ec_lyp.irp.f +++ b/src/dft_utils_one_e/ec_lyp.irp.f @@ -35,11 +35,8 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr include 'constants.include.F' ! Input variables - double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2 - ! Local variables - double precision :: a,b,c,d,c_f,omega,delta double precision :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom double precision :: thr,huge_num,rho_inv @@ -47,8 +44,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr double precision :: tmp1,tmp2,tmp3,tmp4 double precision :: big1,big2,big3 -! Output variables - ! Constants of the LYP correlation functional @@ -57,15 +52,25 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr c = 0.2533d0 d = 0.349d0 - thr = 1d-10 + ec_lyp_88 = 0.d0 + + thr = 1d-15 huge_num = 1.d0/thr + if(dabs(rho_a).lt.thr)then + return + endif + + if(dabs(rho_b).lt.thr)then + return + endif if(rho.lt.0.d0)then print*,'pb !! rho.lt.0.d0' stop endif - rho_13 = rho**(1d0/3d0) - rho_113 = rho**(11d0/3d0) + + rho_13 = rho**(1.d0/3.d0) + rho_113 = rho**(11.d0/3.d0) if(dabs(rho_13) < thr) then rho_inv_13 = huge_num @@ -76,13 +81,13 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr if (dabs(rho_113) < thr) then rho_inv_113 = huge_num else - rho_inv_113 = 1d0/rho_113 + rho_inv_113 = 1.d0/rho_113 endif if (dabs(rho) < thr) then rho_inv = huge_num else - rho_inv = 1d0/rho + rho_inv = 1.d0/rho endif ! Useful quantities to predefine @@ -90,21 +95,21 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr denom = 1d0/(1d0 + d*rho_inv_13) omega = rho_inv_113*exp(-c*rho_inv_13)*denom delta = c*rho_inv_13 + d*rho_inv_13*denom - c_f = 0.3d0*(3d0*pi*pi)**(2d0/3d0) + c_f = 0.3d0*(3.d0*pi*pi)**(2.d0/3.d0) rho_2 = rho *rho rho_a_2 = rho_a*rho_a rho_b_2 = rho_b*rho_b - cst_2_113 = 2d0**(11d0/3d0) - cst_8_3 = 8d0/3d0 + cst_2_113 = 2.d0**(11.d0/3.d0) + cst_8_3 = 8.d0/3.d0 ! first term in the equation (2) of Preuss CPL, 1989 - big1 = 4d0*denom*rho_a*rho_b*rho_inv + big1 = 4.d0*denom*rho_a*rho_b*rho_inv tmp1 = cst_2_113*c_f*(rho_a**cst_8_3 + rho_b**cst_8_3) - tmp2 = (47d0/18d0 - 7d0/18d0*delta)*grad_rho_2 + tmp2 = (47.d0/18.d0 - 7.d0/18.d0*delta)*grad_rho_2 tmp3 = - (5d0/2d0 - 1.d0/18d0*delta)*(grad_rho_a_2 + grad_rho_b_2) tmp4 = - (delta - 11d0)/9d0*(rho_a*rho_inv*grad_rho_a_2 + rho_b*rho_inv*grad_rho_b_2) big2 = rho_a*rho_b*(tmp1 + tmp2 + tmp3 + tmp4) @@ -114,7 +119,6 @@ double precision function ec_lyp_88(rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,gr tmp3 = grad_rho_a_2*(2d0/3d0*rho_2 - rho_b_2) big3 = tmp1 + tmp2 + tmp3 - ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3 end diff --git a/src/dft_utils_one_e/ec_scan.irp.f b/src/dft_utils_one_e/ec_scan.irp.f new file mode 100644 index 00000000..4807b89f --- /dev/null +++ b/src/dft_utils_one_e/ec_scan.irp.f @@ -0,0 +1,100 @@ +double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2) + include 'constants.include.F' + implicit none + double precision, intent(in) :: rho_a,rho_b,tau,grad_rho_2 + double precision :: cst_13,cst_23,cst_43,cst_53,rho_inv,cst_18,cst_3pi2 + double precision :: thr,nup,ndo,xi,s,spin_d,drho,drho2,rho,inv_1alph,e_c_lsda1,h0 + double precision :: rs,t_w,t_unif,ds_xi,alpha,fc_alpha,step_f,cst_1alph,beta_inf + double precision :: c_1c,c_2c,d_c,e_c_ldsa1,h1,phi,t,beta_rs,gama,a,w_1,g_at2,phi_3,e_c_1 + double precision :: b_1c,b_2c,b_3c,dx_xi,gc_xi,e_c_lsda0,w_0,g_inf,cx_xi,x_inf,f0,e_c_0 + thr = 1.d-12 + nup = max(rho_a,thr) + ndo = max(rho_b,thr) + rho = nup + ndo + ec_scan = 0.d0 + if((rho).lt.thr)return + ! constants ... + rho_inv = 1.d0/rho + cst_13 = 1.d0/3.d0 + cst_23 = 2.d0 * cst_13 + cst_43 = 4.d0 * cst_13 + cst_53 = 5.d0 * cst_13 + cst_18 = 1.d0/8.d0 + cst_3pi2 = 3.d0 * pi*pi + drho2 = max(grad_rho_2,thr) + drho = dsqrt(drho2) + if((nup-ndo).gt.0.d0)then + spin_d = max(nup-ndo,thr) + else + spin_d = min(nup-ndo,-thr) + endif + c_1c = 0.64d0 + c_2c = 1.5d0 + d_c = 0.7d0 + b_1c = 0.0285764d0 + b_2c = 0.0889d0 + b_3c = 0.125541d0 + gama = 0.031091d0 + ! correlation energy lsda1 + call ec_only_lda_sr(0.d0,nup,ndo,e_c_lsda1) + + xi = spin_d/rho + rs = (cst_43 * pi * rho)**(-cst_13) + s = drho/( 2.d0 * cst_3pi2**(cst_13) * rho**cst_43 ) + t_w = drho2 * cst_18 * rho_inv + ds_xi = 0.5d0 * ( (1.d0+xi)**cst_53 + (1.d0 - xi)**cst_53) + t_unif = 0.3d0 * (cst_3pi2)**cst_23 * rho**cst_53*ds_xi + t_unif = max(t_unif,thr) + alpha = (tau - t_w)/t_unif + cst_1alph= 1.d0 - alpha + if(cst_1alph.gt.0.d0)then + cst_1alph= max(cst_1alph,thr) + else + cst_1alph= min(cst_1alph,-thr) + endif + inv_1alph= 1.d0/cst_1alph + phi = 0.5d0 * ( (1.d0+xi)**cst_23 + (1.d0 - xi)**cst_23) + phi_3 = phi*phi*phi + t = (cst_3pi2/16.d0)**cst_13 * s / (phi * rs**0.5d0) + w_1 = dexp(-e_c_lsda1/(gama * phi_3)) - 1.d0 + a = beta_rs(rs) /(gama * w_1) + g_at2 = 1.d0/(1.d0 + 4.d0 * a*t*t)**0.25d0 + h1 = gama * phi_3 * dlog(1.d0 + w_1 * (1.d0 - g_at2)) + ! interpolation function + fc_alpha = dexp(-c_1c * alpha * inv_1alph) * step_f(cst_1alph) - d_c * dexp(c_2c * inv_1alph) * step_f(-cst_1alph) + ! first part of the correlation energy + e_c_1 = e_c_lsda1 + h1 + + dx_xi = 0.5d0 * ( (1.d0+xi)**cst_43 + (1.d0 - xi)**cst_43) + gc_xi = (1.d0 - 2.3631d0 * (dx_xi - 1.d0) ) * (1.d0 - xi**12.d0) + e_c_lsda0= - b_1c / (1.d0 + b_2c * rs**0.5d0 + b_3c * rs) + w_0 = dexp(-e_c_lsda0/b_1c) - 1.d0 + beta_inf = 0.066725d0 * 0.1d0 / 0.1778d0 + cx_xi = -3.d0/(4.d0*pi) * (9.d0 * pi/4.d0)**cst_13 * dx_xi + + x_inf = 0.128026d0 + f0 = -0.9d0 + g_inf = 1.d0/(1.d0 + 4.d0 * x_inf * s*s)**0.25d0 + + h0 = b_1c * dlog(1.d0 + w_0 * (1.d0 - g_inf)) + e_c_0 = (e_c_lsda0 + h0) * gc_xi + + ec_scan = e_c_1 + fc_alpha * (e_c_0 - e_c_1) +end + +double precision function step_f(x) + implicit none + double precision, intent(in) :: x + if(x.lt.0.d0)then + step_f = 0.d0 + else + step_f = 1.d0 + endif +end + +double precision function beta_rs(rs) + implicit none + double precision, intent(in) ::rs + beta_rs = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs) + +end diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index f3b7db2c..812cd3d4 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -36,24 +36,24 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file f2.ezfio qp set_frozen_core - run_stoch -199.30496 1.e-4 + run_stoch -199.30486 1.e-4 } @test "NH3" { # 10.6657s qp set_file nh3.ezfio qp set_mo_class --core="[1-4]" --act="[5-72]" - run -56.244753429144986 1.e-5 + run -56.244753429144986 1.e-4 } @test "DHNO" { # 11.4721s qp set_file dhno.ezfio qp set_mo_class --core="[1-7]" --act="[8-64]" - run -130.459020029816 1.e-5 + run -130.459020029816 1.e-4 } @test "HCO" { # 12.2868s qp set_file hco.ezfio - run -113.297494345682 2.e-05 + run -113.297494345682 1.e-4 } @test "H2O2" { # 12.9214s @@ -65,82 +65,82 @@ function run_stoch() { @test "HBO" { # 13.3144s [[ -n $TRAVIS ]] && skip qp set_file hbo.ezfio - run -100.214185815312 1.e-5 + run -100.212829869715 1.e-4 } @test "H2O" { # 11.3727s [[ -n $TRAVIS ]] && skip qp set_file h2o.ezfio - run -76.2359268957699 2.e-5 + run -76.2359268957699 1.e-4 } @test "ClO" { # 13.3755s [[ -n $TRAVIS ]] && skip qp set_file clo.ezfio - run -534.546005867797 5.e-5 + run -534.545881614967 1.e-4 } @test "SO" { # 13.4952s [[ -n $TRAVIS ]] && skip qp set_file so.ezfio - run -26.0124797722154 1.e-5 + run -26.0158153138924 1.e-4 } @test "H2S" { # 13.6745s [[ -n $TRAVIS ]] && skip qp set_file h2s.ezfio - run -398.859480581924 1.e-5 + run -398.859168655255 1.e-4 } @test "OH" { # 13.865s [[ -n $TRAVIS ]] && skip qp set_file oh.ezfio - run -75.6119887538831 1.e-05 + run -75.6120779012574 1.e-4 } @test "SiH2_3B1" { # 13.938ss [[ -n $TRAVIS ]] && skip qp set_file sih2_3b1.ezfio - run -290.017539006762 1.e-5 + run -290.017539006762 1.e-4 } @test "H3COH" { # 14.7299s [[ -n $TRAVIS ]] && skip qp set_file h3coh.ezfio - run -115.205054063687 1.e-5 + run -115.205941463667 1.e-4 } @test "SiH3" { # 15.99s [[ -n $TRAVIS ]] && skip qp set_file sih3.ezfio - run -5.57269434557089 2.e-05 + run -5.57241217753818 1.e-4 } @test "CH4" { # 16.1612s [[ -n $TRAVIS ]] && skip qp set_file ch4.ezfio qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" - run -40.2409059687324 2.e-5 + run -40.2409678239136 1.e-4 } @test "ClF" { # 16.8864s [[ -n $TRAVIS ]] && skip qp set_file clf.ezfio - run -559.170406471496 1.e-5 + run -559.170272077166 1.e-4 } @test "SO2" { # 17.5645s [[ -n $TRAVIS ]] && skip qp set_file so2.ezfio qp set_mo_class --core="[1-8]" --act="[9-87]" - run -41.5746738713298 5.e-5 + run -41.5746738713298 1.e-4 } @test "C2H2" { # 17.6827s [[ -n $TRAVIS ]] && skip qp set_file c2h2.ezfio qp set_mo_class --act="[1-30]" --del="[31-36]" - run -12.3670840202635 2.e-5 + run -12.3656179738175 1.e-4 } @test "N2" { # 18.0198s @@ -154,14 +154,14 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file n2h4.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]" - run -111.367234092521 2.e-5 + run -111.367332681559 1.e-4 } @test "CO2" { # 21.1748s [[ -n $TRAVIS ]] && skip qp set_file co2.ezfio qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]" - run -187.969676381867 1.e-5 + run -187.968599504402 1.e-4 } @@ -169,13 +169,13 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file cu_nh3_4_2plus.ezfio qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]" - run -1862.98610987882 1.e-05 + run -1862.98614665139 1.e-04 } @test "HCN" { # 20.3273s [[ -n $TRAVIS ]] && skip qp set_file hcn.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]" - run -93.0799328685679 2.e-5 + run -93.0728641601823 1.e-4 } diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index 7bcbf013..ebdcfead 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -21,7 +21,7 @@ type: double precision size: (nuclei.nucl_num,3) interface: ezfio -[disk_access_nuclear_repulsion] +[io_nuclear_repulsion] doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ] type: Disk_access interface: ezfio,provider,ocaml diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index c797ff30..88ee24e8 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] END_DOC PROVIDE mpi_master nucl_coord nucl_charge nucl_num - if (disk_access_nuclear_repulsion.EQ.'Read') then + if (io_nuclear_repulsion == 'Read') then logical :: has if (mpi_master) then @@ -194,7 +194,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] call write_time(6) call write_double(6,nuclear_repulsion,'Nuclear repulsion energy') - if (disk_access_nuclear_repulsion.EQ.'Write') then + if (io_nuclear_repulsion == 'Write') then if (mpi_master) then call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) endif diff --git a/src/tools/print_ci_vectors.irp.f b/src/tools/print_ci_vectors.irp.f index 9ba06d9a..97dfdc0b 100644 --- a/src/tools/print_ci_vectors.irp.f +++ b/src/tools/print_ci_vectors.irp.f @@ -24,6 +24,7 @@ subroutine routine implicit none integer :: i,k integer :: degree + call print_energy_components do i = 1, N_det print *, 'Determinant ', i call debug_det(psi_det(1,1,i),N_int)