10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-22 21:22:25 +02:00

Merge branch 'bugfix' of github.com:QuantumPackage/qp2 into bugfix

This commit is contained in:
Anthony Scemama 2019-10-21 14:13:32 +02:00
commit 9ef3e904de
50 changed files with 1561 additions and 469 deletions

View File

@ -143,7 +143,7 @@ IRPF90
to Parameters (IRP) method. to Parameters (IRP) method.
* Download the latest version of IRPF90 * Download the latest version of IRPF90
here : `<https://github.com/scemama/irpf90/releases/latest>`_ and move here : `<https://gitlab.com/scemama/irpf90/-/archive/v1.7.2/irpf90-v1.7.2.tar.gz>`_ and move
the downloaded archive in the :file:`${QP_ROOT}/external` directory the downloaded archive in the :file:`${QP_ROOT}/external` directory
* Extract the archive and go into the :file:`irpf90-*` directory to run * Extract the archive and go into the :file:`irpf90-*` directory to run

View File

@ -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_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_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 potential_sr_xc_beta_ao_pbe --rename=potential_xc_beta_ao_sr_pbe
qp_name disk_access_nuclear_repulsion --rename=io_nuclear_repulsion

View File

@ -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 <address> 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["<address>"]
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)

29
configure vendored
View File

@ -290,6 +290,9 @@ EOF
# Special commands for Travis CI # Special commands for Travis CI
chmod +x "${QP_ROOT}"/external/opam_installer.sh chmod +x "${QP_ROOT}"/external/opam_installer.sh
rm --force ${QP_ROOT}/bin/opam rm --force ${QP_ROOT}/bin/opam
if [[ -n ${NO_CACHE} ]] ; then
rm -rf ${HOME}/.opam
fi
export OPAMROOT=${HOME}/.opam export OPAMROOT=${HOME}/.opam
cat << EOF | bash ${QP_ROOT}/external/opam_installer.sh --no-backup cat << EOF | bash ${QP_ROOT}/external/opam_installer.sh --no-backup
${QP_ROOT}/bin ${QP_ROOT}/bin
@ -310,21 +313,21 @@ EOF
else else
# Conventional commands # Conventional commands
execute << EOF execute << EOF
chmod +x "\${QP_ROOT}"/external/opam_installer.sh chmod +x "${QP_ROOT}"/external/opam_installer.sh
"\${QP_ROOT}"/external/opam_installer.sh --no-backup "${QP_ROOT}"/external/opam_installer.sh --no-backup
EOF EOF
execute << EOF execute << EOF
rm --force \${QP_ROOT}/bin/opam rm --force ${QP_ROOT}/bin/opam
export OPAMROOT=\${OPAMROOT:-\${QP_ROOT}/external/opam} export OPAMROOT=${OPAMROOT:-${QP_ROOT}/external/opam}
echo \${QP_ROOT}/bin \ echo ${QP_ROOT}/bin \
| sh \${QP_ROOT}/external/opam_installer.sh | sh ${QP_ROOT}/external/opam_installer.sh
EOF EOF
rm ${QP_ROOT}/external/opam_installer.sh rm ${QP_ROOT}/external/opam_installer.sh
# source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true # source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
# opam switch create ocaml-base-compiler.4.07.1 || exit 1 # opam switch create ocaml-base-compiler.4.07.1 || exit 1
opam init --verbose --yes --compiler=4.07.1 --disable-sandboxing opam init --verbose --yes --compiler=4.07.1 --disable-sandboxing
eval $(opam env) eval $(opam env)
EOF
execute << EOF execute << EOF
opam install -y \${OCAML_PACKAGES} || exit 1 opam install -y \${OCAML_PACKAGES} || exit 1
EOF EOF
@ -423,18 +426,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
fail fail
fi fi
BWRAP=$(find_exe bwrap)
if [[ ${BWRAP} = $(not_found) ]] ; then
error "Bubblewrap (bwrap) is not installed."
fail
fi
LIBCAP=$(find_lib -lcap) LIBCAP=$(find_lib -lcap)
if [[ ${LIBCAP} = $(not_found) ]] ; then if [[ ${LIBCAP} = $(not_found) ]] ; then
error "Libcap (libcap) is not installed." error "Libcap (libcap) is not installed."
fail fail
fi fi
BWRAP=$(find_exe bwrap)
if [[ ${BWRAP} = $(not_found) ]] ; then
error "Bubblewrap (bwrap) is not installed."
fail
fi
OPAM=$(find_exe opam) OPAM=$(find_exe opam)
if [[ ${OPAM} = $(not_found) ]] ; then if [[ ${OPAM} = $(not_found) ]] ; then
error "OPAM (ocaml) package manager is not installed." error "OPAM (ocaml) package manager is not installed."

View File

@ -8,6 +8,16 @@ S 1
S 1 S 1
1 0.0360000 1.0000000 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 LITHIUM
S 6 S 6
1 642.4189200 0.0021426 1 642.4189200 0.0021426

View File

@ -1,14 +1,27 @@
HYDROGEN HYDROGEN
S 3 S 3
1 18.7311370 0.03349460 1 0.1873113696E+02 0.3349460434E-01
2 2.8253937 0.23472695 2 0.2825394365E+01 0.2347269535E+00
3 0.6401217 0.81375733 3 0.6401216923E+00 0.8137573261E+00
S 1 S 1
1 0.1612778 1.0000000 1 0.1612777588E+00 1.0000000
S 1 S 1
1 0.0360000 1.0000000 1 0.3600000000E-01 0.1000000000E+01
P 1 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 LITHIUM
S 6 S 6

View File

@ -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 LITHIUM
S 6 S 6
1 642.4189200 0.0021426 1 642.4189200 0.0021426

View File

@ -14,6 +14,18 @@ P 1
P 1 P 1
1 0.3750000 1.0000000 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 LITHIUM
S 6 S 6
1 900.4600000 0.00228704 1 900.4600000 0.00228704

View File

@ -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 LITHIUM
S 6 S 6
1 900.4600000 0.00228704 1 900.4600000 0.00228704

View File

@ -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 LITHIUM
S 6 S 6
1 900.4600000 0.00228704 1 900.4600000 0.00228704

View File

@ -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 LITHIUM
S 6 S 6
1 642.4189200 0.0021426 1 642.4189200 0.0021426

Binary file not shown.

Before

Width:  |  Height:  |  Size: 5.9 MiB

After

Width:  |  Height:  |  Size: 351 KiB

BIN
data/qp2_hd.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 5.9 MiB

View File

@ -1,6 +1,16 @@
%%% ARXIV TO BE UPDATED %%% %%% 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, @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}}, title = {{Chemically Accurate Excitation Energies With Small Basis Sets}},
journal = {arXiv}, journal = {arXiv},
year = {2019}, year = {2019},
@ -9,29 +19,66 @@
url = {https://arxiv.org/abs/1907.01245} 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 %%%% 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, @article{Loos_2019,
doi = {10.1021/acs.jpclett.9b01176}, doi = {10.1021/acs.jpclett.9b01176},
url = {https://doi.org/10.1021%2Facs.jpclett.9b01176}, url = {https://doi.org/10.1021%2Facs.jpclett.9b01176},

View File

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

View File

@ -115,7 +115,7 @@ create an |EZFIO| database with the 6-31G basis set:
.. code:: bash .. 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 The EZFIO database now contains data relative to the nuclear coordinates
and the atomic basis set: and the atomic basis set:

139
man/qp_tunnel.1 Normal file
View File

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

View File

@ -34,7 +34,7 @@ level margin: \\n[rst2man-indent\\n[rst2man-indent-level]]
.INDENT 3.5 .INDENT 3.5
Rotates molecular orbitals i and j by combining them as 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 )$ and
$1/sqrt{2} ( phi_i \- phi_j )$. $1/sqrt{2} ( phi_i - phi_j )$.
.sp .sp
Needs: Needs:
.INDENT 0.0 .INDENT 0.0

55
man/test.1 Normal file
View File

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

View File

@ -1,4 +1,4 @@
PKG core ZMQ cryptokit PKG core zmq cryptokit
B _build/ B _build/

View File

@ -7,82 +7,61 @@ Type for bits strings
list of Bits 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 *) (* String representation *)
let to_string b = let to_string b =
let rec do_work accu = function let int64_to_string x =
| [] -> accu String.init 64 (fun i ->
| head :: tail -> if Int64.logand x @@ Int64.shift_left 1L i <> 0L then
let new_accu = (Bit.to_string head) ^ accu '+'
in do_work new_accu tail else
'-')
in in
do_work "" b Array.map int64_to_string b
|> Array.to_list
|> String.concat ""
let of_string ?(zero='0') ?(one='1') s = let of_string ?(zero='0') ?(one='1') s =
List.init (String.length s) (String.get s) let n_int = ( (String.length s - 1) lsr 6 ) + 1 in
|> List.rev_map ( fun c -> let result = Array.make n_int 0L in
if (c = zero) then Bit.Zero String.iteri (fun i c ->
else if (c = one) then Bit.One if c = one then
else (failwith ("Error in bitstring ") ) ) 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 = let of_string_mp = of_string ~zero:'-' ~one:'+'
List.init (String.length s) (String.get s)
|> List.rev_map (function
| '-' -> Bit.Zero
| '+' -> Bit.One
| _ -> failwith ("Error in bitstring ") )
(* Create a bit list from an int64 *) (* Create a bit list from an int64 *)
let of_int64 i = let of_int64 i = [| 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)
(* Create an int64 from a bit list *) (* Create an int64 from a bit list *)
let to_int64 l = let to_int64 = function
assert ( (List.length l) <= 64) ; | [| i |] -> i
let rec do_work accu = function | _ -> failwith "N_int > 1"
| [] -> 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)
(* 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 *) (* Create a bit list from an array of int64 *)
let of_int64_array l = external of_int64_array : int64 array -> t = "%identity"
Array.map of_int64 l external to_int64_array : t -> int64 array = "%identity"
|> Array.to_list
|> List.concat
(* Create a bit list from a list of int64 *)
let of_int64_list l =
Array.of_list l |> of_int64_array
(* Compute n_int *) (* 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 ) 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 *) (* Create an int64 list from a bit list *)
let to_int64_list l = let to_int64_list l =
let rec do_work accu buf counter = function to_int64_array l |> Array.to_list
| [] ->
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
(* 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 *) (* Create a bit list from a list of MO indices *)
let of_mo_number_list n_int l = let of_mo_number_list n_int l =
let n_int = N_int_number.to_int n_int in let result = zero n_int in
let length = n_int*64 in List.iter (fun j ->
let a = Array.make length (Bit.Zero) in let i = (MO_number.to_int j) - 1 in
List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l; let iint = i lsr 6 in (* i / 64 *)
Array.to_list a 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 to_mo_number_list l =
let a = Array.of_list l in let rec aux_one x shift accu = function
let mo_num = MO_number.get_max () in | -1 -> accu
let rec do_work accu = function | i -> if Int64.logand x (Int64.shift_left 1L i) <> 0L then
| 0 -> accu aux_one x shift ( (i+shift) ::accu) (i-1)
| i -> else
begin aux_one x shift accu (i-1)
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
in 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 *) (* logical operations on bit_list *)
let logical_operator2 op a b = let and_operator a b = Array.map2 Int64.logand a b
let rec do_work_binary result a b = let xor_operator a b = Array.map2 Int64.logxor a b
match a, b with let or_operator a b = Array.map2 Int64.logor a b
| [], [] -> result let not_operator b = Array.map Int64.lognot b
| [], _ | _ , [] -> 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 pop_sign =
let mask =
(Int64.pred (Int64.shift_left 1L 63))
in in
List.rev (do_work_binary [] a b) fun x -> Int64.logand mask x
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
let popcnt b = let popcnt b =
List.fold_left (fun accu -> function Array.fold_left (fun accu x ->
| Bit.One -> accu+1 if x >= 0L then
| Bit.Zero -> accu accu + (Z.popcount @@ Z.of_int64 x)
) 0 b else
accu + 1 + (Z.popcount @@ Z.of_int64 (pop_sign x))
) 0 b

View File

@ -1,4 +1,4 @@
type t = Bit.t list type t
(** The zero bit list *) (** The zero bit list *)
val zero : Qptypes.N_int_number.t -> t val zero : Qptypes.N_int_number.t -> t

View File

@ -3,7 +3,8 @@ open Sexplib.Std
type t = int64 array [@@deriving sexp] 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 = let to_alpha_beta x =
@ -24,19 +25,6 @@ let to_bitlist_couple x =
in (xa,xb) 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 = let of_int64_array ~n_int ~alpha ~beta x =
@ -47,38 +35,29 @@ let of_int64_array ~n_int ~alpha ~beta x =
in in
if ( (Bitlist.popcnt a) <> alpha) then if ( (Bitlist.popcnt a) <> alpha) then
begin 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 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; end;
if ( (Bitlist.popcnt b) <> beta ) then if ( (Bitlist.popcnt b) <> beta ) then
begin 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 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; end;
x 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 = let ba, bb =
Bitlist.to_int64_array xa , Bitlist.to_int64_array xa ,
Bitlist.to_int64_array xb 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 in
of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb]) of_int64_array ~n_int ~alpha ~beta (Array.concat [ba;bb])
let to_string ~mo_num x = let to_string ~mo_num x =
let (xa,xb) = to_bitlist_couple x in let (xa,xb) = to_bitlist_couple x in
[ " " ; bitlist_to_string ~mo_num xa ; "\n" ; [ " " ; Bitlist.to_string xa ; "\n" ;
" " ; bitlist_to_string ~mo_num xb ] " " ; Bitlist.to_string xb ]
|> String.concat "" |> String.concat ""

View File

@ -24,7 +24,7 @@ val to_alpha_beta : t -> (int64 array)*(int64 array)
val to_bitlist_couple : t -> Bitlist.t * Bitlist.t val to_bitlist_couple : t -> Bitlist.t * Bitlist.t
(** Create from a bit list *) (** 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 -> alpha:Qptypes.Elec_alpha_number.t ->
beta:Qptypes.Elec_beta_number.t -> beta:Qptypes.Elec_beta_number.t ->
Bitlist.t * Bitlist.t -> t Bitlist.t * Bitlist.t -> t

View File

@ -7,14 +7,14 @@ module Determinants_by_hand : sig
{ n_int : N_int_number.t; { n_int : N_int_number.t;
bit_kind : Bit_kind.t; bit_kind : Bit_kind.t;
n_det : Det_number.t; n_det : Det_number.t;
n_det_qp_edit : Det_number.t;
n_states : States_number.t; n_states : States_number.t;
expected_s2 : Positive_float.t; expected_s2 : Positive_float.t;
psi_coef : Det_coef.t array; psi_coef : Det_coef.t array;
psi_det : Determinant.t array; psi_det : Determinant.t array;
state_average_weight : Positive_float.t array; state_average_weight : Positive_float.t array;
} [@@deriving sexp] } [@@deriving sexp]
val read : unit -> t val read : ?full:bool -> unit -> t option
val read_maybe : unit -> t option
val write : t -> unit val write : t -> unit
val to_string : t -> string val to_string : t -> string
val to_rst : t -> Rst_string.t val to_rst : t -> Rst_string.t
@ -28,6 +28,7 @@ end = struct
{ n_int : N_int_number.t; { n_int : N_int_number.t;
bit_kind : Bit_kind.t; bit_kind : Bit_kind.t;
n_det : Det_number.t; n_det : Det_number.t;
n_det_qp_edit : Det_number.t;
n_states : States_number.t; n_states : States_number.t;
expected_s2 : Positive_float.t; expected_s2 : Positive_float.t;
psi_coef : Det_coef.t array; psi_coef : Det_coef.t array;
@ -38,8 +39,6 @@ end = struct
let get_default = Qpackage.get_ezfio_default "determinants";; let get_default = Qpackage.get_ezfio_default "determinants";;
let n_det_read_max = 10_000 ;;
let read_n_int () = let read_n_int () =
if not (Ezfio.has_determinants_n_int()) then if not (Ezfio.has_determinants_n_int()) then
Ezfio.get_mo_basis_mo_num () Ezfio.get_mo_basis_mo_num ()
@ -80,11 +79,27 @@ end = struct
|> Det_number.of_int |> 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 = let write_n_det n =
Det_number.to_int n Det_number.to_int n
|> Ezfio.set_determinants_n_det |> 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 () = let read_n_states () =
if not (Ezfio.has_determinants_n_states ()) then if not (Ezfio.has_determinants_n_states ()) then
Ezfio.set_determinants_n_states 1 Ezfio.set_determinants_n_states 1
@ -178,7 +193,7 @@ end = struct
|> Ezfio.set_determinants_expected_s2 |> Ezfio.set_determinants_expected_s2
;; ;;
let read_psi_coef () = let read_psi_coef ~read_only () =
if not (Ezfio.has_determinants_psi_coef ()) then if not (Ezfio.has_determinants_psi_coef ()) then
begin begin
let n_states = let n_states =
@ -189,24 +204,33 @@ end = struct
~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. )) ~data:(List.init n_states (fun i -> if (i=0) then 1. else 0. ))
|> Ezfio.set_determinants_psi_coef |> Ezfio.set_determinants_psi_coef
end; 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 |> Ezfio.flattened_ezfio
|> Array.map Det_coef.of_float |> Array.map Det_coef.of_float
;; ;;
let write_psi_coef ~n_det ~n_states c = let write_psi_coef ~n_det ~n_states c =
let n_det = Det_number.to_int n_det let n_det = Det_number.to_int n_det
and c = Array.to_list c and c =
|> List.map Det_coef.to_float Array.map Det_coef.to_float c
|> Array.to_list
and n_states = and n_states =
States_number.to_int n_states States_number.to_int n_states
in in
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c let r =
|> Ezfio.set_determinants_psi_coef 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 () let n_int = read_n_int ()
and n_alpha = Ezfio.get_electrons_elec_alpha_num () and n_alpha = Ezfio.get_electrons_elec_alpha_num ()
|> Elec_alpha_number.of_int |> Elec_alpha_number.of_int
@ -232,19 +256,26 @@ end = struct
|> Ezfio.set_determinants_psi_det ; |> Ezfio.set_determinants_psi_det ;
end ; end ;
let n_int = N_int_number.to_int n_int in 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 let dim = psi_det_array.Ezfio.dim
and data = Ezfio.flattened_ezfio psi_det_array and data = Ezfio.flattened_ezfio psi_det_array
in in
assert (n_int = dim.(0)); assert (n_int = dim.(0));
assert (dim.(1) = 2); assert (dim.(1) = 2);
assert (dim.(2) = (Det_number.to_int (read_n_det ()))); if read_only then
List.init dim.(2) (fun i -> 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) ) 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) ~n_int:(N_int_number.of_int n_int)
~alpha:n_alpha ~beta:n_beta ) ~alpha:n_alpha ~beta:n_beta )
|> Array.of_list
;; ;;
let write_psi_det ~n_int ~n_det d = let write_psi_det ~n_int ~n_det d =
@ -252,40 +283,45 @@ end = struct
|> Array.concat |> Array.concat
|> Array.to_list |> Array.to_list
in 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 let r =
|> Ezfio.set_determinants_psi_det 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 if (Ezfio.has_mo_basis_mo_num ()) then
try
Some
{ n_int = read_n_int () ; { n_int = read_n_int () ;
bit_kind = read_bit_kind () ; bit_kind = read_bit_kind () ;
n_det = read_n_det () ; n_det = read_n_det () ;
n_det_qp_edit = read_n_det_qp_edit () ;
expected_s2 = read_expected_s2 () ; expected_s2 = read_expected_s2 () ;
psi_coef = read_psi_coef () ; psi_coef = read_psi_coef ~read_only () ;
psi_det = read_psi_det () ; psi_det = read_psi_det ~read_only () ;
n_states = read_n_states () ; n_states = read_n_states () ;
state_average_weight = read_state_average_weight () ; state_average_weight = read_state_average_weight () ;
} }
with _ -> None
else else
failwith "No molecular orbitals, so no determinants" (* 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
None None
;; ;;
let write { n_int ; let write { n_int ;
bit_kind ; bit_kind ;
n_det ; n_det ;
n_det_qp_edit ;
expected_s2 ; expected_s2 ;
psi_coef ; psi_coef ;
psi_det ; psi_det ;
@ -297,9 +333,13 @@ end = struct
write_n_det n_det; write_n_det n_det;
write_n_states n_states; write_n_states n_states;
write_expected_s2 expected_s2; write_expected_s2 expected_s2;
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ; if n_det <= n_det_qp_edit then
write_psi_det ~n_int:n_int ~n_det:n_det psi_det; begin
write_state_average_weight state_average_weight; 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 |> States_number.to_int
and ndet = and ndet =
Det_number.to_int b.n_det Det_number.to_int b.n_det
and ndet_qp_edit =
Det_number.to_int b.n_det_qp_edit
in in
let coefs_string i = let coefs_string i =
Array.init nstates (fun j -> Array.init nstates (fun j ->
let ishift = let ishift =
j*ndet j*ndet_qp_edit
in in
if (ishift < Array.length b.psi_coef) then if (ishift < Array.length b.psi_coef) then
b.psi_coef.(i+ishift) b.psi_coef.(i+ishift)
@ -331,7 +373,7 @@ end = struct
) )
|> Array.to_list |> String.concat "\t" |> Array.to_list |> String.concat "\t"
in in
Array.init ndet (fun i -> Array.init ndet_qp_edit (fun i ->
Printf.sprintf " %s\n%s\n" Printf.sprintf " %s\n%s\n"
(coefs_string i) (coefs_string i)
(Determinant.to_string ~mo_num:mo_num b.psi_det.(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.expected_s2 |> Positive_float.to_string)
(b.n_det |> Det_number.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 det_text
|> Rst_string.of_string |> Rst_string.of_string
;; ;;
@ -387,10 +429,10 @@ psi_det = %s
(b.n_states |> States_number.to_string) (b.n_states |> States_number.to_string)
(b.expected_s2 |> Positive_float.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.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 ", ") |> String.concat ", ")
(b.psi_det |> Array.to_list |> List.map (Determinant.to_string (b.psi_det |> Array.map (Determinant.to_string ~mo_num) |> Array.to_list
~mo_num) |> String.concat "\n\n") |> String.concat "\n\n")
;; ;;
let of_rst r = let of_rst r =
@ -472,6 +514,7 @@ psi_det = %s
(* Handle determinants *) (* Handle determinants *)
let psi_det = 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 () let n_alpha = Ezfio.get_electrons_elec_alpha_num ()
|> Elec_alpha_number.of_int |> Elec_alpha_number.of_int
and n_beta = Ezfio.get_electrons_elec_beta_num () and n_beta = Ezfio.get_electrons_elec_beta_num ()
@ -483,8 +526,8 @@ psi_det = %s
begin begin
let newdet = let newdet =
(Bitlist.of_string ~zero:'-' ~one:'+' alpha , (Bitlist.of_string ~zero:'-' ~one:'+' alpha ,
Bitlist.of_string ~zero:'-' ~one:'+' beta) Bitlist.of_string ~zero:'-' ~one:'+' beta)
|> Determinant.of_bitlist_couple ~alpha:n_alpha ~beta:n_beta |> Determinant.of_bitlist_couple ~n_int ~alpha:n_alpha ~beta:n_beta
|> Determinant.sexp_of_t |> Determinant.sexp_of_t
|> Sexplib.Sexp.to_string |> Sexplib.Sexp.to_string
in in
@ -492,9 +535,6 @@ psi_det = %s
end end
| _::tail -> read_dets accu tail | _::tail -> read_dets accu tail
in in
let dets =
List.map String_ext.rev dets
in
let a = let a =
read_dets [] dets read_dets [] dets
|> String.concat "" |> String.concat ""
@ -510,9 +550,11 @@ psi_det = %s
Printf.sprintf "(n_int %d)" (N_int_number.get_max ()) Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
and n_states = and n_states =
Printf.sprintf "(n_states %d)" (States_number.to_int @@ read_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 in
let s = 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 in
@ -527,7 +569,9 @@ psi_det = %s
Det_number.to_int n_det_new Det_number.to_int n_det_new
in in
let det = let det =
read () match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in in
let n_det_old, n_states = let n_det_old, n_states =
Det_number.to_int det.n_det, Det_number.to_int det.n_det,
@ -558,7 +602,9 @@ psi_det = %s
let extract_state istate = let extract_state istate =
Printf.printf "Extracting state %d\n" (States_number.to_int istate); Printf.printf "Extracting state %d\n" (States_number.to_int istate);
let det = let det =
read () match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in in
let n_det, n_states = let n_det, n_states =
Det_number.to_int det.n_det, Det_number.to_int det.n_det,
@ -588,7 +634,9 @@ psi_det = %s
let extract_states range = let extract_states range =
Printf.printf "Extracting states %s\n" (Range.to_string range); Printf.printf "Extracting states %s\n" (Range.to_string range);
let det = let det =
read () match read () with
| Some x -> x
| None -> failwith "No determinants in file"
in in
let n_det, n_states = let n_det, n_states =
Det_number.to_int det.n_det, Det_number.to_int det.n_det,
@ -614,7 +662,8 @@ psi_det = %s
j*n_det j*n_det
in in
for i=0 to (n_det-1) do 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 done
end; end;
state_shift := !state_shift + n_det state_shift := !state_shift + n_det

View File

@ -65,8 +65,15 @@ end = struct
let read_mo_num () = let read_mo_num () =
Ezfio.get_mo_basis_mo_num () let elec_alpha_num =
|> MO_number.of_int 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 () = let read_mo_class () =

View File

@ -43,7 +43,7 @@ $(QP_ROOT)/data/executables: remake_executables element_create_db.byte Qptypes.m
$(QP_ROOT)/ocaml/element_create_db.byte $(QP_ROOT)/ocaml/element_create_db.byte
external_libs: external_libs:
opam install cryptokit core opam install cryptokit sexplib
qpackage.odocl: $(MLIFILES) qpackage.odocl: $(MLIFILES)
ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl

View File

@ -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 true: thread
false: profile false: profile
<*byte> : linkdep(c_bindings.o), custom <*byte> : linkdep(c_bindings.o), custom

View File

@ -44,9 +44,12 @@ let psi_det () =
let psi_det = let psi_det =
Input.Determinants_by_hand.read () Input.Determinants_by_hand.read ()
in in
Input.Determinants_by_hand.to_rst psi_det match psi_det with
|> Rst_string.to_string | Some psi_det ->
|> print_endline Input.Determinants_by_hand.to_rst psi_det
|> Rst_string.to_string
|> print_endline
| None -> ()

View File

@ -112,9 +112,9 @@ let set ~core ~inact ~act ~virt ~del =
and av = Excitation.create_single act virt and av = Excitation.create_single act virt
in in
let single_excitations = [ ia ; aa ; av ] let single_excitations = [ ia ; aa ; av ]
|> List.map (fun x -> |> List.map (fun z ->
let open Excitation in let open Excitation in
match x with match z with
| Single (x,y) -> | Single (x,y) ->
( MO_class.to_bitlist n_int (Hole.to_mo_class x), ( MO_class.to_bitlist n_int (Hole.to_mo_class x),
MO_class.to_bitlist n_int (Particle.to_mo_class y) ) 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 match aa with
| Double _ -> assert false | Double _ -> assert false
| Single (x,y) -> | Single (x,y) ->
( 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) ) ( MO_class.to_bitlist n_int ( Hole.to_mo_class x) ) @
|> Bitlist.to_int64_list Bitlist.to_int64_list
( MO_class.to_bitlist n_int (Particle.to_mo_class y) )
in in
Ezfio.set_bitmasks_n_mask_cas 1; 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 Ezfio.ezfio_array_of_list ~rank:3 ~dim:([| (N_int_number.to_int n_int) ; 2; 1|]) ~data:result

470
ocaml/qp_tunnel.ml Normal file
View File

@ -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 <host> &
%!"
(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"

View File

@ -58,7 +58,7 @@ let input_data = "
* Det_number_max : int * Det_number_max : int
assert (x > 0) ; assert (x > 0) ;
if (x > 50_00_000_000) then if (x > 50_000_000_000) then
warning \"More than 50 billion determinants\"; warning \"More than 50 billion determinants\";
* States_number : int * States_number : int

View File

@ -79,7 +79,7 @@ let get s =
| Ao_basis -> | Ao_basis ->
f Ao_basis.(read, to_rst) f Ao_basis.(read, to_rst)
| Determinants_by_hand -> | Determinants_by_hand ->
f Determinants_by_hand.(read_maybe, to_rst) f Determinants_by_hand.(read ~full:false, to_rst)
{section_to_rst} {section_to_rst}
end end
with with

View File

@ -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: 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_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis
* `ao_nucl_elec_integral` which are the nuclear-elctron operator integrals on the |AO| basis (see :file:`pot_ao_ints.irp.f`) * `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 (see :file:`ao_mono_ints.irp.f`) * `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`.

View File

@ -2,6 +2,29 @@ use bitmasks
integer function number_of_holes(key_in) integer function number_of_holes(key_in)
BEGIN_DOC BEGIN_DOC
! Function that returns the number of holes in the inact space ! 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 END_DOC
implicit none implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2) 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(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,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)) ) + 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 else
do i = 1, N_int do i = 1, N_int
number_of_holes = number_of_holes & number_of_holes = number_of_holes &
@ -108,7 +115,7 @@ integer function number_of_holes(key_in)
) & ) &
) & ) &
), reunion_of_core_inact_bitmask(i,1)) ) & ), 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 enddo
endif endif
end end

View File

@ -141,8 +141,8 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
b%cur=0 b%cur=0
! ! Try to adjust n_tasks around nproc/2 seconds per job ! ! 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 = min(2*n_tasks,int( dble(n_tasks * nproc/2) / (time1 - time0 + 1.d0)))
n_tasks = 1 ! n_tasks = 1
end do end do
integer, external :: disconnect_from_taskserver integer, external :: disconnect_from_taskserver

View File

@ -85,7 +85,7 @@ subroutine run_selection_slave(thread,iproc,energy)
if(ctask > 0) then if(ctask > 0) then
call sort_selection_buffer(buf) call sort_selection_buffer(buf)
! call merge_selection_buffers(buf,buf2) ! 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) call push_selection_results(zmq_socket_push, pt2, variance, norm, buf, task_id(1), ctask)
! buf%mini = buf2%mini ! buf%mini = buf2%mini
pt2(:) = 0d0 pt2(:) = 0d0

View File

@ -66,6 +66,12 @@ interface: ezfio
doc: Number of determinants in the current wave function doc: Number of determinants in the current wave function
type: integer type: integer
[n_det_qp_edit]
interface: ezfio
doc: Number of determinants to print in qp_edit
type: integer
interface: ezfio
[psi_coef] [psi_coef]
interface: ezfio interface: ezfio
doc: Coefficients of the wave function doc: Coefficients of the wave function
@ -78,6 +84,18 @@ doc: Determinants of the variational space
type: integer*8 type: integer*8
size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det) 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] [expected_s2]
interface: ezfio interface: ezfio
doc: Expected value of |S^2| doc: Expected value of |S^2|

View File

@ -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(bit_kind), intent(inout) :: key_in(N_int,2)
integer, intent(out) :: i_ok integer, intent(out) :: i_ok
integer :: k,j,i integer :: k,j,i
integer(bit_kind) :: mask
use bitmasks use bitmasks
ASSERT (i_hole > 0 ) ASSERT (i_hole > 0 )
ASSERT (i_particle <= mo_num) ASSERT (i_particle <= mo_num)
@ -19,31 +20,66 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
! hole ! hole
k = shiftr(i_hole-1,bit_kind_shift)+1 k = shiftr(i_hole-1,bit_kind_shift)+1
j = i_hole-shiftl(k-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 ! 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) key_in(k,ispin) = ibclr(key_in(k,ispin),j)
else else
i_ok= -1 i_ok= -1
return
end if end if
! particle ! particle
k = shiftr(i_particle-1,bit_kind_shift)+1 k = shiftr(i_particle-1,bit_kind_shift)+1
j = i_particle-shiftl(k-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 ! integer :: n_elec_tmp
n_elec_tmp = 0 ! n_elec_tmp = 0
do i = 1, N_int ! do i = 1, N_int
n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2)) ! n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
enddo ! enddo
if(n_elec_tmp .ne. elec_num)then ! if(n_elec_tmp .ne. elec_num)then
!print*, n_elec_tmp,elec_num ! print*, n_elec_tmp,elec_num
!call debug_det(key_in,N_int) ! call debug_det(key_in,N_int)
i_ok = -1 ! stop -1
endif ! endif
end 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) logical function is_spin_flip_possible(key_in,i_flip,ispin)
implicit none implicit none
BEGIN_DOC BEGIN_DOC

View File

@ -269,7 +269,6 @@ subroutine save_natural_mos
END_DOC END_DOC
call set_natural_mos call set_natural_mos
call save_mos call save_mos
end end

View File

@ -44,6 +44,16 @@ BEGIN_PROVIDER [ integer, N_det ]
ASSERT (N_det > 0) ASSERT (N_det > 0)
END_PROVIDER 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] BEGIN_PROVIDER [integer, max_degree_exc]
implicit none implicit none
integer :: i,degree integer :: i,degree
@ -476,7 +486,7 @@ subroutine save_wavefunction
endif endif
if (mpi_master) then if (mpi_master) then
call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
endif endif
end end
@ -504,12 +514,16 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
integer*8, allocatable :: psi_det_save(:,:,:) integer*8, allocatable :: psi_det_save(:,:,:)
double precision, allocatable :: psi_coef_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 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_N_int(N_int)
call ezfio_set_determinants_bit_kind(bit_kind) call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_N_det(ndet) 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_n_states(nstates)
call ezfio_set_determinants_mo_label(mo_label) call ezfio_set_determinants_mo_label(mo_label)
@ -522,10 +536,10 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef)
enddo enddo
enddo enddo
call ezfio_set_determinants_psi_det(psi_det_save) call ezfio_set_determinants_psi_det(psi_det_save)
call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
deallocate (psi_det_save) deallocate (psi_det_save)
allocate (psi_coef_save(ndet,nstates)) allocate (psi_coef_save(ndet,nstates))
double precision :: accu_norm
do k=1,nstates do k=1,nstates
do i=1,ndet do i=1,ndet
psi_coef_save(i,k) = psicoef(i,k) 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) call ezfio_set_determinants_psi_coef(psi_coef_save)
deallocate (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') call write_int(6,ndet,'Saved determinants')
endif endif
end end
@ -559,54 +585,80 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde
integer :: N_int2 integer :: N_int2
equivalence (det_8, det_bk) equivalence (det_8, det_bk)
integer :: i,k integer :: i,j,k, ndet_qp_edit
call ezfio_set_determinants_N_int(N_int) if (mpi_master) then
call ezfio_set_determinants_bit_kind(bit_kind) ndet_qp_edit = min(ndetsave,N_det_qp_edit)
call ezfio_set_determinants_N_det(ndetsave) call ezfio_set_determinants_N_int(N_int)
call ezfio_set_determinants_n_states(nstates) call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_mo_label(mo_label) 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 N_int2 = (N_int*bit_kind)/8
allocate (psi_det_save(N_int2,2,ndetsave)) 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
do i=1,ndetsave do i=1,ndetsave
accu_norm(k) = accu_norm(k) + psicoef(index_det_save(i),k) * psicoef(index_det_save(i),k) do k=1,N_int
psi_coef_save(i,k) = psicoef(index_det_save(i),k) 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
enddo call ezfio_set_determinants_psi_det(psi_det_save)
do k = 1, nstates call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) deallocate (psi_det_save)
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) allocate (psi_coef_save(ndetsave,nstates))
call write_int(6,ndet,'Saved determinants') double precision :: accu_norm(nstates)
deallocate (psi_coef_save) 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 end

View File

@ -53,7 +53,17 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
use bitmasks use bitmasks
implicit none implicit none
BEGIN_DOC 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 END_DOC
integer ,intent(in) :: Nint integer ,intent(in) :: Nint
integer ,intent(in) :: n_alpha ! Number of alpha electrons integer ,intent(in) :: n_alpha ! Number of alpha electrons

View File

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

View File

@ -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,(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 implicit none
BEGIN_DOC 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 ! 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 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) 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 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 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)] BEGIN_PROVIDER[double precision, mos_lapl_in_r_array,(mo_num,n_points_final_grid,3)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC

View File

@ -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' include 'constants.include.F'
! Input variables ! Input variables
double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2 double precision, intent(in) :: rho,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_2
! Local variables ! Local variables
double precision :: a,b,c,d,c_f,omega,delta 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 :: rho_13,rho_inv_13,rho_83,rho_113,rho_inv_113,denom
double precision :: thr,huge_num,rho_inv 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 :: tmp1,tmp2,tmp3,tmp4
double precision :: big1,big2,big3 double precision :: big1,big2,big3
! Output variables
! Constants of the LYP correlation functional ! 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 c = 0.2533d0
d = 0.349d0 d = 0.349d0
thr = 1d-10 ec_lyp_88 = 0.d0
thr = 1d-15
huge_num = 1.d0/thr 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 if(rho.lt.0.d0)then
print*,'pb !! rho.lt.0.d0' print*,'pb !! rho.lt.0.d0'
stop stop
endif 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 if(dabs(rho_13) < thr) then
rho_inv_13 = huge_num 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 if (dabs(rho_113) < thr) then
rho_inv_113 = huge_num rho_inv_113 = huge_num
else else
rho_inv_113 = 1d0/rho_113 rho_inv_113 = 1.d0/rho_113
endif endif
if (dabs(rho) < thr) then if (dabs(rho) < thr) then
rho_inv = huge_num rho_inv = huge_num
else else
rho_inv = 1d0/rho rho_inv = 1.d0/rho
endif endif
! Useful quantities to predefine ! 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) denom = 1d0/(1d0 + d*rho_inv_13)
omega = rho_inv_113*exp(-c*rho_inv_13)*denom omega = rho_inv_113*exp(-c*rho_inv_13)*denom
delta = c*rho_inv_13 + d*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_2 = rho *rho
rho_a_2 = rho_a*rho_a rho_a_2 = rho_a*rho_a
rho_b_2 = rho_b*rho_b rho_b_2 = rho_b*rho_b
cst_2_113 = 2d0**(11d0/3d0) cst_2_113 = 2.d0**(11.d0/3.d0)
cst_8_3 = 8d0/3d0 cst_8_3 = 8.d0/3.d0
! first term in the equation (2) of Preuss CPL, 1989 ! 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) 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) 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) 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) 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) tmp3 = grad_rho_a_2*(2d0/3d0*rho_2 - rho_b_2)
big3 = tmp1 + tmp2 + tmp3 big3 = tmp1 + tmp2 + tmp3
ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3 ec_lyp_88 = -a*big1 -a*b*omega*big2 -a*b*omega*big3
end end

View File

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

View File

@ -36,24 +36,24 @@ function run_stoch() {
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file f2.ezfio qp set_file f2.ezfio
qp set_frozen_core qp set_frozen_core
run_stoch -199.30496 1.e-4 run_stoch -199.30486 1.e-4
} }
@test "NH3" { # 10.6657s @test "NH3" { # 10.6657s
qp set_file nh3.ezfio qp set_file nh3.ezfio
qp set_mo_class --core="[1-4]" --act="[5-72]" 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 @test "DHNO" { # 11.4721s
qp set_file dhno.ezfio qp set_file dhno.ezfio
qp set_mo_class --core="[1-7]" --act="[8-64]" 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 @test "HCO" { # 12.2868s
qp set_file hco.ezfio qp set_file hco.ezfio
run -113.297494345682 2.e-05 run -113.297494345682 1.e-4
} }
@test "H2O2" { # 12.9214s @test "H2O2" { # 12.9214s
@ -65,82 +65,82 @@ function run_stoch() {
@test "HBO" { # 13.3144s @test "HBO" { # 13.3144s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file hbo.ezfio qp set_file hbo.ezfio
run -100.214185815312 1.e-5 run -100.212829869715 1.e-4
} }
@test "H2O" { # 11.3727s @test "H2O" { # 11.3727s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file h2o.ezfio qp set_file h2o.ezfio
run -76.2359268957699 2.e-5 run -76.2359268957699 1.e-4
} }
@test "ClO" { # 13.3755s @test "ClO" { # 13.3755s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file clo.ezfio qp set_file clo.ezfio
run -534.546005867797 5.e-5 run -534.545881614967 1.e-4
} }
@test "SO" { # 13.4952s @test "SO" { # 13.4952s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file so.ezfio qp set_file so.ezfio
run -26.0124797722154 1.e-5 run -26.0158153138924 1.e-4
} }
@test "H2S" { # 13.6745s @test "H2S" { # 13.6745s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file h2s.ezfio qp set_file h2s.ezfio
run -398.859480581924 1.e-5 run -398.859168655255 1.e-4
} }
@test "OH" { # 13.865s @test "OH" { # 13.865s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file oh.ezfio qp set_file oh.ezfio
run -75.6119887538831 1.e-05 run -75.6120779012574 1.e-4
} }
@test "SiH2_3B1" { # 13.938ss @test "SiH2_3B1" { # 13.938ss
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file sih2_3b1.ezfio qp set_file sih2_3b1.ezfio
run -290.017539006762 1.e-5 run -290.017539006762 1.e-4
} }
@test "H3COH" { # 14.7299s @test "H3COH" { # 14.7299s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file h3coh.ezfio qp set_file h3coh.ezfio
run -115.205054063687 1.e-5 run -115.205941463667 1.e-4
} }
@test "SiH3" { # 15.99s @test "SiH3" { # 15.99s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file sih3.ezfio qp set_file sih3.ezfio
run -5.57269434557089 2.e-05 run -5.57241217753818 1.e-4
} }
@test "CH4" { # 16.1612s @test "CH4" { # 16.1612s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file ch4.ezfio qp set_file ch4.ezfio
qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" 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 @test "ClF" { # 16.8864s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file clf.ezfio qp set_file clf.ezfio
run -559.170406471496 1.e-5 run -559.170272077166 1.e-4
} }
@test "SO2" { # 17.5645s @test "SO2" { # 17.5645s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file so2.ezfio qp set_file so2.ezfio
qp set_mo_class --core="[1-8]" --act="[9-87]" 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 @test "C2H2" { # 17.6827s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file c2h2.ezfio qp set_file c2h2.ezfio
qp set_mo_class --act="[1-30]" --del="[31-36]" 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 @test "N2" { # 18.0198s
@ -154,14 +154,14 @@ function run_stoch() {
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file n2h4.ezfio qp set_file n2h4.ezfio
qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-48]" 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 @test "CO2" { # 21.1748s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file co2.ezfio qp set_file co2.ezfio
qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]" 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 [[ -n $TRAVIS ]] && skip
qp set_file cu_nh3_4_2plus.ezfio qp set_file cu_nh3_4_2plus.ezfio
qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]" 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 @test "HCN" { # 20.3273s
[[ -n $TRAVIS ]] && skip [[ -n $TRAVIS ]] && skip
qp set_file hcn.ezfio qp set_file hcn.ezfio
qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]" 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
} }

View File

@ -21,7 +21,7 @@ type: double precision
size: (nuclei.nucl_num,3) size: (nuclei.nucl_num,3)
interface: ezfio interface: ezfio
[disk_access_nuclear_repulsion] [io_nuclear_repulsion]
doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ] doc: Read/Write Nuclear Repulsion from/to disk [ Write | Read | None ]
type: Disk_access type: Disk_access
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml

View File

@ -143,7 +143,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
END_DOC END_DOC
PROVIDE mpi_master nucl_coord nucl_charge nucl_num 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 logical :: has
if (mpi_master) then if (mpi_master) then
@ -194,7 +194,7 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
call write_time(6) call write_time(6)
call write_double(6,nuclear_repulsion,'Nuclear repulsion energy') 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 if (mpi_master) then
call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion) call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)
endif endif

View File

@ -24,6 +24,7 @@ subroutine routine
implicit none implicit none
integer :: i,k integer :: i,k
integer :: degree integer :: degree
call print_energy_components
do i = 1, N_det do i = 1, N_det
print *, 'Determinant ', i print *, 'Determinant ', i
call debug_det(psi_det(1,1,i),N_int) call debug_det(psi_det(1,1,i),N_int)