LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-24 13:23:41 +01:00
Code Issues Releases Wiki Activity

merge with garniroy

Browse Source
This commit is contained in:
Yann Garniron 2018-09-04 20:07:15 +02:00
commit ba0094f5f8
107 changed files with 6419 additions and 1766 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
README.md
View File

@ -163,3 +163,6 @@ It's caused when we call the DGEMM routine of LAPACK.
Set `ulimit -s unlimited`, before runing `qp_run`. It seems to fix the problem.
#### Error: f77zmq not symbol found
In the Makefile of the f77zmq, you should NOT use ar but `libtool -static`

2
config/ifort.cfg
View File

@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags
#
[OPT]
FC : -traceback
FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g
FCFLAGS : -xAVX -O2 -ip -ftz -g
# Profiling flags
#################

22
configure vendored
View File

@ -76,6 +76,11 @@ from collections import namedtuple
Info = namedtuple("Info", ["url", "description", "default_path"])
path_github = {"head": "http://github.com", "tail": "archive/master.tar.gz"}
def path_gitlab(user,project,branch):
"""For example,
https://gitlab.com/scemama/resultsFile/-/archive/master/resultsFile-master.tar.gz
"""
return "http://gitlab.com/%(user)s/%(project)s/-/archive/%(branch)s/%(project)s-%(branch)s.tar.gz"%locals()
ocaml = Info(
url='http://raw.github.com/ocaml/opam/master/shell/opam_installer.sh',
@ -113,7 +118,7 @@ docopt = Info(
default_path=join(QP_ROOT_INSTALL, "docopt"))
resultsFile = Info(
url='{head}/LCPQ/resultsFile/{tail}'.format(**path_github),
url=path_gitlab("scemama","resultsFile","master"),
description=' resultsFile',
default_path=join(QP_ROOT_INSTALL, "resultsFile"))
@ -128,17 +133,17 @@ emsl = Info(
default_path=join(QP_ROOT_INSTALL, "emsl"))
ezfio = Info(
url='{head}/LCPQ/EZFIO/{tail}'.format(**path_github),
url=path_gitlab("scemama","EZFIO","master"),
description=' EZFIO',
default_path=join(QP_ROOT_INSTALL, "EZFIO"))
zeromq = Info(
url='https://github.com/zeromq/zeromq4-1/releases/download/v4.1.5/zeromq-4.1.5.tar.gz',
url='https://github.com/zeromq/libzmq/archive/v4.2.5.tar.gz',
description=' ZeroMQ',
default_path=join(QP_ROOT_LIB, "libzmq.a"))
gmp= Info(
url='https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2',
url='ftp://ftp.gnu.org/gnu/gmp/gmp-6.1.2.tar.bz2',
description=' The GNU Multiple Precision Arithmetic Library',
default_path=join(QP_ROOT_LIB, "libgmp.a"))
@ -487,9 +492,12 @@ def create_ninja_and_rc(l_installed):
'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
'function qp_prepend_export () {',
' #Prepend path $2:${!1}. Add the semicolon only if ${!1} is defined',
' eval "value_1=\"\${$1}\""',
' echo ${value_1:+${value_1}:}${2}',
'eval "value_1="\${$1}""',
'if [[ -z $value_1 ]] ; then',
' echo "${2}:"',
'else',
' echo "${2}:${value_1}"',
'fi',
'}',
'export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")',
'export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)',

635
data/basis/6-31+g Normal file
View File

@ -0,0 +1,635 @@
$DATA
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
2 96.7985150 0.0162089
3 22.0911210 0.0773156
4 6.2010703 0.2457860
5 1.9351177 0.4701890
6 0.6367358 0.3454708
S 3
1 2.3249184 -0.0350917
2 0.6324306 -0.1912328
3 0.0790534 1.0839878
P 3
1 2.3249184 0.0089415
2 0.6324306 0.1410095
3 0.0790534 0.9453637
S 1
1 0.0359620 1.0000000
P 1
1 0.0359620 1.0000000
S 1
1 0.0074000 1.0000000
P 1
1 0.0074000 1.0000000
BERYLLIUM
S 6
1 1264.5857000 0.0019448
2 189.9368100 0.0148351
3 43.1590890 0.0720906
4 12.0986630 0.2371542
5 3.8063232 0.4691987
6 1.2728903 0.3565202
S 3
1 3.1964631 -0.1126487
2 0.7478133 -0.2295064
3 0.2199663 1.1869167
P 3
1 3.1964631 0.0559802
2 0.7478133 0.2615506
3 0.2199663 0.7939723
S 1
1 0.0823099 1.0000000
P 1
1 0.0823099 1.0000000
S 1
1 0.0207000 1.0000000
P 1
1 0.0207000 1.0000000
BORON
S 6
1 2068.8823000 0.0018663
2 310.6495700 0.0142515
3 70.6830330 0.0695516
4 19.8610800 0.2325729
5 6.2993048 0.4670787
6 2.1270270 0.3634314
S 3
1 4.7279710 -0.1303938
2 1.1903377 -0.1307889
3 0.3594117 1.1309444
P 3
1 4.7279710 0.0745976
2 1.1903377 0.3078467
3 0.3594117 0.7434568
S 1
1 0.1267512 1.0000000
P 1
1 0.1267512 1.0000000
S 1
1 0.0315000 1.0000000
P 1
1 0.0315000 1.0000000
CARBON
S 6
1 3047.5249000 0.0018347
2 457.3695100 0.0140373
3 103.9486900 0.0688426
4 29.2101550 0.2321844
5 9.2866630 0.4679413
6 3.1639270 0.3623120
S 3
1 7.8682724 -0.1193324
2 1.8812885 -0.1608542
3 0.5442493 1.1434564
P 3
1 7.8682724 0.0689991
2 1.8812885 0.3164240
3 0.5442493 0.7443083
S 1
1 0.1687144 1.0000000
P 1
1 0.1687144 1.0000000
S 1
1 0.0438000 1.0000000
P 1
1 0.0438000 1.0000000
NITROGEN
S 6
1 4173.5110000 0.0018348
2 627.4579000 0.0139950
3 142.9021000 0.0685870
4 40.2343300 0.2322410
5 12.8202100 0.4690700
6 4.3904370 0.3604550
S 3
1 11.6263580 -0.1149610
2 2.7162800 -0.1691180
3 0.7722180 1.1458520
P 3
1 11.6263580 0.0675800
2 2.7162800 0.3239070
3 0.7722180 0.7408950
S 1
1 0.2120313 1.0000000
P 1
1 0.2120313 1.0000000
S 1
1 0.0639000 1.0000000
P 1
1 0.0639000 1.0000000
OXYGEN
S 6
1 5484.6717000 0.0018311
2 825.2349500 0.0139501
3 188.0469600 0.0684451
4 52.9645000 0.2327143
5 16.8975700 0.4701930
6 5.7996353 0.3585209
S 3
1 15.5396160 -0.1107775
2 3.5999336 -0.1480263
3 1.0137618 1.1307670
P 3
1 15.5396160 0.0708743
2 3.5999336 0.3397528
3 1.0137618 0.7271586
S 1
1 0.2700058 1.0000000
P 1
1 0.2700058 1.0000000
S 1
1 0.0845000 1.0000000
P 1
1 0.0845000 1.0000000
FLUORINE
S 6
1 7001.7130900 0.0018196169
2 1051.3660900 0.0139160796
3 239.2856900 0.0684053245
4 67.3974453 0.233185760
5 21.5199573 0.471267439
6 7.40310130 0.356618546
S 3
1 20.8479528 -0.108506975
2 4.80830834 -0.146451658
3 1.34406986 1.128688580
P 3
1 20.8479528 0.0716287243
2 4.80830834 0.3459121030
3 1.34406986 0.7224699570
S 1
1 0.358151393 1.0000000
P 1
1 0.358151393 1.0000000
S 1
1 0.1076000 1.0000000
P 1
1 0.1076000 1.0000000
NEON
S 6
1 8425.8515300 0.0018843481
2 1268.5194000 0.0143368994
3 289.6214140 0.0701096233
4 81.8590040 0.2373732660
5 26.2515079 0.4730071260
6 9.09472051 0.3484012410
S 3
1 26.5321310 -0.107118287
2 6.10175501 -0.146163821
3 1.69627153 1.127773500
P 3
1 26.5321310 0.0719095885
2 6.10175501 0.3495133720
3 1.69627153 0.7199405120
S 1
1 0.44581870 1.0000000
P 1
1 0.44581870 1.0000000
S 1
1 0.1300000 1.0000000
P 1
1 0.1300000 1.0000000
SODIUM
S 6
1 9993.2000000 0.0019377
2 1499.8900000 0.0148070
3 341.9510000 0.0727060
4 94.6797000 0.2526290
5 29.7345000 0.4932420
6 10.0063000 0.3131690
S 6
1 150.9630000 -0.0035421
2 35.5878000 -0.0439590
3 11.1683000 -0.1097521
4 3.9020100 0.1873980
5 1.3817700 0.6466990
6 0.4663820 0.3060580
P 6
1 150.9630000 0.0050017
2 35.5878000 0.0355110
3 11.1683000 0.1428250
4 3.9020100 0.3386200
5 1.3817700 0.4515790
6 0.4663820 0.2732710
S 3
1 0.4979660 -0.2485030
2 0.0843530 -0.1317040
3 0.0666350 1.2335200
P 3
1 0.4979660 -0.0230230
2 0.0843530 0.9503590
3 0.0666350 0.0598580
S 1
1 0.0259544 1.0000000
P 1
1 0.0259544 1.0000000
S 1
1 0.0076000 1.0000000
P 1
1 0.0076000 1.0000000
MAGNESIUM
S 6
1 11722.8000000 0.0019778
2 1759.9300000 0.0151140
3 400.8460000 0.0739110
4 112.8070000 0.2491910
5 35.9997000 0.4879280
6 12.1828000 0.3196620
S 6
1 189.1800000 -0.0032372
2 45.2119000 -0.0410080
3 14.3563000 -0.1126000
P 6
1 189.1800000 0.0049281
2 45.2119000 0.0349890
3 14.3563000 0.1407250
4 5.1388600 0.1486330 0.3336420
5 1.9065200 0.6164970 0.4449400
6 0.7058870 0.3648290 0.2692540
S 3
1 0.9293400 -0.2122900
2 0.2690350 -0.1079850
3 0.1173790 1.1758400
P 3
1 0.9293400 -0.0224190
2 0.2690350 0.1922700
3 0.1173790 0.8461810
S 1
1 0.0421061 1.0000000
P 1
1 0.0421061 1.0000000
S 1
1 0.0146000 1.0000000
P 1
1 0.0146000 1.0000000
ALUMINUM
S 6
1 13983.1000000 0.00194267
2 2098.7500000 0.0148599
3 477.7050000 0.0728494
4 134.3600000 0.2468300
5 42.8709000 0.4872580
6 14.5189000 0.3234960
S 6
1 239.6680000 -0.00292619
2 57.4419000 -0.0374080
3 18.2859000 -0.1144870
4 6.5991400 0.1156350
5 2.4904900 0.6125950
6 0.9445400 0.3937990
P 6
1 239.6680000 0.00460285
2 57.4419000 0.0331990
3 18.2859000 0.1362820
4 6.5991400 0.3304760
5 2.4904900 0.4491460
6 0.9445400 0.2657040
S 3
1 1.2779000 -0.2276060
2 0.3975900 0.00144583
3 0.1600950 1.0927900
P 3
1 1.2779000 -0.0175130
2 0.3975900 0.2445330
3 0.1600950 0.8049340
S 1
1 0.0556577 1.0000000
P 1
1 0.0556577 1.0000000
S 1
1 0.0318000 1.0000000
P 1
1 0.0318000 1.0000000
SILICON
S 6
1 16115.9000000 0.00195948
2 2425.5800000 0.01492880
3 553.8670000 0.07284780
4 156.3400000 0.24613000
5 50.0683000 0.48591400
6 17.0178000 0.32500200
S 6
1 292.7180000 -0.00278094
2 69.8731000 -0.03571460
3 22.3363000 -0.11498500
4 8.1503900 0.09356340
5 3.1345800 0.60301700
6 1.2254300 0.41895900
P 6
1 292.7180000 0.00443826
2 69.8731000 0.03266790
3 22.3363000 0.13472100
4 8.1503900 0.32867800
5 3.1345800 0.44964000
6 1.2254300 0.26137200
S 3
1 1.7273800 -0.24463000
2 0.5729220 0.00431572
3 0.2221920 1.09818000
P 3
1 1.7273800 -0.01779510
2 0.5729220 0.25353900
3 0.2221920 0.80066900
S 1
1 0.0778369 1.00000000
P 1
1 0.0778369 1.00000000
S 1
1 0.0331000 1.0000000
P 1
1 0.0331000 1.0000000
PHOSPHOROUS
S 6
1 19413.3000000 0.0018516
2 2909.4200000 0.0142062
3 661.3640000 0.0699995
4 185.7590000 0.2400790
5 59.1943000 0.4847620
6 20.0310000 0.3352000
S 6
1 339.4780000 -0.00278217
2 81.0101000 -0.0360499
3 25.8780000 -0.1166310
4 9.4522100 0.0968328
5 3.6656600 0.6144180
6 1.4674600 0.4037980
P 6
1 339.4780000 0.00456462
2 81.0101000 0.03369360
3 25.8780000 0.13975500
4 9.4522100 0.33936200
5 3.6656600 0.45092100
6 1.4674600 0.23858600
S 3
1 2.1562300 -0.2529230
2 0.7489970 0.0328517
3 0.2831450 1.0812500
P 3
1 2.1562300 -0.01776530
2 0.7489970 0.27405800
3 0.2831450 0.78542100
S 1
1 0.0998317 1.0000000
P 1
1 0.0998317 1.0000000
S 1
1 0.0348000 1.0000000
P 1
1 0.0348000 1.0000000
SULFUR
S 6
1 21917.1000000 0.0018690
2 3301.4900000 0.0142300
3 754.1460000 0.0696960
4 212.7110000 0.2384870
5 67.9896000 0.4833070
6 23.0515000 0.3380740
S 6
1 423.7350000 -0.0023767
2 100.7100000 -0.0316930
3 32.1599000 -0.1133170
4 11.8079000 0.0560900
5 4.6311000 0.5922550
6 1.8702500 0.4550060
P 6
1 423.7350000 0.0040610
2 100.7100000 0.0306810
3 32.1599000 0.1304520
4 11.8079000 0.3272050
5 4.6311000 0.4528510
6 1.8702500 0.2560420
S 3
1 2.6158400 -0.2503740
2 0.9221670 0.0669570
3 0.3412870 1.0545100
P 3
1 2.6158400 -0.0145110
2 0.9221670 0.3102630
3 0.3412870 0.7544830
S 1
1 0.1171670 1.0000000
P 1
1 0.1171670 1.0000000
S 1
1 0.0405000 1.0000000
P 1
1 0.0405000 1.0000000
CHLORINE
S 6
1 25180.1000000 0.0018330
2 3780.3500000 0.0140340
3 860.4740000 0.0690970
4 242.1450000 0.2374520
5 77.3349000 0.4830340
6 26.2470000 0.3398560
S 6
1 491.7650000 -0.0022974
2 116.9840000 -0.0307140
3 37.4153000 -0.1125280
4 13.7834000 0.0450160
5 5.4521500 0.5893530
6 2.2258800 0.4652060
P 6
1 491.7650000 0.0039894
2 116.9840000 0.0303180
3 37.4153000 0.1298800
4 13.7834000 0.3279510
5 5.4521500 0.4535270
6 2.2258800 0.2521540
S 3
1 3.1864900 -0.2518300
2 1.1442700 0.0615890
3 0.4203770 1.0601800
P 3
1 3.1864900 -0.0142990
2 1.1442700 0.3235720
3 0.4203770 0.7435070
S 1
1 0.1426570 1.0000000
P 1
1 0.1426570 1.0000000
S 1
1 0.0483000 1.0000000
P 1
1 0.0483000 1.0000000
ARGON
S 6
1 28348.3000000 0.00182526
2 4257.6200000 0.01396860
3 969.8570000 0.06870730
4 273.2630000 0.23620400
5 87.3695000 0.48221400
6 29.6867000 0.34204300
S 6
1 575.8910000 -0.00215972
2 136.8160000 -0.02907750
3 43.8098000 -0.11082700
4 16.2094000 0.02769990
5 6.4608400 0.57761300
6 2.6511400 0.48868800
P 6
1 575.8910000 0.00380665
2 136.8160000 0.02923050
3 43.8098000 0.12646700
4 16.2094000 0.32351000
5 6.4608400 0.45489600
6 2.6511400 0.25663000
S 3
1 3.8602800 -0.2555920
2 1.4137300 0.0378066
3 0.5166460 1.0805600
P 3
1 3.8602800 -0.01591970
2 1.4137300 0.32464600
3 0.5166460 0.74399000
S 1
1 0.1738880 1.0000000
P 1
1 0.1738880 1.0000000
S 1
1 0.0600000 1.0000000
P 1
1 0.0600000 1.0000000
POTASSIUM
S 6
1 31594.4200000 1.828010E-03
2 4744.3300000 1.399403E-02
3 1080.4190000 6.887129E-02
4 304.2338000 2.369760E-01
5 97.2458600 4.829040E-01
6 33.0249500 3.404795E-01
S 6
1 622.7625000 -2.502976E-03
2 147.8839000 -3.315550E-02
3 47.3273500 -1.226387E-01
4 17.5149500 5.353643E-02
5 6.9227220 6.193860E-01
6 2.7682770 4.345878E-01
P 6
1 622.7625000 4.094637E-03
2 147.8839000 3.145199E-02
3 47.3273500 1.351558E-01
4 17.5149500 3.390500E-01
5 6.9227220 4.629455E-01
6 2.7682770 2.242638E-01
S 6
1 11.8480200 1.277689E-02
2 4.0792110 2.098767E-01
3 1.7634810 -3.095274E-03
4 0.7889270 -5.593884E-01
5 0.3503870 -5.134760E-01
6 0.1463440 -6.598035E-02
P 6
1 11.8480200 -1.221377E-02
2 4.0792110 -6.900537E-03
3 1.7634810 2.007466E-01
4 0.7889270 4.281332E-01
5 0.3503870 3.970156E-01
6 0.1463440 1.104718E-01
S 3
1 0.7168010 -5.237772E-02
2 0.2337410 -2.798503E-01
3 0.0386750 1.141547E+00
P 3
1 0.7168010 0.0316430
2 0.2337410 -0.0404616
3 0.0386750 1.0120290
S 1
1 0.0165210 1.000000E+00
P 1
1 0.0165210 1.000000E+00
S 1
1 0.0047000 1.0000000
P 1
1 0.0047000 1.0000000
CALCIUM
S 6
1 35264.8600000 1.813501E-03
2 5295.5030000 1.388493E-02
3 1206.0200000 6.836162E-02
4 339.6839000 2.356188E-01
5 108.6264000 4.820639E-01
6 36.9210300 3.429819E-01
S 6
1 706.3096000 2.448225E-03
2 167.8187000 3.241504E-02
3 53.8255800 1.226219E-01
4 20.0163800 -4.316965E-02
5 7.9702790 -6.126995E-01
6 3.2120590 -4.487540E-01
P 6
1 706.3096000 4.020371E-03
2 167.8187000 3.100601E-02
3 53.8255800 1.337279E-01
4 20.0163800 3.367983E-01
5 7.9702790 4.631281E-01
6 3.2120590 2.257532E-01
S 6
1 14.1951800 1.084500E-02
2 4.8808280 2.088333E-01
3 2.1603900 3.150338E-02
4 0.9878990 -5.526518E-01
5 0.4495170 -5.437997E-01
6 0.1873870 -6.669342E-02
P 6
1 14.1951800 -1.289621E-02
2 4.8808280 -1.025198E-02
3 2.1603900 1.959781E-01
4 0.9878990 4.357933E-01
5 0.4495170 3.996452E-01
6 0.1873870 9.713636E-02
S 3
1 1.0322710 -4.439720E-02
2 0.3811710 -3.284563E-01
3 0.0651310 1.163010E+00
P 3
1 1.0322710 -0.4298621
2 0.3811710 0.006935829
3 0.0651310 0.9705933
S 1
1 0.0260100 1.000000E+00
P 1
1 0.0260100 1.000000E+00
S 1
1 0.0071000 1.0000000
P 1
1 0.0071000 1.0000000
$END

1303
data/basis/ncsu-vdz Normal file
View File

File diff suppressed because it is too large Load Diff

1960
data/basis/ncsu-vtz Normal file
View File

File diff suppressed because it is too large Load Diff

183
data/pseudo/ncsu Normal file
View File

@ -0,0 +1,183 @@
H GEN 0 1
3
-10.851924053 2 21.7769665504
1.0 1 21.2435950826
21.2435950826 3 21.2435950826
1
0.0 2 1.0
C GEN 2 1
3
4.0 1 14.43502
57.74008 3 8.39889
-25.81955 2 7.38188
1
52.13345 2 7.76079
Cl GEN 10 2
3
7.0 1 22.71655173
159.01586213 3 78.57185685
-15.6531065 2 7.47352436
2
6.50888648 2 17.23708573
46.763467 2 4.31148447
2
2.9946477 2 11.38275704
28.0170341 2 3.83218762
Co GEN 10 2
4
17.0 1 24.7400138129
420.580234819 3 23.5426031368
-194.630579018 2 24.0406241364
-2.94301943013 2 10.237411369
2
270.86974114 2 23.0205711168
54.1910212498 2 10.9219568474
2
200.63032558 2 25.3244045243
38.9480947892 2 10.6533915029
Cr GEN 10 2
4
14.0 1 18.2809107439
255.932750414 3 17.0980065531
-132.018263171 2 16.7226727605
-0.773887613451 2 5.02865105891
2
219.481462096 2 16.9007876081
28.079331766 2 7.33662150761
2
139.983968717 2 17.3197451654
19.5483578632 2 6.92409757503
F GEN 2 1
3
7.0 1 11.3954401213
79.7680808491 3 10.49201883
-49.4990068225 2 10.2868054098
1
51.3934743997 2 11.3903478843
Fe GEN 10 2
4
16.0 1 23.2209171361
371.534674178 3 23.5471467972
-181.226034452 2 23.4725634461
-2.3730523614 2 9.85238815041
2
277.500325475 2 22.2106269743
46.2049558527 2 9.51515800919
2
194.998750566 2 24.5700087185
31.6794513291 2 8.86648776669
Mn GEN 10 2
4
15.0 1 21.9061889166
328.592833748 3 21.3460106503
-162.049880237 2 21.2709151562
-1.85679609726 2 7.90771171833
2
244.669998154 2 18.9263045646
33.5399867643 2 8.31114792811
2
162.350195446 2 20.162449313
24.1593874179 2 7.79269955633
Ni GEN 10 2
4
18.0 1 37.839331506
681.107967108 3 23.875701156
-173.162219465 2 19.8803935987
0.34274858261 2 3.56565870195
2
91.6513902571 2 13.545098213
331.659352198 2 27.7907700999
2
7.5147228016 2 6.46792786898
265.586894944 2 23.6921476759
O GEN 2 1
3
6.0 1 12.30997
73.85984 3 14.76962
-47.876 2 13.71419
1
85.86406 2 13.65512
S GEN 2 1
3
14.00000000 1 17.46806994
244.55297916 3 16.40396851
-128.37752591 2 16.71429998
2
30.00006536 2 54.87912854
125.50010056 2 31.32968867
Sc GEN 10 2
4
11.0 1 16.0484863686
176.533350054 3 14.07764439
-83.673420518 2 11.993486653
0.331064789149 2 3.75115298216
2
153.959870288 2 11.4712713921
14.9643185607 2 5.00756742752
2
97.2094454291 2 11.4449481137
10.8162163087 2 4.78509457131
Ti GEN 10 2
4
12.0 1 18.4136620219
220.963944263 3 15.9229241432
-94.2902582468 2 13.6500062314
0.0979114248227 2 5.0955521057
2
173.946572359 2 12.7058061392
18.8376833381 2 6.11178551988
2
111.45672882 2 12.6409192965
11.1770268269 2 5.35437415684
V GEN 10 2
4
13.0 1 20.3216891426
264.181958854 3 19.5969804012
-115.292932083 2 17.3314734817
-0.662887260057 2 5.12320657929
2
195.567138911 2 15.1250215054
22.8864283476 2 6.2989891447
2
126.421195008 2 15.9385511327
16.0359712766 2 5.74006266866
Zn GEN 10 2
4
20.0 1 35.8079761618
716.159523235 3 34.536460837
-204.683933235 2 28.6283017827
0.760266144617 2 7.9623968256
2
95.8764043739 2 14.6349869153
431.708043027 2 35.0214135667
2
74.0127004894 2 14.5742930415
313.577705639 2 42.2297923499
Cu GEN 10 2
4
19.0 1 31.5381126304
599.224139977 3 31.0692553147
-244.689154841 2 30.5903586806
-1.2934952584 2 14.0514106386
2
66.2756081341 2 12.7723591969
370.71371825 2 29.355622426
2
49.7626505709 2 12.5247148487
271.662810283 2 33.5169454376

2
install/scripts/install_f77zmq.sh
View File

@ -14,7 +14,7 @@ function _install()
make -j 8 || exit 1
mv libf77zmq.a "${QP_ROOT}"/lib || exit 1
mv libf77zmq.so "${QP_ROOT}"/lib || exit 1
cp f77_zmq.h "${QP_ROOT}"/src/ZMQ/
cp f77_zmq_free.h "${QP_ROOT}"/src/ZMQ/f77_zmq.h
cd -
return 0
}

2
install/scripts/install_ocaml.sh
View File

@ -5,7 +5,7 @@ QP_ROOT=$PWD
cd -
# Normal installation
PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 ZMQ ppx_sexp_conv ppx_deriving jbuilder.1.0+beta17"
PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 zmq ppx_sexp_conv ppx_deriving jbuilder.1.0+beta17"
# Needed for ZeroMQ
export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}"

3
install/scripts/install_zeromq.sh
View File

@ -11,7 +11,8 @@ function _install()
set -u
ORIG=$(pwd)
cd "${BUILD}"
./configure --prefix=$QP_ROOT --without-libsodium || exit 1
./autogen.sh
./configure --prefix=$QP_ROOT --without-libsodium --disable-libunwind || exit 1
make -j 8 || exit 1
make install || exit 1
cd ${ORIG}

1
lib64
View File

@ -1 +0,0 @@
lib

3
ocaml/Makefile
View File

@ -12,7 +12,8 @@ endif
LIBS=
PKGS=
OCAMLCFLAGS="-g"
OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS)
OCAMLOPTFLAGS="opt -O3 -remove-unused-arguments -rounds 16 -inline 100 -inline-max-unroll 100"
OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) -ocamlopt $(OCAMLOPTFLAGS)
MLLFILES=$(wildcard *.mll)
MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml Input_auto_generated.ml qp_edit.ml
MLIFILES=$(wildcard *.mli) git

40
ocaml/Symmetry.ml
View File

@ -16,30 +16,30 @@ let to_string = function
| L -> "L"
let of_string = function
| "S" -> S
| "P" -> P
| "D" -> D
| "F" -> F
| "G" -> G
| "H" -> H
| "I" -> I
| "J" -> J
| "K" -> K
| "L" -> L
| "S" | "s" -> S
| "P" | "p" -> P
| "D" | "d" -> D
| "F" | "f" -> F
| "G" | "g" -> G
| "H" | "h" -> H
| "I" | "i" -> I
| "J" | "j" -> J
| "K" | "k" -> K
| "L" | "l" -> L
| x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L,
not "^x^"."))
let of_char = function
| 'S' -> S
| 'P' -> P
| 'D' -> D
| 'F' -> F
| 'G' -> G
| 'H' -> H
| 'I' -> I
| 'J' -> J
| 'K' -> K
| 'L' -> L
| 'S' | 's' -> S
| 'P' | 'p' -> P
| 'D' | 'd' -> D
| 'F' | 'f' -> F
| 'G' | 'g' -> G
| 'H' | 'h' -> H
| 'I' | 'i' -> I
| 'J' | 'j' -> J
| 'K' | 'k' -> K
| 'L' | 'l' -> L
| x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L"))
let to_l = function

108
ocaml/TaskServer.ml
View File

@ -47,10 +47,10 @@ let debug str =
let zmq_context =
ZMQ.Context.create ()
Zmq.Context.create ()
let () =
ZMQ.Context.set_io_threads zmq_context 8
Zmq.Context.set_io_threads zmq_context 8
let bind_socket ~socket_type ~socket ~port =
@ -61,7 +61,7 @@ let bind_socket ~socket_type ~socket ~port =
| -1 -> ()
| i ->
try
ZMQ.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port;
Zmq.Socket.bind socket @@ Printf.sprintf "tcp://*:%d" port;
loop (-1)
with
| Unix.Unix_error _ -> (Time.pause @@ Time.Span.of_sec 1. ; loop (i-1) )
@ -105,31 +105,31 @@ let ip_address = lazy (
let reply_ok rep_socket =
Message.Ok_msg.create
|> Message.Ok_msg.to_string
|> ZMQ.Socket.send rep_socket
|> Zmq.Socket.send rep_socket
let reply_wrong_state rep_socket =
Message.Error_msg.create "Wrong state"
|> Message.Error_msg.to_string
|> ZMQ.Socket.send rep_socket
|> Zmq.Socket.send rep_socket
let stop ~port =
debug "STOP";
let req_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.req
Zmq.Socket.create zmq_context Zmq.Socket.req
and address =
Printf.sprintf "tcp://localhost:%d" port
in
ZMQ.Socket.set_linger_period req_socket 1_000_000;
ZMQ.Socket.connect req_socket address;
Zmq.Socket.set_linger_period req_socket 1_000_000;
Zmq.Socket.connect req_socket address;
Message.Terminate (Message.Terminate_msg.create)
|> Message.to_string
|> ZMQ.Socket.send req_socket ;
|> Zmq.Socket.send req_socket ;
let msg =
ZMQ.Socket.recv req_socket
Zmq.Socket.recv req_socket
|> Message.of_string
in
let () =
@ -137,8 +137,8 @@ let stop ~port =
| Message.Ok _ -> ()
| _ -> failwith "Problem in termination"
in
ZMQ.Socket.set_linger_period req_socket 1_000;
ZMQ.Socket.close req_socket
Zmq.Socket.set_linger_period req_socket 1_000;
Zmq.Socket.close req_socket
let new_job msg program_state rep_socket pair_socket =
@ -166,7 +166,7 @@ let new_job msg program_state rep_socket pair_socket =
in
reply_ok rep_socket;
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket ;
|> Zmq.Socket.send pair_socket ;
result
let change_pub_state msg program_state rep_socket pair_socket =
@ -186,7 +186,7 @@ let change_pub_state msg program_state rep_socket pair_socket =
in
reply_ok rep_socket;
string_of_pub_state msg
|> ZMQ.Socket.send pair_socket ;
|> Zmq.Socket.send pair_socket ;
program_state
@ -216,7 +216,7 @@ let end_job msg program_state rep_socket pair_socket =
Printf.sprintf "waiting for %d slaves..." n
|> Message.Error_msg.create
|> Message.Error_msg.to_string
|> ZMQ.Socket.send rep_socket ;
|> Zmq.Socket.send rep_socket ;
program_state
in
@ -227,13 +227,13 @@ let end_job msg program_state rep_socket pair_socket =
if (msg.Message.Endjob_msg.state = force_state) then
begin
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket ;
|> Zmq.Socket.send pair_socket ;
success ()
end
else if (msg.Message.Endjob_msg.state = state) then
begin
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket ;
|> Zmq.Socket.send pair_socket ;
if (Queuing_system.number_of_clients program_state.queue = 0) then
success ()
else
@ -280,7 +280,7 @@ let connect msg program_state rep_socket =
Message.ConnectReply (Message.ConnectReply_msg.create
~state:state ~client_id ~push_address)
|> Message.to_string
|> ZMQ.Socket.send rep_socket ;
|> Zmq.Socket.send rep_socket ;
{ program_state with
queue = new_queue
}
@ -306,7 +306,7 @@ let disconnect msg program_state rep_socket =
in
Message.DisconnectReply (Message.DisconnectReply_msg.create ~state)
|> Message.to_string
|> ZMQ.Socket.send rep_socket ;
|> Zmq.Socket.send rep_socket ;
new_program_state
in
@ -352,7 +352,7 @@ let del_task msg program_state rep_socket =
in
Message.DelTaskReply (Message.DelTaskReply_msg.create ~task_ids ~more)
|> Message.to_string
|> ZMQ.Socket.send ~block:true rep_socket ; (** /!\ Has to be blocking *)
|> Zmq.Socket.send ~block:true rep_socket ; (** /!\ Has to be blocking *)
new_program_state
in
@ -426,10 +426,10 @@ let get_task msg program_state rep_socket pair_socket =
if no_task then
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket
|> Zmq.Socket.send pair_socket
else
string_of_pub_state (Running (Message.State.to_string state))
|> ZMQ.Socket.send pair_socket;
|> Zmq.Socket.send pair_socket;
let new_program_state =
{ program_state with
@ -440,7 +440,7 @@ let get_task msg program_state rep_socket pair_socket =
Message.GetTaskReply (Message.GetTaskReply_msg.create ~task ~task_id)
|> Message.to_string
|> ZMQ.Socket.send rep_socket ;
|> Zmq.Socket.send rep_socket ;
new_program_state
in
@ -498,10 +498,10 @@ let get_tasks msg program_state rep_socket pair_socket =
if no_task then
string_of_pub_state Waiting
|> ZMQ.Socket.send pair_socket
|> Zmq.Socket.send pair_socket
else
string_of_pub_state (Running (Message.State.to_string state))
|> ZMQ.Socket.send pair_socket;
|> Zmq.Socket.send pair_socket;
let new_program_state =
{ program_state with
@ -512,7 +512,7 @@ let get_tasks msg program_state rep_socket pair_socket =
Message.GetTasksReply (Message.GetTasksReply_msg.create result)
|> Message.to_string_list
|> ZMQ.Socket.send_all rep_socket ;
|> Zmq.Socket.send_all rep_socket ;
new_program_state
in
@ -596,7 +596,7 @@ let put_data msg rest_of_msg program_state rep_socket =
StringHashtbl.set program_state.data ~key ~data:value ;
Message.PutDataReply (Message.PutDataReply_msg.create ())
|> Message.to_string
|> ZMQ.Socket.send rep_socket;
|> Zmq.Socket.send rep_socket;
program_state
and failure () =
@ -625,11 +625,11 @@ let get_data msg program_state rep_socket =
let value =
match StringHashtbl.find program_state.data key with
| Some value -> value
| None -> ""
| None -> "\000"
in
Message.GetDataReply (Message.GetDataReply_msg.create ~value)
|> Message.to_string_list
|> ZMQ.Socket.send_all rep_socket;
|> Zmq.Socket.send_all rep_socket;
program_state
and failure () =
@ -699,7 +699,7 @@ let abort program_state rep_socket =
let error msg program_state rep_socket =
Message.Error (Message.Error_msg.create msg)
|> Message.to_string
|> ZMQ.Socket.send rep_socket ;
|> Zmq.Socket.send rep_socket ;
program_state
let start_pub_thread ~port =
@ -709,54 +709,54 @@ let start_pub_thread ~port =
in
let pair_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.pair
Zmq.Socket.create zmq_context Zmq.Socket.pair
and address =
"inproc://pair"
in
ZMQ.Socket.connect pair_socket address;
Zmq.Socket.connect pair_socket address;
let pub_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.pub
Zmq.Socket.create zmq_context Zmq.Socket.pub
in
bind_socket ~socket_type:"PUB" ~socket:pub_socket ~port;
let pollitem =
ZMQ.Poll.mask_of
[| (pair_socket, ZMQ.Poll.In) |]
Zmq.Poll.mask_of
[| (pair_socket, Zmq.Poll.In) |]
in
let rec run state =
let new_state =
let polling =
ZMQ.Poll.poll ~timeout pollitem
Zmq.Poll.poll ~timeout pollitem
in
if (polling.(0) = Some ZMQ.Poll.In) then
ZMQ.Socket.recv ~block:false pair_socket
if (polling.(0) = Some Zmq.Poll.In) then
Zmq.Socket.recv ~block:false pair_socket
|> pub_state_of_string
else
state
in
ZMQ.Socket.send pub_socket @@ string_of_pub_state new_state;
Zmq.Socket.send pub_socket @@ string_of_pub_state new_state;
match state with
| Stopped -> ()
| _ -> run new_state
in
run Waiting;
ZMQ.Socket.set_linger_period pair_socket 1000 ;
ZMQ.Socket.close pair_socket;
ZMQ.Socket.set_linger_period pub_socket 1000 ;
ZMQ.Socket.close pub_socket;
Zmq.Socket.set_linger_period pair_socket 1000 ;
Zmq.Socket.close pair_socket;
Zmq.Socket.set_linger_period pub_socket 1000 ;
Zmq.Socket.close pub_socket;
)
let run ~port =
(** Bind inproc socket for changing state of pub *)
let pair_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.pair
Zmq.Socket.create zmq_context Zmq.Socket.pair
and address =
"inproc://pair"
in
ZMQ.Socket.bind pair_socket address;
Zmq.Socket.bind pair_socket address;
let pub_thread =
start_pub_thread ~port:(port+1) ()
@ -764,9 +764,9 @@ let run ~port =
(** Bind REP socket *)
let rep_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.rep
Zmq.Socket.create zmq_context Zmq.Socket.rep
in
ZMQ.Socket.set_linger_period rep_socket 1_000_000;
Zmq.Socket.set_linger_period rep_socket 1_000_000;
bind_socket "REP" rep_socket port;
let initial_program_state =
@ -783,8 +783,8 @@ let run ~port =
(** ZMR polling item *)
let pollitem =
ZMQ.Poll.mask_of
[| (rep_socket, ZMQ.Poll.In) |]
Zmq.Poll.mask_of
[| (rep_socket, Zmq.Poll.In) |]
in
let address =
@ -798,9 +798,9 @@ let run ~port =
| false -> ()
| true ->
let polling =
ZMQ.Poll.poll ~timeout:1000 pollitem
Zmq.Poll.poll ~timeout:1000 pollitem
in
if (polling.(0) <> Some ZMQ.Poll.In) then
if (polling.(0) <> Some Zmq.Poll.In) then
main_loop program_state true
else
begin
@ -818,7 +818,7 @@ let run ~port =
(** Extract message *)
let raw_message, rest =
match ZMQ.Socket.recv_all rep_socket with
match Zmq.Socket.recv_all rep_socket with
| x :: rest -> x, rest
| [] -> failwith "Badly formed message"
in
@ -874,9 +874,9 @@ let run ~port =
end
in main_loop initial_program_state true;
ZMQ.Socket.send pair_socket @@ string_of_pub_state Stopped;
Zmq.Socket.send pair_socket @@ string_of_pub_state Stopped;
Thread.join pub_thread;
ZMQ.Socket.close rep_socket
Zmq.Socket.close rep_socket

32
ocaml/TaskServer.mli
View File

@ -19,14 +19,14 @@ val debug_env : bool
(** Print a debug message *)
val debug : string -> unit
(** {1} ZMQ *)
(** {1} Zmq *)
(** ZeroMQ context *)
val zmq_context : ZMQ.Context.t
val zmq_context : Zmq.Context.t
(** Bind a ZMQ socket to a TCP port and to an IPC file /tmp/qp_run.<port> *)
(** Bind a Zmq socket to a TCP port and to an IPC file /tmp/qp_run.<port> *)
val bind_socket :
socket_type:string -> socket:'a ZMQ.Socket.t -> port:int -> unit
socket_type:string -> socket:'a Zmq.Socket.t -> port:int -> unit
(** Name of the host on which the server runs *)
val hostname : string lazy_t
@ -35,8 +35,8 @@ val hostname : string lazy_t
val ip_address : string lazy_t
(** Standard messages *)
val reply_ok : [> `Req ] ZMQ.Socket.t -> unit
val reply_wrong_state : [> `Req ] ZMQ.Socket.t -> unit
val reply_ok : [> `Req ] Zmq.Socket.t -> unit
val reply_wrong_state : [> `Req ] Zmq.Socket.t -> unit
(** Stop server *)
val stop : port:int -> unit
@ -44,34 +44,34 @@ val stop : port:int -> unit
(** {1} Server functions *)
(** Create a new job *)
val new_job : Message.Newjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> [> `Pair] ZMQ.Socket.t -> t
val new_job : Message.Newjob_msg.t -> t -> [> `Req ] Zmq.Socket.t -> [> `Pair] Zmq.Socket.t -> t
(** Finish a running job *)
val end_job : Message.Endjob_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> [> `Pair] ZMQ.Socket.t -> t
val end_job : Message.Endjob_msg.t -> t -> [> `Req ] Zmq.Socket.t -> [> `Pair] Zmq.Socket.t -> t
(** Connect a client *)
val connect: Message.Connect_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
val connect: Message.Connect_msg.t -> t -> [> `Req ] Zmq.Socket.t -> t
(** Disconnect a client *)
val disconnect: Message.Disconnect_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
val disconnect: Message.Disconnect_msg.t -> t -> [> `Req ] Zmq.Socket.t -> t
(** Add a task to the pool *)
val add_task: Message.AddTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
val add_task: Message.AddTask_msg.t -> t -> [> `Req ] Zmq.Socket.t -> t
(** Mark the task as done by the client *)
val task_done: Message.TaskDone_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
val task_done: Message.TaskDone_msg.t -> t -> [> `Req ] Zmq.Socket.t -> t
(** Delete a task when it has been pulled by the collector *)
val del_task: Message.DelTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> t
val del_task: Message.DelTask_msg.t -> t -> [> `Req ] Zmq.Socket.t -> t
(** The client get a new task to execute *)
val get_task: Message.GetTask_msg.t -> t -> [> `Req ] ZMQ.Socket.t -> [> `Pair] ZMQ.Socket.t -> t
val get_task: Message.GetTask_msg.t -> t -> [> `Req ] Zmq.Socket.t -> [> `Pair] Zmq.Socket.t -> t
(** Terminate server *)
val terminate : t -> [> `Req ] ZMQ.Socket.t -> t
val terminate : t -> [> `Req ] Zmq.Socket.t -> t
(** Reply an Error message *)
val error : string -> t -> [> `Req ] ZMQ.Socket.t -> t
val error : string -> t -> [> `Req ] Zmq.Socket.t -> t
(** Run server *)
val run : port:int -> unit

2
ocaml/_tags
View File

@ -1,3 +1,3 @@
true: package(core,cryptokit,ZMQ,str,ppx_sexp_conv,ppx_deriving)
true: package(core,cryptokit,zmq,str,ppx_sexp_conv,ppx_deriving)
true: thread
false: profile

77
ocaml/qp_find_pi_space.ml
View File

@ -2,7 +2,7 @@ open Core
open Qputils
open Qptypes
let run ?(sym="None") ezfio_filename =
let run ?(sym="None") ?(apply="no") fc ezfio_filename =
Ezfio.set_file ezfio_filename ;
let aos =
@ -72,16 +72,61 @@ let run ?(sym="None") ezfio_filename =
| _ -> ([],[])
in
match sym with
| "x" | "X" | "y" | "Y" | "z" | "Z" ->
begin
Printf.printf "Pi: [";
List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi;
Printf.printf "\b]\n\nSigma: [";
List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma;
Printf.printf "\b]\n"
end
| _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result
begin
match sym with
| "x" | "X" | "y" | "Y" | "z" | "Z" ->
if apply = "cas" then
begin
let ne = Ezfio.get_electrons_elec_alpha_num () in
let command =
"qp_set_mo_class -act \"[" ^
(String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^
"]\" -core \"[" ^
(String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i <= ne) sigma) ^
"]\" -del \"[" ^
(String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > ne) sigma) ^
"]\" " ^ ezfio_filename
in
print_endline command;
if Sys.command command <> 0 then
failwith "Command failed"
end
else if apply = "cassd" then
begin
let ne = Ezfio.get_electrons_elec_alpha_num () in
let nfc =
if fc then
let ic =
Unix.open_process_in ("qp_set_frozen_core.py -q "^ezfio_filename)
in
let result = int_of_string @@ input_line ic in
close_in ic; result
else 0
in
let command =
"qp_set_mo_class" ^
( if nfc > 0 then
Printf.sprintf " -core \"[1-%d]\"" nfc else "") ^
" -inact \"[" ^
(String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > nfc && mo_i <= ne) sigma) ^ "]\"" ^
" -act \"[" ^ (String.concat ~sep:"," @@ List.map ~f:string_of_int pi) ^ "]\"" ^
" -virt \"[" ^
(String.concat ~sep:"," @@ List.map ~f:string_of_int @@ List.filter ~f:(fun mo_i -> mo_i > ne) sigma) ^ "]\" " ^ ezfio_filename
in
print_endline command;
if Sys.command command <> 0 then
failwith "Command failed"
end
else
begin
Printf.printf "Pi: [";
List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) pi;
Printf.printf "\b]\n\nSigma: [";
List.iter ~f:(fun mo_i -> Printf.printf "%d," mo_i) sigma;
Printf.printf "\b]\n"
end
| _ -> List.iter ~f:(fun (mo_i,x,y,z) -> Printf.printf "%d: (%f,%f,%f)\n" mo_i x y z) result
end
@ -90,10 +135,16 @@ let run ?(sym="None") ezfio_filename =
let spec =
let open Command.Spec in
empty
+> flag "sym" (optional string) ~doc:"{x,y,z} Axis perpendicular to the plane"
+> flag "sym" (optional string)
~doc:"{x,y,z} Axis perpendicular to the plane"
+> flag "apply" (optional string)
~doc:"[cas|cassd] Set the pi space as the active space, all other MOs frozen. If [cassd] is chosen, set inactive and virtual mos."
+> flag "fc" no_arg
~doc:"freeze core MOs with -apply cassd"
+> anon ("ezfio_filename" %: string)
let command =
Command.basic_spec
~summary: "Quantum Package command"
@ -101,7 +152,7 @@ let command =
"Find all the pi molecular orbitals to create a pi space.
")
spec
(fun sym ezfio_filename () -> run ?sym ezfio_filename )
(fun sym apply fc ezfio_filename () -> run ?sym ?apply fc ezfio_filename)
let () =

12
ocaml/qp_run.ml
View File

@ -20,10 +20,10 @@ let run slave exe ezfio_file =
(** Check availability of the ports *)
let port_number =
let zmq_context =
ZMQ.Context.create ()
Zmq.Context.create ()
in
let dummy_socket =
ZMQ.Socket.create zmq_context ZMQ.Socket.rep
Zmq.Socket.create zmq_context Zmq.Socket.rep
in
let rec try_new_port port_number =
try
@ -31,8 +31,8 @@ let run slave exe ezfio_file =
let address =
Printf.sprintf "tcp://%s:%d" (Lazy.force TaskServer.ip_address) (port_number+i)
in
ZMQ.Socket.bind dummy_socket address;
ZMQ.Socket.unbind dummy_socket address;
Zmq.Socket.bind dummy_socket address;
Zmq.Socket.unbind dummy_socket address;
);
port_number
with
@ -41,8 +41,8 @@ let run slave exe ezfio_file =
let result =
try_new_port 41279
in
ZMQ.Socket.close dummy_socket;
ZMQ.Context.terminate zmq_context;
Zmq.Socket.close dummy_socket;
Zmq.Context.terminate zmq_context;
result
in

4
ocaml/qptypes_generator.ml
View File

@ -84,8 +84,8 @@ let input_data = "
* MO_coef : float
* MO_occ : float
if (x < 0.) || (x > 2.) then
raise (Invalid_argument (Printf.sprintf \"MO_occ : (0. <= x <= 2.) : x=%f\" x));
if x < 0. then 0. else
if x > 2. then 2. else
* AO_coef : float

6
plugins/Bk/EZFIO.cfg
View File

@ -15,3 +15,9 @@ doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 10
[h0_type]
type: Perturbation
doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
interface: ezfio,provider,ocaml
default: EN

2
plugins/Bk/NEEDED_CHILDREN_MODULES
View File

@ -1,2 +1,2 @@
Bitmask dress_zmq DavidsonDressed
Bitmask dress_zmq DavidsonDressed Generators_full Selectors_full

3
plugins/Bk/dressing.irp.f
View File

@ -53,7 +53,9 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
do i_state=1,N_states
hdress = c_alpha(i_state) * i_h_alpha
sdress = c_alpha(i_state) * i_s_alpha
!$OMP ATOMIC
delta_ij_loc(i_state,j,1) = delta_ij_loc(i_state,j,1) + hdress
!$OMP ATOMIC
delta_ij_loc(i_state,j,2) = delta_ij_loc(i_state,j,2) + sdress
enddo
enddo
@ -62,4 +64,3 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
end subroutine

58
plugins/Bk/extra_functions.irp.f Normal file
View File

@ -0,0 +1,58 @@
BEGIN_PROVIDER [ integer, N_dress_int_buffer ]
&BEGIN_PROVIDER [ integer, N_dress_double_buffer ]
&BEGIN_PROVIDER [ integer, N_dress_det_buffer ]
implicit none
N_dress_int_buffer = 1
N_dress_double_buffer = 1
N_dress_det_buffer = 1
END_PROVIDER
subroutine delta_ij_done()
BEGIN_DOC
! This subroutine is executed on the master when the dressing has been computed,
! before the diagonalization.
END_DOC
end
subroutine dress_pulled(ind, int_buf, double_buf, det_buf, N_buf)
use bitmasks
implicit none
BEGIN_DOC
! Dress the contributions pulled from the slave.
END_DOC
integer, intent(in) :: ind, N_buf(3)
integer, intent(in) :: int_buf(*)
double precision, intent(in) :: double_buf(*)
integer(bit_kind), intent(in) :: det_buf(N_int,2,*)
end
subroutine generator_start(i_gen, iproc)
implicit none
BEGIN_DOC
! This subroutine is executed on the slave before computing the contribution of a generator.
END_DOC
integer, intent(in) :: i_gen, iproc
integer :: i
end
subroutine generator_done(i_gen, int_buf, double_buf, det_buf, N_buf, iproc)
implicit none
BEGIN_DOC
! This subroutine is executed on the slave after computing the contribution of a generator.
END_DOC
integer, intent(in) :: i_gen, iproc
integer, intent(out) :: int_buf(N_dress_int_buffer), N_buf(3)
double precision, intent(out) :: double_buf(N_dress_double_buffer)
integer(bit_kind), intent(out) :: det_buf(N_int, 2, N_dress_det_buffer)
N_buf(:) = 1
int_buf(:) = 0
double_buf(:) = 0.d0
det_buf(:,:,:) = 0
end

76
plugins/CAS_SD_ZMQ/cassd_zmq.irp.f
View File

@ -5,16 +5,28 @@ program cassd_zmq
integer :: degree
integer :: n_det_before, to_select
double precision :: threshold_davidson_in
double precision :: error(N_states)
allocate (pt2(N_states))
double precision :: hf_energy_ref
logical :: has
integer :: N_states_p
character*(512) :: fmt
character*(8) :: pt2_string
pt2 = -huge(1.d0)
error = 0.d0
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
if (do_pt2) then
pt2_string = ' '
else
pt2_string = '(approx)'
endif
call diagonalize_CI
call save_wavefunction
@ -45,7 +57,6 @@ program cassd_zmq
double precision :: E_CI_before(N_states)
print*,'Beginning the selection ...'
if (.True.) then ! Avoid pre-calculation of CI_energy
E_CI_before(1:N_states) = CI_energy(1:N_states)
endif
@ -60,6 +71,8 @@ program cassd_zmq
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. &
(correlation_energy_ratio <= correlation_energy_ratio_max) &
)
write(*,'(A)') '--------------------------------------------------------------------------------'
correlation_energy_ratio = (CI_energy(1) - hf_energy_ref) / &
(E_CI_before(1) + pt2(1) - hf_energy_ref)
@ -98,6 +111,67 @@ program cassd_zmq
to_select = min(to_select, N_det_max-n_det_before)
call ZMQ_selection(to_select, pt2)
N_states_p = min(N_det,N_states)
print *, ''
print '(A,I12)', 'Summary at N_det = ', N_det
print '(A)', '-----------------------------------'
print *, ''
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
print *, ''
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
write(*,fmt) ('State',k, k=1,N_states_p)
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
write(*,fmt) '# E ', E_CI_before(1:N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', E_CI_before(1:N_states_p)-E_CI_before(1)
write(*,fmt) '# Excit. (eV)', (E_CI_before(1:N_states_p)-E_CI_before(1))*27.211396641308d0
endif
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (E_CI_before(k)+pt2(k),error(k), k=1,N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1)), &
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
write(*,fmt) '# Excit. (eV)', ( (E_CI_before(k)+pt2(k)-E_CI_before(1)-pt2(1))*27.211396641308d0, &
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
endif
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
print *, ''
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print*, 'correlation_ratio = ', correlation_energy_ratio
do k=1, N_states_p
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', E_CI_before(k)
print *, 'E+PT2'//pt2_string//' = ', E_CI_before(k)+pt2(k), ' +/- ', error(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print *, 'Variational Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (E_CI_before(i) - E_CI_before(1)), &
(E_CI_before(i) - E_CI_before(1)) * 27.211396641308d0
enddo
print *, '-----'
print*, 'Variational + perturbative Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))), &
(E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1))) * 27.211396641308d0
enddo
endif
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted

4
plugins/CISD/H_apply.irp.f
View File

@ -3,7 +3,9 @@
BEGIN_SHELL [ /usr/bin/env python2 ]
from generate_h_apply import H_apply
H = H_apply("cisd")
H = H_apply("cisd",do_double_exc=True)
print H
H = H_apply("cis",do_double_exc=False)
print H
END_SHELL

18
plugins/CISD/cisd.irp.f
View File

@ -1,20 +1,22 @@
program cisd
program cis
implicit none
read_wf = .False.
SOFT_TOUCH read_wf
call run
end
subroutine run
implicit none
integer :: i
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
N_det = 1
touch psi_det psi_coef N_det
call H_apply_cisd
print *, 'N_det = ', N_det
do i = 1,N_states
print *, 'energy = ',CI_energy(i)
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
enddo
call ezfio_set_cisd_energy(CI_energy(1))
psi_coef = ci_eigenvectors
SOFT_TOUCH psi_coef
call save_wavefunction
call save_wavefunction_truncated(1.d-12)
end

4
plugins/DavidsonDressed/diagonalize_CI.irp.f
View File

@ -63,8 +63,8 @@ END_PROVIDER
call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,&
size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,&
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1)
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors_dressed,1))
! call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
! N_states_diag,size(CI_eigenvectors_dressed,1))
else if (diag_algorithm == "Lapack") then

14
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -22,8 +22,8 @@ program fci_zmq
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
call diagonalize_CI
call save_wavefunction
! call diagonalize_CI
! call save_wavefunction
call ezfio_has_hartree_fock_energy(has)
if (has) then
@ -52,6 +52,11 @@ program fci_zmq
double precision :: error(N_states)
correlation_energy_ratio = 0.d0
if (do_pt2) then
pt2_string = ' '
else
pt2_string = '(approx)'
endif
if (.True.) then ! Avoid pre-calculation of CI_energy
do while ( &
@ -63,17 +68,14 @@ program fci_zmq
if (do_pt2) then
pt2_string = ' '
pt2 = 0.d0
threshold_selectors = 1.d0
threshold_generators = 1d0
threshold_generators = 1.d0
SOFT_TOUCH threshold_selectors threshold_generators
call ZMQ_pt2(CI_energy, pt2,relative_error,absolute_error,error) ! Stochastic PT2
threshold_selectors = threshold_selectors_save
threshold_generators = threshold_generators_save
SOFT_TOUCH threshold_selectors threshold_generators
else
pt2_string = '(approx)'
endif

80
plugins/Full_CI_ZMQ/pt2_slave.irp.f
View File

@ -1,80 +0,0 @@
program pt2_slave
implicit none
BEGIN_DOC
! Helper program to compute the PT2 in distributed mode.
END_DOC
read_wf = .False.
distributed_davidson = .False.
SOFT_TOUCH read_wf distributed_davidson
call provide_everything
call switch_qp_run_to_master
call run_wf
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
end
subroutine run_wf
use f77_zmq
implicit none
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: states(1)
integer :: rc, i
integer, external :: zmq_get_dvector
integer, external :: zmq_get_psi
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'pt2'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,1)
if(trim(zmq_state) == 'Stopped') then
exit
else if (trim(zmq_state) == 'pt2') then
! Selection
! ---------
print *, 'PT2'
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call pt2_slave_tcp(i, energy)
!$OMP END PARALLEL
print *, 'PT2 done'
endif
end do
end
subroutine pt2_slave_tcp(i,energy)
implicit none
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: i
logical :: lstop
lstop = .False.
call run_pt2_slave(0,i,energy,lstop)
end

12
plugins/Full_CI_ZMQ/pt2_stoch.irp.f
View File

@ -3,6 +3,7 @@ program pt2_stoch
read_wf = .True.
SOFT_TOUCH read_wf
PROVIDE mo_bielec_integrals_in_map
PROVIDE psi_energy
call run
end
@ -19,14 +20,13 @@ subroutine run
double precision :: E_CI_before, relative_error, absolute_error, eqt
allocate (pt2(N_states))
call diagonalize_CI()
pt2 = 0.d0
pt2(:) = 0.d0
E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion
E_CI_before = psi_energy(1) + nuclear_repulsion
threshold_selectors = 1.d0
threshold_generators = 1d0
relative_error = 1.d-9
absolute_error = 1.d-9
threshold_generators = 1.d0
relative_error=PT2_relative_error
absolute_error=PT2_absolute_error
call ZMQ_pt2(E_CI_before, pt2, relative_error, absolute_error, eqt)
print *, 'Final step'

31
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

@ -95,10 +95,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
state_average_weight_save(:) = state_average_weight(:)
do pt2_stoch_istate=1,N_states
SOFT_TOUCH pt2_stoch_istate
state_average_weight(:) = 0.d0
state_average_weight(pt2_stoch_istate) = 1.d0
TOUCH state_average_weight
TOUCH state_average_weight pt2_stoch_istate
provide nproc pt2_F mo_bielec_integrals_in_map mo_mono_elec_integral pt2_w psi_selectors
@ -112,6 +111,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
integer, external :: zmq_put_N_det_generators
integer, external :: zmq_put_N_det_selectors
integer, external :: zmq_put_dvector
integer, external :: zmq_put_ivector
if (zmq_put_psi(zmq_to_qp_run_socket,1) == -1) then
stop 'Unable to put psi on ZMQ server'
endif
@ -124,6 +124,19 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then
stop 'Unable to put energy on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
stop 'Unable to put state_average_weight on ZMQ server'
endif
if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then
stop 'Unable to put pt2_stoch_istate on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) then
stop 'Unable to put threshold_selectors on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then
stop 'Unable to put threshold_generators on ZMQ server'
endif
integer, external :: add_task_to_taskserver
@ -143,7 +156,17 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
endif
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc+1) &
integer :: nproc_target
nproc_target = nproc
double precision :: mem
mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3)
call write_double(6,mem,'Estimated memory/thread (Gb)')
if (qp_max_mem > 0) then
nproc_target = max(1,int(dble(qp_max_mem)/mem))
nproc_target = min(nproc_target,nproc)
endif
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) &
!$OMP PRIVATE(i)
i = omp_get_thread_num()
if (i==0) then
@ -255,7 +278,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, absolute_error, pt2
end do
avg = S(t) / dble(c)
eqt = (S2(t) / c) - (S(t)/c)**2
eqt = sqrt(eqt / dble(c-1))
eqt = sqrt(eqt / dble(c-1+1))
pt2(pt2_stoch_istate) = E0-E+avg
error(pt2_stoch_istate) = eqt
time = omp_get_wtime()

23
plugins/Full_CI_ZMQ/run_pt2_slave.irp.f
View File

@ -39,13 +39,14 @@ subroutine run_pt2_slave(thread,iproc,energy)
zmq_socket_push = new_zmq_push_socket(thread)
buf%N = 0
n_tasks = 0
n_tasks = 1
call create_selection_buffer(0, 0, buf)
done = .False.
do while (.not.done)
n_tasks = min(n_tasks+1,pt2_n_tasks_max)
n_tasks = max(1,n_tasks)
n_tasks = min(n_tasks,pt2_n_tasks_max)
integer, external :: get_tasks_from_taskserver
if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then
@ -59,25 +60,35 @@ subroutine run_pt2_slave(thread,iproc,energy)
read (task(k),*) subset(k), i_generator(k)
enddo
double precision :: time0, time1
call wall_time(time0)
do k=1,n_tasks
pt2(:,k) = 0.d0
buf%cur = 0
call select_connected(i_generator(k),energy,pt2(1,k),buf,subset(k),pt2_F(i_generator(k)))
enddo
call wall_time(time1)
integer, external :: tasks_done_to_taskserver
if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then
done = .true.
endif
call push_pt2_results(zmq_socket_push, i_generator, pt2, task_id, n_tasks)
! Try to adjust n_tasks around 1 second per job
n_tasks = min(n_tasks,int( 1.d0*dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
! n_tasks = n_tasks+1
end do
integer, external :: disconnect_from_taskserver
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
continue
endif
do i=1,300
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) /= -2) exit
call sleep(1)
print *, 'Retry disconnect...'
end do
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call end_zmq_push_socket(zmq_socket_push,thread)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call delete_selection_buffer(buf)
end subroutine

125
plugins/Full_CI_ZMQ/run_selection_slave.irp.f
View File

@ -1,5 +1,122 @@
subroutine run_selection_slave(thread,iproc,energy)
implicit none
integer, intent(in) :: thread, iproc
double precision, intent(in) :: energy(N_states_diag)
call run_selection_slave_new(thread,iproc,energy)
end
subroutine run_selection_slave_new(thread,iproc,energy)
use f77_zmq
use selection_types
implicit none
integer, intent(in) :: thread, iproc
double precision, intent(in) :: energy(N_states_diag)
integer :: rc, i, N
logical :: buffer_ready
integer :: worker_id, ltask
character*(512), allocatable :: task(:)
integer, allocatable :: task_id(:)
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), external :: new_zmq_push_socket
integer(ZMQ_PTR) :: zmq_socket_push
type(selection_buffer) :: buf, buf2
logical :: done
double precision,allocatable :: pt2(:,:)
integer :: n_tasks, k, n_tasks_max
integer, allocatable :: i_generator(:), subset(:)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order
buffer_ready = .False.
n_tasks_max = N_det_generators/100+1
allocate(task_id(n_tasks_max), task(n_tasks_max))
allocate(pt2(N_states,n_tasks_max), i_generator(n_tasks_max), subset(n_tasks_max))
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
integer, external :: connect_to_taskserver
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
return
endif
zmq_socket_push = new_zmq_push_socket(thread)
buf%N = 0
n_tasks = 1
call create_selection_buffer(0, 0, buf)
done = .False.
do while (.not.done)
n_tasks = max(1,n_tasks)
n_tasks = min(n_tasks,n_tasks_max)
integer, external :: get_tasks_from_taskserver
if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then
exit
endif
done = task_id(n_tasks) == 0
if (done) n_tasks = n_tasks-1
if (n_tasks == 0) exit
do k=1,n_tasks
read (task(k),*) subset(k), i_generator(k), N
enddo
if(buf%N == 0) then
! Only first time
call create_selection_buffer(N, N*2, buf)
call create_selection_buffer(N, N*2, buf2)
buffer_ready = .True.
endif
double precision :: time0, time1
call wall_time(time0)
do k=1,n_tasks
pt2(:,k) = 0.d0
buf%cur = 0
call select_connected(i_generator(k),energy,pt2(1,k),buf,subset(k))
enddo
call wall_time(time1)
integer, external :: tasks_done_to_taskserver
if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then
done = .true.
endif
call sort_selection_buffer(buf)
call merge_selection_buffers(buf,buf2)
call push_selection_results(zmq_socket_push, pt2, buf, task_id, n_tasks)
buf%mini = buf2%mini
pt2(:,:) = 0d0
buf%cur = 0
! ! Try to adjust n_tasks around 5 second per job
! n_tasks = min(n_tasks,int( 5.d0 * dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
n_tasks = n_tasks+1
end do
integer, external :: disconnect_from_taskserver
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
continue
endif
call end_zmq_push_socket(zmq_socket_push,thread)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
call delete_selection_buffer(buf)
end
subroutine run_selection_slave_old(thread,iproc,energy)
use f77_zmq
use selection_types
implicit none
@ -50,8 +167,8 @@ subroutine run_selection_slave(thread,iproc,energy)
if (done) then
ctask = ctask - 1
else
integer :: i_generator, N
read(task,*) i_generator, N
integer :: i_generator, N, subset
read(task,*) subset, i_generator, N
if(buf%N == 0) then
! Only first time
call create_selection_buffer(N, N*2, buf)
@ -60,7 +177,7 @@ subroutine run_selection_slave(thread,iproc,energy)
else
ASSERT (N == buf%N)
end if
call select_connected(i_generator,energy,pt2,buf,1,1)
call select_connected(i_generator,energy,pt2,buf,subset,fragment_count)
endif
integer, external :: task_done_to_taskserver

176
plugins/Full_CI_ZMQ/selection.irp.f
View File

@ -137,7 +137,7 @@ subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
puti = p(j, sp)
if(bannedOrb(puti)) cycle
pmob = p(turn2(j), sp)
hij = mo_bielec_integral(pfix, pmob, hfix, hmob)
hij = mo_bielec_integral(pmob, pfix, hmob, hfix)
hij *= get_phase_bi(phasemask, sp, sfix, hmob, pmob, hfix, pfix)
vect(:, puti) += hij * coefs
end do
@ -191,17 +191,26 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
p2 = p(2, sp)
lbanned(p2) = .true.
double precision :: hij_cache(mo_tot_num,2)
call get_mo_bielec_integrals(hole,p1,p2,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(hole,p2,p1,mo_tot_num,hij_cache(1,2),mo_integrals_map)
do i=1,hole-1
if(lbanned(i)) cycle
hij = (mo_bielec_integral(p1, p2, i, hole) - mo_bielec_integral(p2, p1, i, hole))
hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
hij = hij_cache(i,1)-hij_cache(i,2)
if (hij /= 0.d0) then
hij *= get_phase_bi(phasemask, sp, sp, i, p1, hole, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
endif
end do
do i=hole+1,mo_tot_num
if(lbanned(i)) cycle
hij = (mo_bielec_integral(p1, p2, hole, i) - mo_bielec_integral(p2, p1, hole, i))
hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
hij = hij_cache(i,2)-hij_cache(i,1)
if (hij /= 0.d0) then
hij *= get_phase_bi(phasemask, sp, sp, hole, p1, i, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
endif
end do
call apply_particle(mask, sp, p2, det, ok, N_int)
@ -209,11 +218,14 @@ subroutine get_m1(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
vect(1:N_states, p2) += hij * coefs(1:N_states)
else
p2 = p(1, sh)
call get_mo_bielec_integrals(hole,p1,p2,mo_tot_num,hij_cache(1,1),mo_integrals_map)
do i=1,mo_tot_num
if(lbanned(i)) cycle
hij = mo_bielec_integral(p1, p2, i, hole)
hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
hij = hij_cache(i,1)
if (hij /= 0.d0) then
hij *= get_phase_bi(phasemask, sp, sh, i, p1, hole, p2)
vect(1:N_states,i) += hij * coefs(1:N_states)
endif
end do
end if
deallocate(lbanned)
@ -594,7 +606,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
logical :: ok
integer :: s1, s2, p1, p2, ib, j, istate
integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
double precision :: e_pert, delta_E, val, Hii, min_e_pert,tmp
double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp
double precision, external :: diag_H_mat_elem_fock
logical, external :: detEq
@ -621,7 +633,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
min_e_pert = 0d0
sum_e_pert = 0d0
do istate=1,N_states
delta_E = E0(istate) - Hii
@ -632,11 +644,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
endif
e_pert = 0.5d0 * (tmp - delta_E)
pt2(istate) = pt2(istate) + e_pert
min_e_pert = min(e_pert,min_e_pert)
sum_e_pert = sum_e_pert + e_pert * state_average_weight(istate)
end do
if(min_e_pert <= buf%mini) then
call add_to_selection_buffer(buf, det, min_e_pert)
if(sum_e_pert <= buf%mini) then
call add_to_selection_buffer(buf, det, sum_e_pert)
end if
end do
end do
@ -690,8 +702,8 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
if (interesting(i) == i_gen) then
if(sp == 3) then
do j=1,mo_tot_num
do k=1,mo_tot_num
do k=1,mo_tot_num
do j=1,mo_tot_num
banned(j,k,2) = banned(k,j,1)
enddo
enddo
@ -707,17 +719,17 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
call bitstring_to_list_in_selection(mobMask(1,2), p(1,2), p(0,2), N_int)
perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
do j=2,N_int
perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
end do
call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
if (interesting(i) >= i_gen) then
perMask(1,1) = iand(mask(1,1), not(det(1,1,i)))
perMask(1,2) = iand(mask(1,2), not(det(1,2,i)))
do j=2,N_int
perMask(j,1) = iand(mask(j,1), not(det(j,1,i)))
perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
end do
call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
call bitstring_to_list_in_selection(perMask(1,2), h(1,2), h(0,2), N_int)
call get_mask_phase(psi_det_sorted(1,1,interesting(i)), phasemask)
if(nt == 4) then
call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
@ -731,6 +743,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
if(nt == 3) call past_d1(bannedOrb, p)
end if
end do
end
@ -881,9 +894,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
double precision, allocatable :: hij_cache(:,:)
PROVIDE mo_integrals_map
allocate (lbanned(mo_tot_num, 2))
allocate (hij_cache(mo_tot_num,2))
lbanned = bannedOrb
do i=1, p(0,1)
@ -907,16 +922,26 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
call get_mo_bielec_integrals(hfix,p1,p2,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(hfix,p2,p1,mo_tot_num,hij_cache(1,2),mo_integrals_map)
tmp_row = 0d0
do putj=1, hfix-1
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
endif
end do
do putj=hfix+1, mo_tot_num
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states)
endif
end do
if(ma == 1) then
@ -930,19 +955,27 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
pfix = p(1,mi)
tmp_row = 0d0
tmp_row2 = 0d0
call get_mo_bielec_integrals(hfix,pfix,p1,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(hfix,pfix,p2,mo_tot_num,hij_cache(1,2),mo_integrals_map)
do puti=1,mo_tot_num
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
if(.not. banned(putj,puti,bant)) then
hij = mo_bielec_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
tmp_row(:,puti) += hij * coefs(:)
hij = hij_cache(puti,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix)
tmp_row(:,puti) += hij * coefs(:)
endif
end if
putj = p2
if(.not. banned(putj,puti,bant)) then
hij = mo_bielec_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
tmp_row2(:,puti) += hij * coefs(:)
hij = hij_cache(puti,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix)
tmp_row2(:,puti) += hij * coefs(:)
endif
end if
end do
@ -953,23 +986,35 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
mat(:,p1,:) += tmp_row(:,:)
mat(:,p2,:) += tmp_row2(:,:)
end if
else
else ! sp /= 3
if(p(0,ma) == 3) then
do i=1,3
hfix = h(1,ma)
puti = p(i, ma)
p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma)
call get_mo_bielec_integrals(hfix,p1,p2,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(hfix,p2,p1,mo_tot_num,hij_cache(1,2),mo_integrals_map)
tmp_row = 0d0
do putj=1,hfix-1
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_bielec_integral(p1, p2, putj, hfix)-mo_bielec_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
tmp_row(:,putj) += hij * coefs(:)
if(lbanned(putj,ma)) cycle
if(banned(putj,puti,1)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2)
tmp_row(:,putj) += hij * coefs(:)
endif
end do
do putj=hfix+1,mo_tot_num
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_bielec_integral(p1, p2, hfix, putj)-mo_bielec_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
tmp_row(:,putj) += hij * coefs(:)
if(lbanned(putj,ma)) cycle
if(banned(putj,puti,1)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2)
tmp_row(:,putj) += hij * coefs(:)
endif
end do
mat(:, :puti-1, puti) += tmp_row(:,:puti-1)
@ -982,18 +1027,26 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(2,ma)
tmp_row = 0d0
tmp_row2 = 0d0
call get_mo_bielec_integrals(hfix,p1,pfix,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(hfix,p2,pfix,mo_tot_num,hij_cache(1,2),mo_integrals_map)
do puti=1,mo_tot_num
if(lbanned(puti,ma)) cycle
putj = p2
if(.not. banned(puti,putj,1)) then
hij = mo_bielec_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
tmp_row(:,puti) += hij * coefs(:)
hij = hij_cache(puti,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1)
tmp_row(:,puti) += hij * coefs(:)
endif
end if
putj = p1
if(.not. banned(puti,putj,1)) then
hij = mo_bielec_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
tmp_row2(:,puti) += hij * coefs(:)
hij = hij_cache(puti,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2)
tmp_row2(:,puti) += hij * coefs(:)
endif
end if
end do
mat(:,:p2-1,p2) += tmp_row(:,:p2-1)
@ -1002,7 +1055,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
mat(:,p1,p1:) += tmp_row2(:,p1:)
end if
end if
deallocate(lbanned)
deallocate(lbanned,hij_cache)
!! MONO
if(sp == 3) then
@ -1016,10 +1069,11 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
do i1=1,p(0,s1)
ib = 1
if(s1 == s2) ib = i1+1
p1 = p(i1,s1)
if(bannedOrb(p1, s1)) cycle
do i2=ib,p(0,s2)
p1 = p(i1,s1)
p2 = p(i2,s2)
if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
if(bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
mat(:, p1, p2) += coefs(:) * hij
@ -1048,6 +1102,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
logical :: ok
integer :: bant
double precision, allocatable :: hij_cache(:,:)
allocate (hij_cache(mo_tot_num,2))
bant = 1
@ -1056,6 +1112,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
h2 = p(1,2)
do p1=1, mo_tot_num
if(bannedOrb(p1, 1)) cycle
call get_mo_bielec_integrals(p1,h2,h1,mo_tot_num,hij_cache(1,1),mo_integrals_map)
do p2=1, mo_tot_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
@ -1064,7 +1121,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
call i_h_j(gen, det, N_int, hij)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2)
hij = mo_bielec_integral(p1, p2, h1, h2) * phase
hij = hij_cache(p2,1) * phase
end if
mat(:, p1, p2) += coefs(:) * hij
end do
@ -1074,19 +1131,28 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p2 = p(2,sp)
do puti=1, mo_tot_num
if(bannedOrb(puti, sp)) cycle
call get_mo_bielec_integrals(puti,p2,p1,mo_tot_num,hij_cache(1,1),mo_integrals_map)
call get_mo_bielec_integrals(puti,p1,p2,mo_tot_num,hij_cache(1,2),mo_integrals_map)
do putj=puti+1, mo_tot_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
if (hij /= 0.d0) then
mat(:, puti, putj) += coefs(:) * hij
endif
else
hij = (mo_bielec_integral(p1, p2, puti, putj) - mo_bielec_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
hij = hij_cache(putj,1) - hij_cache(putj,2)
if (hij /= 0.d0) then
hij *= get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2)
mat(:, puti, putj) += coefs(:) * hij
endif
end if
mat(:, puti, putj) += coefs(:) * hij
end do
end do
end if
deallocate(hij_cache)
end

98
plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f
View File

@ -30,10 +30,12 @@ subroutine run_wf
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states)
character*(64) :: states(3)
character*(64) :: old_state
integer :: rc, i, ierr
double precision :: t0, t1
integer, external :: zmq_get_dvector, zmq_get_N_det_generators
integer, external :: zmq_get_ivector
integer, external :: zmq_get_psi, zmq_get_N_det_selectors
integer, external :: zmq_get_N_states_diag
@ -43,29 +45,59 @@ subroutine run_wf
states(1) = 'selection'
states(2) = 'davidson'
states(3) = 'pt2'
old_state = 'Waiting'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,size(states))
print *, trim(zmq_state)
if (mpi_master) then
call wait_for_states(states,zmq_state,size(states))
if (zmq_state(1:64) == old_state(1:64)) then
call sleep(1)
cycle
else
old_state(1:64) = zmq_state(1:64)
endif
print *, trim(zmq_state)
endif
IRP_IF MPI
call MPI_BCAST (zmq_state, 128, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in broadcast of zmq_state'
endif
IRP_ENDIF
if(zmq_state(1:7) == 'Stopped') then
exit
endif
else if (zmq_state(1:9) == 'selection') then
if (zmq_state(1:9) == 'selection') then
! Selection
! ---------
call wall_time(t0)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy state_average_weight threshold_selectors threshold_generators
if (mpi_master) then
print *, 'N_det', N_det
print *, 'N_det_generators', N_det_generators
print *, 'N_det_selectors', N_det_selectors
print *, 'psi_energy', psi_energy
print *, 'pt2_stoch_istate', pt2_stoch_istate
print *, 'state_average_weight', state_average_weight
endif
call wall_time(t1)
call write_double(6,(t1-t0),'Broadcast time')
@ -74,37 +106,65 @@ subroutine run_wf
call run_selection_slave(0,i,energy)
!$OMP END PARALLEL
print *, 'Selection done'
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
print *, 'All selection done'
else if (zmq_state(1:8) == 'davidson') then
! Davidson
! --------
print *, 'Davidson'
call wall_time(t0)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states_diag) == -1) cycle
call wall_time(t1)
call write_double(6,(t1-t0),'Broadcast time')
if (mpi_master) then
call write_double(6,(t1-t0),'Broadcast time')
endif
call omp_set_nested(.True.)
call davidson_slave_tcp(0)
call omp_set_nested(.False.)
print *, 'Davidson done'
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
print *, 'All Davidson done'
else if (zmq_state(1:3) == 'pt2') then
! PT2
! ---
print *, 'PT2'
call wall_time(t0)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy state_average_weight pt2_stoch_istate threshold_selectors threshold_generators
if (mpi_master) then
print *, 'N_det', N_det
print *, 'N_det_generators', N_det_generators
print *, 'N_det_selectors', N_det_selectors
print *, 'psi_energy', psi_energy
print *, 'pt2_stoch_istate', pt2_stoch_istate
print *, 'state_average_weight', state_average_weight
endif
call wall_time(t1)
call write_double(6,(t1-t0),'Broadcast time')
@ -113,22 +173,24 @@ subroutine run_wf
lstop = .False.
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_pt2_slave(0,i,energy,lstop)
call run_pt2_slave(0,i,pt2_e0_denominator)
!$OMP END PARALLEL
print *, 'PT2 done'
FREE state_average_weight
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
print *, 'All PT2 done'
endif
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
end do
IRP_IF MPI
call MPI_finalize(i)
call MPI_finalize(ierr)
IRP_ENDIF
end

110
plugins/Full_CI_ZMQ/selection_slave.irp.f
View File

@ -13,23 +13,31 @@ program selection_slave
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context n_states_diag
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count ci_energy mpi_master zmq_state zmq_context
PROVIDE psi_det psi_coef
end
subroutine run_wf
use f77_zmq
implicit none
IRP_IF MPI
include 'mpif.h'
IRP_ENDIF
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states)
character*(64) :: states(4)
character*(64) :: states(3)
character*(64) :: old_state
integer :: rc, i, ierr
double precision :: t0, t1
integer, external :: zmq_get_dvector
integer, external :: zmq_get_psi
integer, external :: zmq_get_dvector, zmq_get_N_det_generators
integer, external :: zmq_get_ivector
integer, external :: zmq_get_psi, zmq_get_N_det_selectors
integer, external :: zmq_get_N_states_diag
call provide_everything
@ -37,52 +45,112 @@ subroutine run_wf
states(1) = 'selection'
states(2) = 'davidson'
states(3) = 'pt2'
old_state = 'Waiting'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,3)
if (mpi_master) then
call wait_for_states(states,zmq_state,size(states))
if (zmq_state(1:64) == old_state(1:64)) then
call sleep(1)
cycle
else
old_state(1:64) = zmq_state(1:64)
endif
print *, trim(zmq_state)
endif
if(trim(zmq_state) == 'Stopped') then
IRP_IF MPI
call MPI_BCAST (zmq_state, 128, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in broadcast of zmq_state'
endif
IRP_ENDIF
if(zmq_state(1:7) == 'Stopped') then
exit
endif
else if (trim(zmq_state) == 'selection') then
if (zmq_state(1:9) == 'selection') then
! Selection
! ---------
print *, 'Selection'
call wall_time(t0)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy state_average_weight threshold_selectors threshold_generators
if (mpi_master) then
print *, 'N_det', N_det
print *, 'N_det_generators', N_det_generators
print *, 'N_det_selectors', N_det_selectors
print *, 'psi_energy', psi_energy
print *, 'pt2_stoch_istate', pt2_stoch_istate
print *, 'state_average_weight', state_average_weight
endif
call wall_time(t1)
call write_double(6,(t1-t0),'Broadcast time')
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_selection_slave(0,i,energy)
!$OMP END PARALLEL
print *, 'Selection done'
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
print *, 'All selection done'
if (N_det < 100000) then
exit
endif
else if (trim(zmq_state) == 'pt2') then
else if (zmq_state(1:8) == 'davidson') then
! PT2
! ---
! Davidson
! --------
print *, 'PT2'
call wall_time(t0)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states_diag) == -1) cycle
logical :: lstop
lstop = .False.
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_pt2_slave(0,i,energy,lstop)
!$OMP END PARALLEL
print *, 'PT2 done'
call wall_time(t1)
if (mpi_master) then
call write_double(6,(t1-t0),'Broadcast time')
endif
call omp_set_nested(.True.)
call davidson_slave_tcp(0)
call omp_set_nested(.False.)
print *, 'Davidson done'
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
print *, 'All Davidson done'
exit
endif
end do
IRP_IF MPI
call MPI_finalize(ierr)
IRP_ENDIF
end

70
plugins/Full_CI_ZMQ/zmq_selection.irp.f
View File

@ -10,7 +10,6 @@ subroutine ZMQ_selection(N_in, pt2)
integer :: i, N
integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states)
integer, parameter :: maxtasks=10000
PROVIDE fragment_count
@ -21,7 +20,7 @@ subroutine ZMQ_selection(N_in, pt2)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order
PROVIDE psi_bilinear_matrix_transp_order fragment_count
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
@ -42,31 +41,55 @@ subroutine ZMQ_selection(N_in, pt2)
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'energy',pt2_e0_denominator,size(pt2_e0_denominator)) == -1) then
stop 'Unable to put energy on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) then
stop 'Unable to put threshold_selectors on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
stop 'Unable to put state_average_weight on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then
stop 'Unable to put threshold_generators on ZMQ server'
endif
call create_selection_buffer(N, N*2, b)
endif
integer, external :: add_task_to_taskserver
character*(20*maxtasks) :: task
character(len=64000) :: task
integer :: j,k,ipos
ipos=1
task = ' '
integer :: k
k=0
do i= 1, N_det_generators
k = k+1
write(task(20*(k-1)+1:20*k),'(I9,1X,I9,''|'')') i, N
if (k>=maxtasks) then
k=0
if (add_task_to_taskserver(zmq_to_qp_run_socket,task) == -1) then
stop 'Unable to add task to task server'
endif
endif
end do
if (k > 0) then
if (add_task_to_taskserver(zmq_to_qp_run_socket,task) == -1) then
! /!\ Fragments don't work
! if (i>-ishft(N_det_generators,-2)) then
write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') 0, i, N
ipos += 30
if (ipos > 63970) then
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
stop 'Unable to add task to task server'
endif
ipos=1
endif
! else
! do j=1,fragment_count
! write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, i, N
! ipos += 30
! if (ipos > 63970) then
! if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
! stop 'Unable to add task to task server'
! endif
! ipos=1
! endif
! end do
! endif
enddo
if (ipos > 1) then
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
stop 'Unable to add task to task server'
endif
endif
ASSERT (associated(b%det))
ASSERT (associated(b%val))
@ -75,7 +98,17 @@ subroutine ZMQ_selection(N_in, pt2)
print *, irp_here, ': Failed in zmq_set_running'
endif
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1)
integer :: nproc_target
nproc_target = nproc
double precision :: mem
mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3)
call write_double(6,mem,'Estimated memory/thread (Gb)')
if (qp_max_mem > 0) then
nproc_target = max(1,int(dble(qp_max_mem)/mem))
nproc_target = min(nproc_target,nproc)
endif
!$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc_target+1)
i = omp_get_thread_num()
if (i==0) then
call selection_collector(zmq_socket_pull, b, N, pt2)
@ -131,6 +164,9 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2)
integer, allocatable :: task_id(:)
type(selection_buffer) :: b2
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
call create_selection_buffer(N, N*2, b2)
allocate(task_id(N_det_generators))

1
plugins/Full_CI_ZMQ_GPI2/NEEDED_CHILDREN_MODULES
View File

@ -1 +0,0 @@
Full_CI_ZMQ GPI2

15
plugins/Full_CI_ZMQ_GPI2/README.rst
View File

@ -1,15 +0,0 @@
================
Full_CI_ZMQ_GPI2
================
GPI2 Slave for Full_CI with ZMQ. There should be one instance of the slave
per compute node.
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

105
plugins/Full_CI_ZMQ_GPI2/selection_davidson_slave_gpi2.irp.f
View File

@ -1,105 +0,0 @@
program selection_slave
implicit none
BEGIN_DOC
! Helper program to compute the PT2 in distributed mode.
END_DOC
read_wf = .False.
distributed_davidson = .False.
SOFT_TOUCH read_wf distributed_davidson
call provide_everything
call switch_qp_run_to_master
call run_wf
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count GASPI_is_Initialized
end
subroutine run_wf
use f77_zmq
implicit none
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states)
character*(64) :: states(4)
integer :: rc, i, ierr
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'selection'
states(2) = 'davidson'
states(3) = 'pt2'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,3)
if(trim(zmq_state) == 'Stopped') then
exit
else if (trim(zmq_state) == 'selection') then
! Selection
! ---------
print *, 'Selection'
if (is_gaspi_master) then
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
endif
call broadcast_wf(energy)
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_selection_slave(0,i,energy)
!$OMP END PARALLEL
print *, 'Selection done'
else if (trim(zmq_state) == 'davidson') then
! Davidson
! --------
print *, 'Davidson'
if (is_gaspi_master) then
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
endif
call broadcast_wf(energy)
call omp_set_nested(.True.)
call davidson_slave_tcp(0)
call omp_set_nested(.False.)
print *, 'Davidson done'
else if (trim(zmq_state) == 'pt2') then
! PT2
! ---
print *, 'PT2'
if (is_gaspi_master) then
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
endif
call broadcast_wf(energy)
logical :: lstop
lstop = .False.
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call run_pt2_slave(0,i,energy,lstop)
!$OMP END PARALLEL
print *, 'PT2 done'
endif
end do
end

254
plugins/GPI2/broadcast.irp.f
View File

@ -1,254 +0,0 @@
subroutine broadcast_wf(energy)
implicit none
BEGIN_DOC
! Segment corresponding to the wave function. This is segment 0.
END_DOC
use bitmasks
use GASPI
use ISO_C_BINDING
double precision, intent(inout) :: energy(N_states)
integer(gaspi_return_t) :: res
if (is_gaspi_master) then
call broadcast_wf_put(energy)
else
call broadcast_wf_get(energy)
endif
res = gaspi_barrier(GASPI_GROUP_ALL, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_barrier failed"
stop -1
end if
integer(gaspi_segment_id_t) :: seg_id
do seg_id=0,3
res = gaspi_segment_delete(seg_id)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_segment_delete failed", seg_id
stop -1
end if
end do
end
subroutine broadcast_wf_put(energy)
implicit none
BEGIN_DOC
! Initiates the broadcast of the wave function
END_DOC
use bitmasks
use GASPI
use ISO_C_BINDING
double precision, intent(in) :: energy(N_states)
integer(gaspi_segment_id_t) :: seg_id
integer(gaspi_alloc_t) :: seg_alloc_policy
integer(gaspi_size_t) :: seg_size(0:3)
type(c_ptr) :: seg_ptr(0:3)
integer, pointer :: params_int(:) ! Segment 0
double precision, pointer :: psi_coef_tmp(:,:) ! Segment 1
integer(bit_kind), pointer :: psi_det_tmp(:,:,:) ! Segment 2
double precision, pointer :: params_double(:) ! Segment 3
integer(gaspi_return_t) :: res
seg_alloc_policy = GASPI_MEM_UNINITIALIZED
seg_size(0) = 4 * 5
seg_id=0
res = gaspi_segment_create(seg_id, seg_size(seg_id), GASPI_GROUP_ALL, &
GASPI_BLOCK, seg_alloc_policy)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_create_segment failed", gaspi_rank, seg_id
stop -1
end if
res = gaspi_segment_ptr(seg_id, seg_ptr(seg_id))
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_segment_ptr failed", gaspi_rank
stop -1
end if
call c_f_pointer(seg_ptr(0), params_int, shape=(/ 5 /))
params_int(1) = N_states
params_int(2) = N_det
params_int(3) = psi_det_size
res = gaspi_barrier(GASPI_GROUP_ALL, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_barrier failed", gaspi_rank
stop -1
end if
seg_size(1) = 8 * psi_det_size * N_states
seg_size(2) = bit_kind * psi_det_size * 2 * N_int
seg_size(3) = 8 * N_states
do seg_id=1, 3
res = gaspi_segment_create(seg_id, seg_size(seg_id), GASPI_GROUP_ALL, &
GASPI_BLOCK, seg_alloc_policy)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_create_segment failed", gaspi_rank, seg_id
stop -1
end if
res = gaspi_segment_ptr(seg_id, seg_ptr(seg_id))
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_segment_ptr failed", gaspi_rank
stop -1
end if
end do
call c_f_pointer(seg_ptr(1), psi_coef_tmp, shape=shape(psi_coef))
call c_f_pointer(seg_ptr(2), psi_det_tmp, shape=shape(psi_det))
call c_f_pointer(seg_ptr(3), params_double, shape=(/ N_states /))
psi_coef_tmp = psi_coef
psi_det_tmp = psi_det
params_double = energy
res = gaspi_barrier(GASPI_GROUP_ALL, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_barrier failed", gaspi_rank
stop -1
end if
end
subroutine broadcast_wf_get(energy)
implicit none
BEGIN_DOC
! Gets the broadcasted wave function
END_DOC
use bitmasks
use GASPI
use ISO_C_BINDING
double precision, intent(out) :: energy(N_states)
integer(gaspi_segment_id_t) :: seg_id
integer(gaspi_alloc_t) :: seg_alloc_policy
integer(gaspi_size_t) :: seg_size(0:3)
type(c_ptr) :: seg_ptr(0:3)
integer, pointer :: params_int(:) ! Segment 0
double precision, pointer :: psi_coef_tmp(:,:) ! Segment 1
integer(bit_kind), pointer :: psi_det_tmp(:,:,:) ! Segment 2
double precision, pointer :: params_double(:) ! Segment 3
integer(gaspi_return_t) :: res
seg_alloc_policy = GASPI_MEM_UNINITIALIZED
seg_size(0) = 4 * 5
seg_id=0
res = gaspi_segment_create(seg_id, seg_size(seg_id), GASPI_GROUP_ALL,&
GASPI_BLOCK, seg_alloc_policy)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_create_segment failed"
stop -1
end if
res = gaspi_segment_ptr(seg_id, seg_ptr(seg_id))
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_segment_ptr failed"
stop -1
end if
res = gaspi_barrier(GASPI_GROUP_ALL, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_barrier failed"
stop -1
end if
integer(gaspi_offset_t) :: localOff, remoteOff
integer(gaspi_rank_t) :: remoteRank
integer(gaspi_queue_id_t) :: queue
localOff = 0
remoteRank = 0
queue = 0
res = gaspi_read(seg_id, localOff, remoteRank, &
seg_id, remoteOff, seg_size(seg_id), queue, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_read failed"
stop -1
end if
res = gaspi_wait(queue, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_wait failed"
stop -1
end if
call c_f_pointer(seg_ptr(0), params_int, shape=shape( (/ 5 /) ))
N_states = params_int(1)
N_det = params_int(2)
psi_det_size = params_int(3)
TOUCH N_states N_det psi_det_size
seg_size(1) = 8 * psi_det_size * N_states
seg_size(2) = bit_kind * psi_det_size * 2 * N_int
seg_size(3) = 8 * N_states
do seg_id=1, 3
res = gaspi_segment_create(seg_id, seg_size(seg_id), GASPI_GROUP_ALL, &
GASPI_BLOCK, seg_alloc_policy)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_create_segment failed"
stop -1
end if
res = gaspi_segment_ptr(seg_id, seg_ptr(seg_id))
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_segment_ptr failed"
stop -1
end if
end do
res = gaspi_barrier(GASPI_GROUP_ALL, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_barrier failed"
stop -1
end if
do seg_id=1, 3
res = gaspi_read(seg_id, localOff, remoteRank, &
seg_id, remoteOff, seg_size(seg_id), queue, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_read failed"
stop -1
end if
res = gaspi_wait(queue, GASPI_BLOCK)
if(res .ne. GASPI_SUCCESS) then
write(*,*) "gaspi_wait failed"
stop -1
end if
end do
call c_f_pointer(seg_ptr(1), psi_coef_tmp, shape=shape(psi_coef))
call c_f_pointer(seg_ptr(2), psi_det_tmp, shape=shape(psi_det))
call c_f_pointer(seg_ptr(3), params_double, shape=shape(energy))
psi_coef = psi_coef_tmp
psi_det = psi_det_tmp
energy = params_double
end

2
plugins/Generators_CAS/generators.irp.f
View File

@ -35,10 +35,10 @@ END_PROVIDER
integer, external :: number_of_holes,number_of_particles
integer, allocatable :: nongen(:)
integer :: inongen
inongen = 0
allocate(nongen(N_det))
inongen = 0
m=0
do i=1,N_det
good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 )

6
plugins/Generators_full/generators.irp.f
View File

@ -29,7 +29,6 @@ END_PROVIDER
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
integer :: i, k
psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
@ -44,12 +43,9 @@ END_PROVIDER
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
integer :: i, k
psi_det_sorted_gen = psi_det_sorted
psi_coef_sorted_gen = psi_coef_sorted
!do i=1,N_det_generators
psi_det_sorted_gen_order = psi_det_sorted_order
!end do
psi_det_sorted_gen_order = psi_det_sorted_order
END_PROVIDER

2
plugins/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (ao_num) ]
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_tot_num) ]
implicit none
BEGIN_DOC

2
plugins/Hartree_Fock/localize_mos.irp.f
View File

@ -21,7 +21,7 @@ program localize_mos
mo_coef(1,1),size(mo_coef,1),1.d-6,rank)
print *, rank
if (elec_alpha_num>elec_alpha_num) then
if (elec_alpha_num>elec_beta_num) then
W = 0.d0
do k=elec_beta_num+1,elec_alpha_num
do j=1,ao_num

2
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -1,6 +1,6 @@
use bitmasks
BEGIN_PROVIDER [ integer, mrmode ]
BEGIN_PROVIDER [ integer, mrmode ]
mrmode = 0
END_PROVIDER

80
plugins/MRPT_Utils/mrpt_utils.irp.f
View File

@ -16,23 +16,23 @@
integer :: i,j,m
integer :: i_state
double precision :: accu(N_states)
double precision, allocatable :: delta_ij_tmp(:,:,:)
double precision, allocatable :: delta_ij_local(:,:,:)
delta_ij_mrpt = 0.d0
allocate (delta_ij_tmp(N_det,N_det,N_states))
allocate (delta_ij_local(N_det,N_det,N_states))
! 1h
delta_ij_tmp = 0.d0
call H_apply_mrpt_1h(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call H_apply_mrpt_1h(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_1h(i_state) = accu(i_state)
@ -40,14 +40,14 @@
print*, '1h = ',accu
! 1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1p(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call H_apply_mrpt_1p(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_1p(i_state) = accu(i_state)
@ -55,15 +55,15 @@
print*, '1p = ',accu
! 1h1p
delta_ij_tmp = 0.d0
call H_apply_mrpt_1h1p(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call H_apply_mrpt_1h1p(delta_ij_local,N_det)
double precision :: e_corr_from_1h1p_singles(N_states)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
@ -72,14 +72,14 @@
! 1h1p third order
if(do_third_order_1h1p)then
delta_ij_tmp = 0.d0
call give_1h1p_sec_order_singles_contrib(delta_ij_tmp)
delta_ij_local = 0.d0
call give_1h1p_sec_order_singles_contrib(delta_ij_local)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_1h1p(i_state) = accu(i_state)
@ -88,14 +88,14 @@
endif
! 2h
delta_ij_tmp = 0.d0
call H_apply_mrpt_2h(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call H_apply_mrpt_2h(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_2h(i_state) = accu(i_state)
@ -103,14 +103,14 @@
print*, '2h = ',accu
! 2p
delta_ij_tmp = 0.d0
call H_apply_mrpt_2p(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call H_apply_mrpt_2p(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_2p(i_state) = accu(i_state)
@ -118,15 +118,15 @@
print*, '2p = ',accu
! 1h2p
delta_ij_tmp = 0.d0
call give_1h2p_contrib(delta_ij_tmp)
call H_apply_mrpt_1h2p(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call give_1h2p_contrib(delta_ij_local)
call H_apply_mrpt_1h2p(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_1h2p(i_state) = accu(i_state)
@ -134,15 +134,15 @@
print*, '1h2p = ',accu
! 2h1p
delta_ij_tmp = 0.d0
call give_2h1p_contrib(delta_ij_tmp)
call H_apply_mrpt_2h1p(delta_ij_tmp,N_det)
delta_ij_local = 0.d0
call give_2h1p_contrib(delta_ij_local)
call H_apply_mrpt_2h1p(delta_ij_local,N_det)
accu = 0.d0
do i_state = 1, N_states
do i = 1, N_det
do j = 1, N_det
accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
enddo
enddo
second_order_pt_new_2h1p(i_state) = accu(i_state)
@ -150,14 +150,14 @@
print*, '2h1p = ',accu
! 2h2p
!delta_ij_tmp = 0.d0
!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det)
!delta_ij_local = 0.d0
!call H_apply_mrpt_2h2p(delta_ij_local,N_det)
!accu = 0.d0
!do i_state = 1, N_states
!do i = 1, N_det
! do j = 1, N_det
! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij_mrpt(j,i,i_state) += delta_ij_tmp(j,i,i_state)
! accu(i_state) += delta_ij_local(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state)
! delta_ij_mrpt(j,i,i_state) += delta_ij_local(j,i,i_state)
! enddo
!enddo
!second_order_pt_new_2h2p(i_state) = accu(i_state)

4
plugins/Psiref_Utils/psi_ref_utils.irp.f
View File

@ -9,7 +9,7 @@ use bitmasks
! function.
END_DOC
call sort_dets_by_det_search_key(N_det_ref, psi_ref, psi_ref_coef, &
psi_ref_sorted_bit, psi_ref_coef_sorted_bit)
psi_ref_sorted_bit, psi_ref_coef_sorted_bit, N_states)
END_PROVIDER
@ -152,7 +152,7 @@ END_PROVIDER
! function.
END_DOC
call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, &
psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit)
psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit, N_states)
END_PROVIDER

4
plugins/QMC/qp_convert_qmcpack_to_ezfio.py
View File

@ -60,7 +60,7 @@ beta = ezfio.get_electrons_elec_beta_num()
print "elec_alpha_num", alpha
print "elec_beta_num", beta
print "elec_tot_num", alpha + beta
print "spin_multiplicity", 2 * (alpha - beta) + 1
print "spin_multiplicity", (alpha - beta) + 1
l_label = ezfio.get_nuclei_nucl_label()
l_charge = ezfio.get_nuclei_nucl_charge()
@ -133,7 +133,7 @@ d_gms_order ={ 0:["s"],
1:[ "x", "y", "z" ],
2:[ "xx", "yy", "zz", "xy", "xz", "yz" ],
3:[ "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz"],
4: ["xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy", "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz","zzxy"] }
4:[ "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy"] }
def compare_gamess_style(item1, item2):
n1,n2 = map(len,(item1,item2))

83
plugins/QMC/truncate_wf_spin.irp.f
View File

@ -65,53 +65,52 @@ subroutine run
enddo
!$OMP END PARALLEL DO
endif
if (ci_threshold > norm_sort(j)) then
cycle
if (ci_threshold <= norm_sort(j)) then
exit
endif
u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_0(1:N_det,1:N_states) = 0.d0
u_t(1:N_states,1:N_det) = 0.d0
v_t(1:N_states,1:N_det) = 0.d0
s_t(1:N_states,1:N_det) = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
print *, 'Computing H|Psi> ...'
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
print *, 'Done'
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
print *, 'E = ', e_0(i) + nuclear_repulsion
enddo
m = 0
do k=1,n_det
if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
m = m+1
endif
enddo
do k=1,N_states
E = E_0(k) + nuclear_repulsion
enddo
print *, 'Number of determinants:', m
exit
enddo
m = 0
do k=1,n_det
if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
m = m+1
endif
enddo
do k=1,N_states
E = E_0(k) + nuclear_repulsion
enddo
print *, 'Number of determinants:', m
call wf_of_psi_bilinear_matrix(.True.)
call save_wavefunction
u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_0(1:N_det,1:N_states) = 0.d0
u_t(1:N_states,1:N_det) = 0.d0
v_t(1:N_states,1:N_det) = 0.d0
s_t(1:N_states,1:N_det) = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
print *, 'Computing H|Psi> ...'
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
print *, 'Done'
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
print *, 'E(',i,') = ', e_0(i) + nuclear_repulsion
enddo
deallocate (iorder, norm_sort)
end

28
plugins/Symmetry/find_sym.irp.f
View File

@ -6,34 +6,6 @@ BEGIN_PROVIDER [ logical, molecule_is_linear ]
molecule_is_linear = (minval(inertia_tensor_eigenvalues) < 1.d-5)
END_PROVIDER
BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ]
implicit none
BEGIN_DOC
! If true, there is a center of inversion in the WF
END_DOC
molecule_has_center_of_inversion = .True.
integer :: i,j,k
double precision :: point(3)
logical :: found
double precision, external :: u_dot_u
do i=1,nucl_num
found = .False.
do j=1,nucl_num
if (nucl_charge(i) /= nucl_charge(j)) cycle
point(:) = nucl_coord_sym_transp(:,i) + nucl_coord_sym_transp(:,j)
if (u_dot_u(point,3) < 1.d-5) then
found = .True.
exit
endif
enddo
if (.not.found) then
molecule_has_center_of_inversion = .False.
exit
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, sym_rotation_axis, (3) ]
implicit none

17
plugins/dress_zmq/EZFIO.cfg
View File

@ -1,17 +0,0 @@
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[thresh_dressed_ci]
type: Threshold
doc: Threshold on the convergence of the dressed CI energy
interface: ezfio,provider,ocaml
default: 1.e-5
[n_it_max_dressed_ci]
type: Strictly_positive_int
doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 10

2
plugins/dress_zmq/alpha_factory.irp.f
View File

@ -637,7 +637,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, indexes, ab
integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
logical, intent(inout) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num)
integer, intent(inout) :: abuf(*)
integer, intent(inout) :: abuf(*)
integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt, s
integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
integer :: phasemask(2,N_int*bit_kind_size)

33
plugins/dress_zmq/dress_general.irp.f
View File

@ -28,36 +28,41 @@ subroutine run_dressing(N_st,energy)
E_new = 0.d0
delta_E = 1.d0
iteration = 0
do while (delta_E > thresh_dress)
N_det_delta_ij = N_det
touch N_det_delta_ij
iteration += 1
do iteration=1,n_it_dress_max
print *, '==============================================='
print *, 'Iteration', iteration, '/', n_it_dress_max
print *, '==============================================='
print *, ''
E_old = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
!print *, "DELTA IJ", delta_ij(1,1,1)
!if(.true.) provide delta_ij_tmp
if(.true.) call delta_ij_done()
E_old = sum(psi_energy(:))
print *, 'Variational energy <Psi|H|Psi>'
do i=1,N_st
if(.true.) call write_double(6,ci_energy_dressed(i),"Energy")
print *, i, psi_energy(i)+nuclear_repulsion
enddo
print *, 'Dressed energy <Psi|H+Delta|Psi>'
do i=1,N_st
print *, i, ci_energy_dressed(i)
enddo
call diagonalize_ci_dressed
E_new = dress_e0_denominator(1) !sum(ci_energy_dressed(1:N_states))
E_new = sum(psi_energy(:))
delta_E = (E_new - E_old)/dble(N_states)
print *, ''
call write_double(6,thresh_dress,"thresh_dress")
call write_double(6,delta_E,"delta_E")
call write_double(6,delta_E,"delta_E (undressed)")
delta_E = dabs(delta_E)
call save_wavefunction
! call ezfio_set_dress_zmq_energy(ci_energy_dressed(1))
if (iteration >= n_it_dress_max) then
if (delta_E < thresh_dress) then
exit
endif
enddo
if(.true.) call write_double(6,ci_energy_dressed(1),"Final energy")
print *, 'Variational energy <Psi|H|Psi>'
do i=1,N_st
print *, i, psi_energy(i)+nuclear_repulsion
enddo
print *, 'Dressed energy <Psi|H+Delta|Psi>'
do i=1,N_st
print *, i, ci_energy_dressed(i)+nuclear_repulsion
enddo
endif
if(.true.) energy(1:N_st) = 0d0 ! ci_energy_dressed(1:N_st)

10
plugins/dress_zmq/dress_slave.irp.f
View File

@ -28,7 +28,7 @@ subroutine run_wf
double precision :: energy(N_states_diag)
character*(64) :: states(1)
integer :: rc, i
integer, external :: zmq_get_dvector, zmq_get_N_det_generators
integer, external :: zmq_get_dvector, zmq_get_N_det_generators
integer, external :: zmq_get_psi, zmq_get_N_det_selectors
integer, external :: zmq_get_N_states_diag
double precision :: tmp
@ -50,19 +50,15 @@ integer, external :: zmq_get_dvector, zmq_get_N_det_generators
else if (zmq_state(:5) == 'dress') then
! Dress
! ---------
!call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
!TOUCH psi_det
if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'dress_stoch_istate',tmp,1) == -1) cycle
dress_stoch_istate = int(tmp)
TOUCH dress_stoch_istate
TOUCH state_average_weight
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy dress_stoch_istate state_average_weight
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_gen_order psi_bilinear_matrix_order

58
plugins/dress_zmq/dress_stoch_routines.irp.f
View File

@ -274,16 +274,27 @@ subroutine ZMQ_dress(E, dress, delta_out, delta_s2_out, relative_error)
print *, irp_here, ': Failed in zmq_set_running'
endif
!!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc) &
! !$OMP PRIVATE(i)
!i = omp_get_thread_num()
!if (i==0) then
call omp_set_nested(.true.)
if (.true.) then !! TODO
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) &
!$OMP PRIVATE(i)
i = omp_get_thread_num()
if (i==0) then
call dress_collector(zmq_socket_pull,E, relative_error, delta, delta_s2, dress,&
dress_stoch_istate)
!else
! call dress_slave_inproc(i)
!endif
!!$OMP END PARALLEL
else
call dress_slave_inproc(i)
endif
!$OMP END PARALLEL
else
call dress_collector(zmq_socket_pull,E, relative_error, delta, delta_s2, dress,&
dress_stoch_istate)
endif
call omp_set_nested(.false.)
delta_out(dress_stoch_istate,1:N_det) = delta(dress_stoch_istate,1:N_det)
delta_s2_out(dress_stoch_istate,1:N_det) = delta_s2(dress_stoch_istate,1:N_det)
@ -451,13 +462,18 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
end do
t = dress_dot_t(m)
avg = S(t) / dble(c)
eqt = (S2(t) / c) - (S(t)/c)**2
eqt = sqrt(eqt / dble(c-1))
error = eqt
time = omp_get_wtime()
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', c, avg+E0+E(dress_stoch_istate), eqt, time-time0, ''
if (c > 1) then
eqt = (S2(t) / c) - (S(t)/c)**2
eqt = sqrt(eqt / dble(c-1))
error = eqt
time = omp_get_wtime()
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', c, avg+E0+E(dress_stoch_istate), eqt, time-time0, ''
else
eqt = 1.d0
error = eqt
endif
m += 1
if(eqt <= relative_error) then
if(eqt <=0d0* relative_error) then
integer, external :: zmq_put_dvector
i= zmq_put_dvector(zmq_to_qp_run_socket, worker_id, "ending", dble(m-1), 1)
found = .true.
@ -480,6 +496,7 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
end if
end do
do i=1,n_tasks
if(edI(1, edI_index(i)) /= 0d0) stop "NIN M"
edI(:, edI_index(i)) += edI_task(:, i)
end do
dot_f(m_task) -= f
@ -515,6 +532,16 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
end do
dress(istate) = E(istate)+E0+avg
if(ff /= 0) stop "WRONG NUMBER OF FRAGMENTS COLLECTED"
!double precision :: tmp
!tmp = 0d0
!do i=1,N_det
! if(edi(1,i) == 0d0) stop "EMPTY"
! tmp += psi_coef(i, 1) * delta(1, i)
!end do
!print *, "SUM", E(1)+sum(edi(1,:))
!print *, "DOT", E(1)+tmp
call disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id)
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
end subroutine
@ -601,6 +628,3 @@ end function
END_PROVIDER

2
plugins/dress_zmq/dress_zmq_routines.irp.f
View File

@ -6,7 +6,7 @@ subroutine dress_zmq()
threshold_generators = 1d0
read_wf = .True.
SOFT_TOUCH read_wf
SOFT_TOUCH read_wf threshold_generators threshold_selectors
if (.True.) then
integer :: i,j

16
plugins/dress_zmq/dressing.irp.f
View File

@ -65,12 +65,14 @@ END_PROVIDER
BEGIN_PROVIDER [ integer , N_det_delta_ij ]
implicit none
!N_det_delta_ij = 0!N_det
N_det_delta_ij = N_det
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_states, N_det, 2) ]
implicit none
if(.true.) delta_ij(:,:N_det_delta_ij, :) = delta_ij_tmp(:,:,:)
if(.true.) then
delta_ij(:,:N_det_delta_ij, :) = delta_ij_tmp(:,:,:)
endif
delta_ij(:,N_det_delta_ij+1:,:) = 0d0
END_PROVIDER
@ -84,15 +86,6 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ]
double precision :: E_CI_before(N_states), relative_error
integer :: cnt = 0
! prevents re-providing if delta_ij_tmp is
! just being copied
!if(N_det_delta_ij /= N_det) return
!cnt += 1
!if(mod(cnt,2) == 0) return
if(.true.) then
allocate(dress(N_states), del(N_states, N_det_delta_ij), del_s2(N_states, N_det_delta_ij))
delta_ij_tmp = 0d0
@ -108,7 +101,6 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ]
deallocate(dress, del, del_s2)
end if
END_PROVIDER

7
plugins/dress_zmq/dressing_vector.irp.f
View File

@ -18,14 +18,7 @@
do j = 1, n_det
dressing_column_h(j,k) = delta_ij(k,j,1)
dressing_column_s(j,k) = delta_ij(k,j,2)
! print *, j, delta_ij(k,j,:)
enddo
! tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det) &
! - dressing_column_h(l,k) * psi_coef(l,k)
! dressing_column_h(l,k) -= tmp * f
! tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det) &
! - dressing_column_s(l,k) * psi_coef(l,k)
! dressing_column_s(l,k) -= tmp * f
enddo
END_PROVIDER

19
plugins/dress_zmq/energy.irp.f
View File

@ -13,13 +13,24 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ]
END_DOC
integer :: i
if (initialize_dress_E0_denominator) then
call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
do i=N_det+1,N_states
dress_E0_denominator(i) = 0.d0
enddo
if (h0_type == "EN") then
dress_E0_denominator(1:N_states) = psi_energy(1:N_states)
else if (h0_type == "Barycentric") then
! dress_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
dress_E0_denominator(1:N_states) = minval(diagonal_H_matrix_on_psi_det(1:N_det))
else
print *, h0_type, ' not implemented'
stop
endif
! call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
! do i=N_det+1,N_states
! dress_E0_denominator(i) = 0.d0
! enddo
call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator')
else
dress_E0_denominator = -huge(1.d0)
endif
END_PROVIDER

58
plugins/dress_zmq/extra_functions.irp.f.example Normal file
View File

@ -0,0 +1,58 @@
BEGIN_PROVIDER [ integer, N_dress_int_buffer ]
&BEGIN_PROVIDER [ integer, N_dress_double_buffer ]
&BEGIN_PROVIDER [ integer, N_dress_det_buffer ]
implicit none
N_dress_int_buffer = 1
N_dress_double_buffer = 1
N_dress_det_buffer = 1
END_PROVIDER
subroutine delta_ij_done()
BEGIN_DOC
! This subroutine is executed on the master when the dressing has been computed,
! before the diagonalization.
END_DOC
end
subroutine dress_pulled(ind, int_buf, double_buf, det_buf, N_buf)
use bitmasks
implicit none
BEGIN_DOC
! Dress the contributions pulled from the slave.
END_DOC
integer, intent(in) :: ind, N_buf(3)
integer, intent(in) :: int_buf(*)
double precision, intent(in) :: double_buf(*)
integer(bit_kind), intent(in) :: det_buf(N_int,2,*)
end
subroutine generator_start(i_gen, iproc)
implicit none
BEGIN_DOC
! This subroutine is executed on the slave before computing the contribution of a generator.
END_DOC
integer, intent(in) :: i_gen, iproc
integer :: i
end
subroutine generator_done(i_gen, int_buf, double_buf, det_buf, N_buf, iproc)
implicit none
BEGIN_DOC
! This subroutine is executed on the slave after computing the contribution of a generator.
END_DOC
integer, intent(in) :: i_gen, iproc
integer, intent(out) :: int_buf(N_dress_int_buffer), N_buf(3)
double precision, intent(out) :: double_buf(N_dress_double_buffer)
integer(bit_kind), intent(out) :: det_buf(N_int, 2, N_dress_det_buffer)
N_buf(:) = 1
int_buf(:) = 0
double_buf(:) = 0.d0
det_buf(:,:,:) = 0
end

5
plugins/dress_zmq/run_dress_slave.irp.f
View File

@ -41,7 +41,7 @@ subroutine run_dress_slave(thread,iproce,energy)
! double precision, external :: omp_get_wtime
double precision :: time, time0
integer :: ntask_tbd, task_tbd(Nproc), i_gen_tbd(Nproc), subset_tbd(Nproc)
if(iproce /= 0) stop "RUN DRESS SLAVE is OMP"
! if(iproce /= 0) stop "RUN DRESS SLAVE is OMP"
allocate(delta_det(N_states, N_det, 0:pt2_N_teeth+1, 2))
allocate(cp(N_states, N_det, dress_N_cp, 2))
@ -67,9 +67,8 @@ double precision :: time, time0
will_send = 0
double precision :: hij, sij, tmp
logical :: purge
purge_task_id = 0
hij = E0_denominator(1) !PROVIDE BEFORE OMP PARALLEL
provide psi_energy
ending(1) = dble(dress_N_cp+1)
ntask_tbd = 0
!$OMP PARALLEL DEFAULT(SHARED) &

1
plugins/mrcc_sto/NEEDED_CHILDREN_MODULES
View File

@ -1 +0,0 @@
dress_zmq DavidsonDressed Psiref_CAS MRPT_Utils Perturbation MRCC_Utils

12
plugins/mrcc_sto/README.rst
View File

@ -1,12 +0,0 @@
========
mrcc_sto
========
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.

240
plugins/mrcc_sto/mrcc_sto.irp.f
View File

@ -1,240 +0,0 @@
program mrcc_sto
implicit none
BEGIN_DOC
! TODO
END_DOC
call dress_zmq()
call ezfio_set_mrcc_sto_energy(ci_energy_dressed(1))
end
BEGIN_PROVIDER [ double precision, hij_cache_, (N_det,Nproc) ]
&BEGIN_PROVIDER [ double precision, sij_cache_, (N_det,Nproc) ]
&BEGIN_PROVIDER [ double precision, dIa_hla_, (N_states,N_det,Nproc) ]
&BEGIN_PROVIDER [ double precision, dIa_sla_, (N_states,N_det,Nproc) ]
&BEGIN_PROVIDER [ integer, excs_ , (0:2,2,2,N_det,Nproc) ]
&BEGIN_PROVIDER [ double precision, phases_, (N_det, Nproc) ]
BEGIN_DOC
! temporay arrays for dress_with_alpha_buffer. Avoids reallocation.
END_DOC
END_PROVIDER
subroutine dress_with_alpha_buffer(delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc)
use bitmasks
implicit none
BEGIN_DOC
!delta_ij_loc(:,:,1) : dressing column for H
!delta_ij_loc(:,:,2) : dressing column for S2
!minilist : indices of determinants connected to alpha ( in psi_det_sorted )
!n_minilist : size of minilist
!alpha : alpha determinant
END_DOC
integer(bit_kind), intent(in) :: alpha(N_int,2), det_minilist(N_int, 2, n_minilist)
integer,intent(in) :: minilist(n_minilist), n_minilist, iproc, i_gen
double precision, intent(inout) :: delta_ij_loc(N_states,N_det,2)
integer :: i,j,k,l,m
integer :: degree1, degree2, degree
double precision :: hIk, hla, hIl, sla, dIk(N_states), dka(N_states), dIa(N_states), hka
double precision :: phase, phase2
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(N_int,2), ctrl
integer :: i_state, k_sd, l_sd, m_sd, ll_sd, i_I
double precision :: Delta_E_inv(N_states)
double precision :: sdress, hdress
logical :: ok, ok2
integer :: canbediamond
PROVIDE mo_class
if(n_minilist == 1) return
do i=1,n_minilist
if(idx_non_ref_rev(minilist(i)) == 0) return
end do
if (perturbative_triples) then
PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
endif
canbediamond = 0
do l_sd=1,n_minilist
call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int)
call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. &
(mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V')
if(ok .and. degree1 == 2) then
ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. &
(mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V')
end if
if(ok) then
canbediamond += 1
excs_(:,:,:,l_sd,iproc) = exc(:,:,:)
phases_(l_sd, iproc) = phase
else
phases_(l_sd, iproc) = 0d0
end if
!call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc))
!call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc))
call i_h_j_s2(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc), sij_cache_(l_sd,iproc))
enddo
if(canbediamond <= 1) return
do i_I=1,N_det_ref
call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,N_int)
if (degree1 > 4) then
cycle
endif
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
do k_sd=1,n_minilist
if(phases_(k_sd,iproc) == 0d0) cycle
call get_excitation_degree(psi_ref(1,1,i_I),det_minilist(1,1,k_sd),degree,N_int)
if (degree > 2) then
cycle
endif
!call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int)
phase = phases_(k_sd, iproc)
exc(:,:,:) = excs_(:,:,:,k_sd,iproc)
degree2 = exc(0,1,1) + exc(0,1,2)
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int)
if((.not. ok) .and. (.not. perturbative_triples)) cycle
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
ok2 = .false.
!do i_state=1,N_states
! !if(dka(i_state) == 0) cycle
! dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state)
! if(dIk(i_state) /= 0d0) then
! ok2 = .true.
! endif
!enddo
!if(.not. ok2) cycle
if (ok) then
phase2 = 0d0
do l_sd=k_sd+1,n_minilist
if(phases_(l_sd, iproc) == 0d0) cycle
call get_excitation_degree(tmp_det,det_minilist(1,1,l_sd),degree,N_int)
if (degree == 0) then
do i_state=1,N_states
dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state)
if(dIk(i_state) /= 0d0) then
if(phase2 == 0d0) call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,l_sd),exc,degree,phase2,N_int)
dka(i_state) = dij(i_I, idx_non_ref_rev(minilist(l_sd)), i_state) * phase * phase2
end if
end do
!call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,l_sd),exc,degree,phase2,N_int)
!do i_state=1,N_states
! if(dIk(i_state) /= 0d0) dka(i_state) = dij(i_I, idx_non_ref_rev(minilist(l_sd)), i_state) * phase * phase2
!enddo
exit
endif
enddo
else if (perturbative_triples) then
hka = hij_cache_(k_sd,iproc)
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif
if (perturbative_triples.and. (degree2 == 1) ) then
call i_h_j(psi_ref(1,1,i_I),tmp_det,N_int,hka)
hka = hij_cache_(k_sd,iproc) - hka
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv)
do i_state=1,N_states
ASSERT (Delta_E_inv(i_state) < 0.d0)
dka(i_state) = hka / Delta_E_inv(i_state)
enddo
endif
endif
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
enddo
ok2 = .false.
do i_state=1,N_states
if(dIa(i_state) /= 0d0) ok2 = .true.
enddo
if(.not. ok2) cycle
do l_sd=1,n_minilist
k_sd = minilist(l_sd)
hla = hij_cache_(l_sd,iproc)
sla = sij_cache_(l_sd,iproc)
do i_state=1,N_states
hdress = dIa(i_state) * hla * psi_ref_coef(i_I,i_state)
sdress = dIa(i_state) * sla * psi_ref_coef(i_I,i_state)
!!!$OMP ATOMIC
delta_ij_loc(i_state,k_sd,1) += hdress
!!!$OMP ATOMIC
delta_ij_loc(i_state,k_sd,2) += sdress
enddo
enddo
enddo
end subroutine
!! TESTS MINILIST
subroutine test_minilist(minilist, n_minilist, alpha)
use bitmasks
implicit none
integer, intent(in) :: n_minilist
integer(bit_kind),intent(in) :: alpha(N_int, 2)
integer, intent(in) :: minilist(n_minilist)
integer :: a, i, deg
integer :: refc(N_det), testc(N_det)
refc = 0
testc = 0
do i=1,N_det
call get_excitation_degree(psi_det(1,1,i), alpha, deg, N_int)
if(deg <= 2) refc(i) = refc(i) + 1
end do
do i=1,n_minilist
call get_excitation_degree(psi_det(1,1,minilist(i)), alpha, deg, N_int)
if(deg <= 2) then
testc(minilist(i)) += 1
else
stop "NON LINKED IN MINILIST"
end if
end do
do i=1,N_det
if(refc(i) /= testc(i)) then
print *, "MINILIST FAIL ", sum(refc), sum(testc), n_minilist
exit
end if
end do
end subroutine

4
plugins/mrcepa0/EZFIO.cfg
View File

@ -18,7 +18,7 @@ interface: ezfio
type: logical
doc: Compute perturbative contribution of the Triples
interface: ezfio,provider,ocaml
default: false
default: true
[energy]
type: double precision
@ -35,5 +35,5 @@ default: 1.e-5
type: Strictly_positive_int
doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 10
default: 30

318
plugins/mrcepa0/dressing.irp.f
View File

@ -74,118 +74,6 @@ BEGIN_PROVIDER [ double precision, mrcc_norm_acc, (0:N_det_non_ref, N_states) ]
END_PROVIDER
! BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref) ]
!&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref) ]
! use bitmasks
! implicit none
! integer :: gen, h, p, n, t, i, j, h1, h2, p1, p2, s1, s2, iproc
! integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2)
! integer(bit_kind),allocatable :: buf(:,:,:)
! logical :: ok
! logical, external :: detEq
! integer, external :: omp_get_thread_num
! double precision :: coefs(N_det_non_ref), myCoef
! integer :: n_in_teeth
! double precision :: contrib(N_states), curn, in_teeth_step, curlim, curnorm
!
! contrib = 0d0
! read(*,*) n_in_teeth
! !n_in_teeth = 2
! in_teeth_step = 1d0 / dfloat(n_in_teeth)
! !double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref)
! !double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref)
!
! coefs = 0d0
! coefs(:mrcc_teeth(1,1)-1) = 1d0
!
! do i=1,N_mrcc_teeth
! print *, "TEETH SIZE", i, mrcc_teeth(i+1,1)-mrcc_teeth(i,1)
! if(mrcc_teeth(i+1,1) - mrcc_teeth(i,1) <= n_in_teeth) then
! coefs(mrcc_teeth(i,1):mrcc_teeth(i+1,1)-1) = 1d0
! else if(.false.) then
! curnorm = 0d0
! curn = 0.5d0
! curlim = curn / dfloat(n_in_teeth)
! do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
! if(mrcc_norm_acc(j,1) >= curlim) then
! coefs(j) = 1d0
! curnorm += mrcc_norm(j,1)
! do while(mrcc_norm_acc(j,1) > curlim)
! curn += 1d0
! curlim = curn / dfloat(n_in_teeth)
! end do
! end if
! end do
! do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
! coefs(j) = coefs(j) / curnorm ! 1d0 / norm computed in teeth
! end do
! else if(.true.) then
! coefs(mrcc_teeth(i,1):mrcc_teeth(i,1)+n_in_teeth-1) = 1d0 / mrcc_norm_acc(mrcc_teeth(i,1)+n_in_teeth-1, 1)
! else
! curnorm = 0d0
! n = mrcc_teeth(i+1,1) - mrcc_teeth(i,1)
! do j=1,n_in_teeth
! t = int((dfloat(j)-0.5d0) * dfloat(n) / dfloat(n_in_teeth)) + 1 + mrcc_teeth(i,1) - 1
! curnorm += mrcc_norm(t,1)
! coefs(t) = 1d0
! end do
! do j=mrcc_teeth(i,1), mrcc_teeth(i+1,1)-1
! coefs(j) = coefs(j) / curnorm ! 1d0 / norm computed in teeth
! end do
! end if
! !coefs(mrcc_teeth(i,1)) =
! end do
!
! !coefs = coefs * dfloat(N_det_generators)
!
!
! delta_ij_mrcc_sto = 0d0
! delta_ij_s2_mrcc_sto = 0d0
! PROVIDE dij
! provide hh_shortcut psi_det_size! lambda_mrcc
! !$OMP PARALLEL DO default(none) schedule(dynamic) &
! !$OMP shared(psi_ref, psi_non_ref, hh_exists, pp_exists, N_int, hh_shortcut) &
! !$OMP shared(N_det_generators, coefs,N_det_non_ref, delta_ij_mrcc_sto) &
! !$OMP shared(contrib,psi_det_generators, delta_ij_s2_mrcc_sto) &
! !$OMP private(i,j,curnorm,myCoef, h, n, mask, omask, buf, ok, iproc)
! do gen= 1,N_det_generators
! if(coefs(gen) == 0d0) cycle
! myCoef = coefs(gen)
! allocate(buf(N_int, 2, N_det_non_ref))
! iproc = omp_get_thread_num() + 1
! if(mod(gen, 1000) == 0) print *, "mrcc_sto ", gen, "/", N_det_generators
!
! do h=1, hh_shortcut(0)
! call apply_hole_local(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int)
! if(.not. ok) cycle
! omask = 0_bit_kind
! if(hh_exists(1, h) /= 0) omask = mask
! n = 1
! do p=hh_shortcut(h), hh_shortcut(h+1)-1
! call apply_particle_local(mask, pp_exists(1, p), buf(1,1,n), ok, N_int)
! if(ok) n = n + 1
! if(n > N_det_non_ref) stop "Buffer too small in MRCC..."
! end do
! n = n - 1
! if(n /= 0) then
! call mrcc_part_dress(delta_ij_mrcc_sto, delta_ij_s2_mrcc_sto, &
! gen,n,buf,N_int,omask,myCoef,contrib)
! endif
! end do
! deallocate(buf)
! end do
! !$OMP END PARALLEL DO
!
!
!
! curnorm = 0d0
! do j=1,N_det_non_ref
! curnorm += delta_ij_mrcc_sto(1,j)*delta_ij_mrcc_sto(1,j)
! end do
! print *, "NORM DELTA ", dsqrt(curnorm)
!
!END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref) ]
@ -251,7 +139,7 @@ END_PROVIDER
&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref) ]
use bitmasks
implicit none
integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc
integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2
integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2)
integer(bit_kind),allocatable :: buf(:,:,:)
logical :: ok
@ -266,13 +154,15 @@ END_PROVIDER
delta_ij_s2_mrcc = 0d0
!$OMP PARALLEL DO default(none) schedule(dynamic) &
!$OMP PARALLEL default(none) &
!$OMP shared(contrib,psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) &
!$OMP shared(N_det_non_ref, N_det_ref, delta_ij_mrcc, delta_ij_s2_mrcc) &
!$OMP private(h, n, mask, omask, buf, ok, iproc)
!$OMP private(h, n, mask, omask, buf, ok,gen)
allocate(buf(N_int, 2, N_det_non_ref))
!$OMP DO schedule(dynamic)
do gen= 1, N_det_generators
allocate(buf(N_int, 2, N_det_non_ref))
iproc = omp_get_thread_num() + 1
if(mod(gen, 1000) == 0) print *, "mrcc ", gen, "/", N_det_generators
do h=1, hh_shortcut(0)
call apply_hole_local(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int)
@ -292,9 +182,12 @@ END_PROVIDER
endif
end do
deallocate(buf)
end do
!$OMP END PARALLEL DO
!$OMP END DO
deallocate(buf)
!$OMP END PARALLEL
END_PROVIDER
@ -323,11 +216,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
double precision :: hIk, hla, hIl, sla, dIk(N_states), dka(N_states), dIa(N_states), hka
double precision, allocatable :: dIa_hla(:,:), dIa_sla(:,:)
double precision :: haj, phase, phase2
double precision :: f(N_states), ci_inv(N_states)
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(Nint,2)
integer :: iint, ipos
integer :: i_state, k_sd, l_sd, i_I, i_alpha
integer(bit_kind),allocatable :: miniList(:,:,:)
@ -345,6 +236,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
double precision, intent(inout) :: contrib(N_states)
double precision :: sdress, hdress
PROVIDE n_act_orb elec_num
if (perturbative_triples) then
PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
@ -529,9 +421,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
enddo
enddo
do i_state=1,N_states
ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state)
enddo
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
hla = hij_cache(k_sd)
@ -542,13 +431,13 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
enddo
enddo
do i_state=1,N_states
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state)
sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state)
!!$OMP ATOMIC
!$OMP ATOMIC
contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state)
contrib(i_state) += hdress * psi_non_ref_coef(k_sd, i_state)
!$OMP ATOMIC
delta_ij_(i_state,k_sd) += hdress
!$OMP ATOMIC
@ -581,7 +470,7 @@ END_PROVIDER
double precision, allocatable :: mrcc(:)
double precision :: E_CI_before!, relative_error
double precision, save :: target_error = 0d0
double precision :: target_error
allocate(mrcc(N_states))
@ -593,12 +482,7 @@ END_PROVIDER
E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion
threshold_selectors = 1.d0
threshold_generators = 1d0
if(target_error /= 0d0) then
target_error = target_error / 2d0 ! (-mrcc_E0_denominator(1) + mrcc_previous_E(1)) / 1d1
else
target_error = 1d-4
end if
target_error = 0d0
target_error = thresh_dressed_ci * 5.d-2
call ZMQ_mrcc(E_CI_before, mrcc, delta_ij_mrcc_zmq, delta_ij_s2_mrcc_zmq, abs(target_error))
mrcc_previous_E(:) = mrcc_E0_denominator(:)
@ -610,21 +494,7 @@ END_PROVIDER
use bitmasks
implicit none
integer :: i, j, i_state
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc, 4=stoch
! if(mrmode == 4) then
! do j = 1, N_det_non_ref
! do i_state = 1, N_states
! delta_ij(i_state,j) = delta_ij_mrcc_sto(i_state,j)
! delta_ij_s2(i_state,j) = delta_ij_s2_mrcc_sto(i_state,j)
! enddo
! end do
! else if(mrmode == 10) then
! do j = 1, N_det_non_ref
! do i_state = 1, N_states
! delta_ij(i_state,j) = delta_ij_mrsc2(i_state,j)
! delta_ij_s2(i_state,j) = delta_ij_s2_mrsc2(i_state,j)
! enddo
! end do
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc 5=mrcc_stoch
if(mrmode == 5) then
do j = 1, N_det_non_ref
do i_state = 1, N_states
@ -657,13 +527,6 @@ END_PROVIDER
stop "invalid mrmode"
end if
!if(mrmode == 2 .or. mrmode == 3) then
! do j = 1, N_det_non_ref
! do i_state = 1, N_states
! delta_ij(i_state,j) += delta_ij_cancel(i_state,j)
! enddo
! end do
!end if
END_PROVIDER
@ -1159,7 +1022,7 @@ subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_m
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then
N_tq += 1
do k=1,N_int
do k=1,Nint
tq(k,1,N_tq) = det_buffer(k,1,i)
tq(k,2,N_tq) = det_buffer(k,2,i)
enddo
@ -1247,3 +1110,146 @@ end
subroutine get_cc_coef(tq,c_alpha)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: tq(N_int,2)
double precision, intent(out) :: c_alpha(N_states)
integer :: k
integer :: degree1, degree2, degree
double precision :: hla, hka, dIk(N_states), dka(N_states), dIa(N_states)
double precision :: phase, phase2
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind) :: tmp_det(N_int,2)
integer :: i_state, k_sd, l_sd, i_I
logical :: ok
PROVIDE n_act_orb elec_num
if (perturbative_triples) then
PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
endif
c_alpha(1:N_states) = 0.d0
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq,psi_ref(1,1,i_I),degree1,N_int)
if (degree1 < 3) then
return
endif
enddo
! |I>
do i_I=1,N_det_ref
! Find triples and quadruple grand parents
call get_excitation_degree(tq,psi_ref(1,1,i_I),degree1,N_int)
if (degree1 > 4) then
cycle
endif
if ( (degree1 < 3).or.(degree1 > 4) ) stop 'bug'
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
! <I| <> |alpha>
do k_sd=1,N_det_non_ref
if (maxval(abs(psi_non_ref_coef(k_sd,1:N_states))) < 1.d-10) then
cycle
endif
call get_excitation_degree(tq,psi_non_ref(1,1,k_sd),degree,N_int)
if (degree > 2) then
cycle
endif
call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,k_sd),degree,N_int)
if (degree > 2) then
cycle
endif
do i_state=1,N_states
dIK(i_state) = dij(i_I, k_sd, i_state)
enddo
if (maxval(abs(dIk)) < 1.d-10) then
cycle
endif
! <I| /k\ |alpha>
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_ref(1,1,k_sd),tq,exc,degree2,phase,N_int)
call decode_exc(exc,degree2,h1,p1,h2,p2,s1,s2)
tmp_det(1:N_int,1:2) = psi_ref(1:N_int,1:2,i_I)
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int)
! <I| \l/ |alpha>
dka(1:N_states) = 0.d0
if (ok) then
do l_sd=k_sd+1,N_det_non_ref
if (maxval(abs(psi_non_ref_coef(l_sd,1:N_states))) < 1.d-10) then
cycle
endif
call get_excitation_degree(tmp_det,psi_non_ref(1,1,l_sd),degree,N_int)
if (degree == 0) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,l_sd),exc,degree,phase2,N_int)
do i_state=1,N_states
dka(i_state) = dij(i_I, l_sd, i_state) * phase * phase2
enddo
exit
endif
enddo
else if (perturbative_triples) then
! Linked
call i_h_j(tq,psi_non_ref(1,1,k_sd),N_int,hka)
if (dabs(hka) > 1.d-12) then
double precision :: Delta_E(N_states)
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq,Delta_E)
do i_state=1,N_states
ASSERT (Delta_E(i_state) < 0.d0)
dka(i_state) = hka / Delta_E(i_state)
enddo
endif
endif
if (perturbative_triples.and. (degree2 == 1) ) then
call i_h_j(psi_ref(1,1,i_I),tmp_det,N_int,hka)
call i_h_j(tq,psi_non_ref(1,1,k_sd),N_int,hla)
hka = hla - hka
if (dabs(hka) > 1.d-12) then
call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq,Delta_E)
do i_state=1,N_states
ASSERT (Delta_E(i_state) < 0.d0)
dka(i_state) = hka / Delta_E(i_state)
enddo
endif
endif
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
enddo
enddo
do i_state=1,N_states
c_alpha(i_state) += dIa(i_state) * psi_ref_coef(i_I,i_state)
enddo
enddo
end

8
plugins/mrcepa0/dressing_slave.irp.f
View File

@ -450,15 +450,15 @@ subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ij_,delta_ij_s2_)
do l=1, n(1)
do i_state=1,N_states
delta_ij_(i_state,idx(l,1)) += delta(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state)
delta_ij_s2_(i_state,idx(l,1)) += delta_s2(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state)
delta_ij_(i_state,idx(l,1)) += delta(i_state,l,1) * psi_ref_coef(i_I,i_state)
delta_ij_s2_(i_state,idx(l,1)) += delta_s2(i_state,l,1) * psi_ref_coef(i_I,i_state)
end do
end do
do l=1, n(2)
do i_state=1,N_states
delta_ij_(i_state,idx(l,2)) += delta(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state)
delta_ij_s2_(i_state,idx(l,2)) += delta_s2(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state)
delta_ij_(i_state,idx(l,2)) += delta(i_state,l,2) * psi_ref_coef(J,i_state)
delta_ij_s2_(i_state,idx(l,2)) += delta_s2(i_state,l,2) * psi_ref_coef(J,i_state)
end do
end do

21
plugins/mrcepa0/dressing_vector.irp.f
View File

@ -6,26 +6,25 @@
BEGIN_DOC
! Null dressing vectors
END_DOC
dressing_column_h(:,:) = 0.d0
dressing_column_s(:,:) = 0.d0
integer :: i,ii,k,j,jj, l
double precision :: f, tmp
double precision, external :: u_dot_v
dressing_column_h(:,:) = 0.d0
dressing_column_s(:,:) = 0.d0
do k=1,N_states
l = dressed_column_idx(k)
f = 1.d0/psi_coef(l,k)
do jj = 1, n_det_non_ref
f = -1.d0/psi_coef(l,k)
do jj=1,N_det_non_ref
j = idx_non_ref(jj)
dressing_column_h(j,k) = delta_ij (k,jj) * f
dressing_column_s(j,k) = delta_ij_s2(k,jj) * f
dressing_column_h(j,k) = 2.d0*delta_ij (k,jj)
dressing_column_s(j,k) = 2.d0*delta_ij_s2(k,jj)
dressing_column_h(l,k) += psi_coef(j,k) * delta_ij(k,jj)
dressing_column_s(l,k) += psi_coef(j,k) * delta_ij_s2(k,jj)
enddo
tmp = u_dot_v(dressing_column_h(1,k), psi_coef(1,k), N_det)
dressing_column_h(l,k) -= tmp * f
tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
dressing_column_s(l,k) -= tmp * f
dressing_column_h(l,k) *= f
dressing_column_s(l,k) *= f
enddo
! stop
END_PROVIDER

27
plugins/mrcepa0/mrcc_omp.irp.f
View File

@ -1,27 +0,0 @@
program mrsc2sub
implicit none
double precision, allocatable :: energy(:)
allocate (energy(N_states))
!!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 4
read_wf = .True.
SOFT_TOUCH read_wf
call set_generators_bitmasks_as_holes_and_particles
if (.True.) then
integer :: i,j
do j=1,N_states
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors(i,j)
enddo
enddo
SOFT_TOUCH psi_coef
endif
call run(N_states,energy)
if(do_pt2)then
call run_pt2(N_states,energy)
endif
deallocate(energy)
end

22
plugins/mrcepa0/mrcc_stoch_routines.irp.f
View File

@ -43,6 +43,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error)
call write_double(6,relative_error,"Target relative error")
print *, '========== ================= ================= ================='
print *, ' Samples Energy Stat. Error Seconds '
print *, '========== ================= ================= ================='
@ -177,7 +178,6 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m
usf = 0d0
num = 0d0
print *, "TARGET ERROR :", relative_error
delta = 0d0
delta_s2 = 0d0
allocate(delta_det(N_states, N_det_non_ref, 0:comb_teeth+1, 2))
@ -310,7 +310,6 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m
end if
end do
if(cur_cp == 0) then
print *, "no checkpoint reached so far..."
cycle pullLoop
end if
!!!!!!!!!!!!
@ -337,7 +336,7 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m
print '(I5,F15.7,E12.4,F10.2)', cur_cp, E(mrcc_stoch_istate)+E0+avg, eqt, time-timeInit
end if
if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30) .or. total_computed == N_det_generators) then
if (( (dabs(eqt/(E(mrcc_stoch_istate)+E0+avg)) < relative_error) .and. (cps_N(cur_cp) >= 10) ) .or. total_computed == N_det_generators) then
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
call sleep(1)
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
@ -401,12 +400,13 @@ end function
&BEGIN_PROVIDER [ integer, comb_teeth ]
&BEGIN_PROVIDER [ integer, N_cps_max ]
implicit none
integer :: comb_per_cp
comb_teeth = 16
N_cps_max = 64
!comb_per_cp = 64
! comb_per_cp = 64
gen_per_cp = (N_det_generators / N_cps_max) + 1
N_cps_max += 1
!N_cps_max = N_det_generators / comb_per_cp + 1
! N_cps_max += 1
N_cps_max = N_det_generators / gen_per_cp + 1
END_PROVIDER
@ -526,6 +526,11 @@ subroutine get_comb_val(stato, detail, cur_cp, val)
val = 0d0
first = cp_first_tooth(cur_cp)
!TODO : check
if (first == 0) then
return
endif
do j = comb_teeth, first, -1
!DIR$ FORCEINLINE
k = mrcc_find(curs, mrcc_cweight,size(mrcc_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
@ -562,7 +567,7 @@ subroutine add_comb(com, computed, cp, N, tbc)
implicit none
double precision, intent(in) :: com
integer, intent(inout) :: N
double precision, intent(inout) :: cp(N_det_non_ref)
double precision, intent(inout) :: cp(N_det_generators)
logical, intent(inout) :: computed(N_det_generators)
integer, intent(inout) :: tbc(N_det_generators)
integer :: i, k, l, dets(comb_teeth)
@ -591,7 +596,7 @@ end subroutine
&BEGIN_PROVIDER [ double precision, mrcc_cweight_cache, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, fractage, (comb_teeth) ]
&BEGIN_PROVIDER [ double precision, comb_step ]
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (0:comb_teeth+1) ]
&BEGIN_PROVIDER [ integer, first_det_of_comb ]
&BEGIN_PROVIDER [ integer, tooth_of_det, (N_det_generators) ]
implicit none
@ -650,6 +655,7 @@ end subroutine
end do
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
first_det_of_teeth(1) = first_det_of_comb
first_det_of_teeth(0) = 1
if(first_det_of_teeth(1) /= first_det_of_comb) then

4
plugins/mrcepa0/mrcc_zmq.irp.f
View File

@ -6,6 +6,10 @@ program mrsc2sub
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
mrmode = 5
threshold_generators = 1.d0
threshold_selectors = 1.d0
TOUCH threshold_generators threshold_selectors
read_wf = .True.
SOFT_TOUCH read_wf
call set_generators_bitmasks_as_holes_and_particles

58
plugins/mrcepa0/save_mrcc_wavefunction.irp.f Normal file
View File

@ -0,0 +1,58 @@
program save_mrcc_wf
implicit none
threshold_generators = 1.d0
threshold_selectors = 1.d0
PROVIDE N_int psi_det
TOUCH threshold_generators threshold_selectors
mrmode=5
read_wf = .True.
SOFT_TOUCH read_wf mrmode
call generate_all_alpha_beta_det_products
call run1
call run2
end
subroutine run1
implicit none
integer :: k
double precision :: c_alpha(N_states)
call set_generators_bitmasks_as_holes_and_particles
call get_cc_coef(psi_det(1,1,1), c_alpha)
!$OMP PARALLEL DO DEFAULT(SHARED) &
!$OMP PRIVATE(k,c_alpha) SCHEDULE(static,64)
do k=1,N_det
! if (maxval(abs(psi_coef(k,1:N_states))) == 0.d0) then
if (iand(k,1023) == 0) then
print *, k, '/', N_det
endif
call get_cc_coef(psi_det(1,1,k), c_alpha)
psi_coef(k,1:N_states) = c_alpha(1:N_states)
! endif
enddo
!$OMP END PARALLEL DO
SOFT_TOUCH psi_coef
end
subroutine run2
implicit none
integer :: k
double precision :: c_alpha(N_states)
psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
psi_coef(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
! do k=N_det,1,-1
! if (maxval(abs(psi_coef(k,1:N_states))) > 0.d0) then
! exit
! endif
! enddo
! N_det = k
SOFT_TOUCH N_det psi_coef psi_det
call save_wavefunction
end

1
plugins/pyscf/PyscfToQp.py
View File

@ -1,3 +1,4 @@
#!/usr/bin/env python2
import numpy,re,sys

17
plugins/read_integral/Gen_Ezfio_from_integral.sh
View File

@ -1,17 +0,0 @@
#!/bin/bash
ezfio=$1
# Create the integral
echo 'Create Integral'
echo 'Create EZFIO'
read nel nmo natom <<< $(cat param)
read e_nucl <<< $(cat e_nuc)
./create_ezfio.py $ezfio $nel $natom $nmo $e_nucl
#Handle the orbital consitensy check
qp_edit -c $ezfio &> /dev/null
cp $ezfio/{ao,mo}_basis/ao_md5
#Read the integral
echo 'Read Integral'
qp_run read_integrals_mo $ezfio

23
plugins/shiftedbk/EZFIO.cfg
View File

@ -1,5 +1,28 @@
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[energy_pt2]
type: double precision
doc: Calculated energy with pt2 contribution
interface: ezfio
[thresh_dressed_ci]
type: Threshold
doc: Threshold on the convergence of the dressed CI energy
interface: ezfio,provider,ocaml
default: 1.e-5
[n_it_max_dressed_ci]
type: Strictly_positive_int
doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 30
[h0_type]
type: Perturbation
doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
interface: ezfio,provider,ocaml
default: EN

86
plugins/shiftedbk/selection_buffer.irp.f
View File

@ -156,3 +156,89 @@ end subroutine
subroutine unique_selection_buffer(b)
use selection_types
implicit none
BEGIN_DOC
! Removes duplicate determinants in the selection buffer
END_DOC
type(selection_buffer), intent(inout) :: b
integer, allocatable :: iorder(:)
integer(bit_kind), pointer :: detmp(:,:,:)
double precision, pointer :: val(:)
integer :: i,j,k
integer(bit_kind), allocatable :: bit_tmp(:)
logical,allocatable :: duplicate(:)
logical :: found_duplicates
integer*8, external :: det_search_key
if (b%N == 0 .or. b%cur == 0) return
allocate (duplicate(b%cur), val(size(b%val)), detmp(N_int, 2, size(b%val)), bit_tmp(b%cur))
call sort_dets_by_det_search_key(b%cur, b%det, b%val, detmp, val, 1)
deallocate(b%det, b%val)
do i=b%cur+1,b%N
val(i) = 0.d0
detmp(1:N_int,1:2,i) = 0_bit_kind
enddo
b%det => detmp
b%val => val
do i=1,b%cur
bit_tmp(i) = det_search_key(b%det(1,1,i),N_int)
duplicate(i) = .False.
enddo
do i=1,b%cur-1
if (duplicate(i)) then
cycle
endif
j = i+1
do while (bit_tmp(j)==bit_tmp(i))
if (duplicate(j)) then
j += 1
if (j > b%cur) then
exit
else
cycle
endif
endif
duplicate(j) = .True.
do k=1,N_int
if ( (b%det(k,1,i) /= b%det(k,1,j) ) &
.or. (b%det(k,2,i) /= b%det(k,2,j) ) ) then
duplicate(j) = .False.
exit
endif
enddo
j += 1
if (j > b%cur) then
exit
endif
enddo
enddo
found_duplicates = .False.
do i=1,b%cur
if (duplicate(i)) then
found_duplicates = .True.
exit
endif
enddo
if (found_duplicates) then
k=0
do i=1,N_det
if (.not.duplicate(i)) then
k += 1
b%det(:,:,k) = b%det(:,:,i)
b%val(k) = b%val(i)
endif
enddo
b%cur = k
endif
deallocate (duplicate,bit_tmp)
end

159
plugins/shiftedbk/shifted_bk_iter.irp.f Normal file
View File

@ -0,0 +1,159 @@
program shifted_bk
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:)
integer :: degree
integer :: n_det_before
double precision :: threshold_davidson_in
allocate (pt2(N_states))
double precision :: hf_energy_ref
logical :: has
double precision :: relative_error, absolute_error
integer :: N_states_p
character*(512) :: fmt
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order
pt2 = -huge(1.e0)
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
call diagonalize_CI_dressed
call save_wavefunction
call ezfio_has_hartree_fock_energy(has)
if (has) then
call ezfio_get_hartree_fock_energy(hf_energy_ref)
else
hf_energy_ref = ref_bitmask_energy
endif
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
call diagonalize_CI_dressed
call save_wavefunction
N_states_p = min(N_det,N_states)
endif
n_det_before = 0
character*(8) :: pt2_string
double precision :: threshold_selectors_save, threshold_generators_save
threshold_selectors_save = threshold_selectors
threshold_generators_save = threshold_generators
double precision :: error(N_states), energy(N_states)
error = 0.d0
threshold_selectors = 1.d0
threshold_generators = 1d0
if (.True.) then
pt2_string = '(sh-Bk) '
do while ( (N_det < N_det_max) )
write(*,'(A)') '--------------------------------------------------------------------------------'
N_det_delta_ij = N_det
do i=1,N_states
energy(i) = psi_energy(i)+nuclear_repulsion
enddo
PROVIDE delta_ij_tmp
call delta_ij_done()
call diagonalize_ci_dressed
do i=1,N_states
pt2(i) = ci_energy_dressed(i) - energy(i)
enddo
N_states_p = min(N_det,N_states)
print *, ''
print '(A,I12)', 'Summary at N_det = ', N_det
print '(A)', '-----------------------------------'
print *, ''
print *, ''
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
write(*,fmt) ('State',k, k=1,N_states_p)
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
write(*,fmt) '# E ', energy(1:N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', energy(1:N_states_p)-energy(1)
write(*,fmt) '# Excit. (eV)', (energy(1:N_states_p)-energy(1))*27.211396641308d0
endif
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (energy(k)+pt2(k),error(k), k=1,N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', ( (energy(k)+pt2(k)-energy(1)-pt2(1)), &
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
write(*,fmt) '# Excit. (eV)', ( (energy(k)+pt2(k)-energy(1)-pt2(1))*27.211396641308d0, &
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
endif
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
print *, ''
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1, N_states_p
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', energy(k)
print *, 'E+PT2'//pt2_string//' = ', energy(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print *, 'Variational Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (energy(i) - energy(1)), &
(energy(i) - energy(1)) * 27.211396641308d0
enddo
print *, '-----'
print*, 'Variational + perturbative Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (energy(i)+ pt2(i) - (energy(1) + pt2(1))), &
(energy(i)+ pt2(i) - (energy(1) + pt2(1))) * 27.211396641308d0
enddo
endif
call ezfio_set_shiftedbk_energy_pt2(energy(1)+pt2(1))
! call dump_fci_iterations_value(N_det,energy,pt2)
n_det_before = N_det
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det >= N_det_max) then
threshold_davidson = threshold_davidson_in
end if
call save_wavefunction
call ezfio_set_shiftedbk_energy(energy(1))
call ezfio_set_shiftedbk_energy_pt2(ci_energy_dressed(1))
enddo
endif
end

51
plugins/shiftedbk/shifted_bk_routines.irp.f
View File

@ -16,21 +16,22 @@ END_PROVIDER
&BEGIN_PROVIDER [ type(selection_buffer), global_sb ]
&BEGIN_PROVIDER [ type(selection_buffer), mini_sb ]
&BEGIN_PROVIDER [ double precision, N_det_increase_factor ]
implicit none
integer :: i
implicit none
fock_diag_tmp_(:,:,:) = 0.d0
integer :: i
N_det_increase_factor = 1d0
N_det_increase_factor = dble(N_states)
n_det_add = max(1, int(float(N_det) * N_det_increase_factor))
call create_selection_buffer(n_det_add, n_det_add*2, global_sb)
call create_selection_buffer(n_det_add, n_det_add*2, mini_sb)
do i=1,Nproc
call create_selection_buffer(n_det_add, n_det_add*2, sb(i))
end do
a_h_i = 0d0
a_s2_i = 0d0
END_PROVIDER
n_det_add = max(1, int(float(N_det) * N_det_increase_factor))
call create_selection_buffer(n_det_add, n_det_add*2, global_sb)
call create_selection_buffer(n_det_add, n_det_add*2, mini_sb)
do i=1,Nproc
call create_selection_buffer(n_det_add, n_det_add*2, sb(i))
end do
a_h_i = 0d0
a_s2_i = 0d0
END_PROVIDER
BEGIN_PROVIDER [ integer, N_dress_int_buffer ]
@ -50,6 +51,7 @@ subroutine generator_done(i_gen, int_buf, double_buf, det_buf, N_buf, iproc)
double precision, intent(out) :: double_buf(N_dress_double_buffer)
integer(bit_kind), intent(out) :: det_buf(N_int, 2, N_dress_det_buffer)
integer :: i
int_buf(:) = 0
call sort_selection_buffer(sb(iproc))
@ -118,15 +120,17 @@ subroutine delta_ij_done()
old_det_gen = N_det_generators
! Add buffer only when the last state is computed
call unique_selection_buffer(global_sb)
call sort_selection_buffer(global_sb)
call fill_H_apply_buffer_no_selection(global_sb%cur,global_sb%det,N_int,0)
call copy_H_apply_buffer_to_wf()
if (s2_eig.or.(N_states > 1) ) then
call make_s2_eigenfunction
endif
call undress_with_alpha(old_generators, old_det_gen, psi_det(1,1,N_det_delta_ij+1), N_det-N_det_delta_ij)
call save_wavefunction
end subroutine
@ -253,7 +257,7 @@ subroutine dress_with_alpha_(Nstates,Ndet,Nint,delta_ij_loc,minilist, det_minili
contrib = 0d0
do i=1,Nstates
de = E0_denominator(i) - haa
de = dress_E0_denominator(i) - haa
if(DABS(de) < 1D-5) cycle
c_alpha(i) = a_h_psi(i) / de
@ -312,23 +316,4 @@ BEGIN_PROVIDER [ logical, initialize_E0_denominator ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
implicit none
BEGIN_DOC
! E0 in the denominator of the PT2
END_DOC
if (initialize_E0_denominator) then
if (h0_type == "EN") then
E0_denominator(1:N_states) = psi_energy(1:N_states)
else if (h0_type == "Barycentric") then
E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states)
else
print *, h0_type, ' not implemented'
stop
endif
else
E0_denominator = -huge(1.d0)
endif
END_PROVIDER

209
scripts/ezfio_interface/qp_convert_output_to_ezfio.py
View File

@ -28,7 +28,8 @@ except:
sys.exit(1)
else:
sys.path = [ QP_ROOT + "/install/EZFIO/Python",
QP_ROOT + "/resultsFile",
QP_ROOT + "/install/resultsFile",
QP_ROOT + "/install",
QP_ROOT + "/scripts"] + sys.path
# ~#~#~#~#~#~ #
@ -92,6 +93,7 @@ def write_ezfio(res, filename):
coord_y.append(a.coord[1] / a0)
coord_z.append(a.coord[2] / a0)
# ~#~#~#~#~ #
# W r i t e #
# ~#~#~#~#~ #
@ -125,9 +127,7 @@ def write_ezfio(res, filename):
coefficient = []
exponent = []
res.clean_contractions()
res.convert_to_cartesian()
# ~#~#~#~#~#~#~ #
# P a r s i n g #
# ~#~#~#~#~#~#~ #
@ -262,159 +262,74 @@ def write_ezfio(res, filename):
ezfio.set_mo_basis_mo_occ(OccNum)
ezfio.set_mo_basis_mo_coef(MoMatrix)
# ______ _
# | ___ \ | |
# | |_/ /__ ___ _ _ __| | ___
# | __/ __|/ _ \ | | |/ _` |/ _ \
# | | \__ \ __/ |_| | (_| | (_) |
# \_| |___/\___|\__,_|\__,_|\___/
#
# INPUT
# {% for lanel,zcore, l_block in l_atom $}
# #local l_block l=0}
# {label} GEN {zcore} {len(l_block)-1 #lmax_block}
# {% for l_param in l_block%}
# {len(l_param) # list of parameter aka n_max_bock_max(n)}
# {% for coef,n,zeta for l_param}
# {coef,n, zeta}
# OUTPUT
# Local are 1 array padded by max(n_max_block) when l == 0 (output:k_loc_max)
# v_k[n-2][atom] = value
#Non Local are 2 array padded with max of lmax_block when l!=0 (output:lmax+1) and max(n_max_block)whem l !=0 (kmax)
# v_kl[l][n-2][atom] = value
def pad(array, size, value=0):
new_array = array
for add in xrange(len(array), size):
new_array.append(value)
return new_array
def parse_str(pseudo_str):
'''Return 4d array atom,l,n, attribute (attribute is coef, n, zeta)'''
matrix = []
array_l_max_block = []
array_z_remove = []
for block in [b for b in pseudo_str.split('\n\n') if b]:
#First element is header, the rest are l_param
array_party = [i for i in re.split(r"\n\d+\n", block) if i]
z_remove, l_max_block = map(int, array_party[0].split()[-2:])
array_l_max_block.append(l_max_block)
array_z_remove.append(z_remove)
x = []
for l in array_party[1:]:
y = []
for coef_n_zeta in l.split('\n'):
z = coef_n_zeta.split()
if z : y.append(z)
x.append(y)
matrix.append(x)
return (matrix, array_l_max_block, array_z_remove)
def get_local_stuff(matrix):
matrix_local_unpad = [atom[0] for atom in matrix]
k_loc_max = max(len(i) for i in matrix_local_unpad)
matrix_local = [ pad(ll, k_loc_max, [0., 2, 0.]) for ll in matrix_local_unpad]
m_coef = [[float(i[0]) for i in atom] for atom in matrix_local]
m_n = [[int(i[1]) - 2 for i in atom] for atom in matrix_local]
m_zeta = [[float(i[2]) for i in atom] for atom in matrix_local]
return (k_loc_max, m_coef, m_n, m_zeta)
def get_non_local_stuff(matrix):
matrix_unlocal_unpad = [atom[1:] for atom in matrix]
l_max_block = max(len(i) for i in matrix_unlocal_unpad)
k_max = max([len(item) for row in matrix_unlocal_unpad for item in row])
matrix_unlocal_semipaded = [[pad(item, k_max, [0., 2, 0.]) for item in row] for row in matrix_unlocal_unpad]
empty_row = [[0., 2, 0.] for k in range(l_max_block)]
matrix_unlocal = [ pad(ll, l_max_block, empty_row) for ll in matrix_unlocal_semipaded ]
m_coef_noloc = [[[float(k[0]) for k in j] for j in i] for i in matrix_unlocal]
m_n_noloc = [[[int(k[1]) - 2 for k in j] for j in i] for i in matrix_unlocal]
m_zeta_noloc = [[[float(k[2]) for k in j] for j in i] for i in matrix_unlocal]
return (l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc)
try:
pseudo_str = []
label = ezfio.get_nuclei_nucl_label()
lmax = 0
nucl_charge_remove = []
klocmax = 0
kmax = 0
nucl_num = len(res.geometry)
for ecp in res.pseudo:
pseudo_str += [ "%(label)s GEN %(zcore)d %(lmax)d" % { "label": label[ ecp["atom"]-1 ],
"zcore": ecp["zcore"], "lmax": ecp["lmax"] } ]
lmax = ecp["lmax"]
for l in [lmax] + list(range(0,lmax)):
pseudo_str += [ "%d"%len(ecp[str(l)]) ]
for t in ecp[str(l)]:
pseudo_str += [ "%f %d %f"%t ]
pseudo_str += [""]
pseudo_str = "\n".join(pseudo_str)
lmax_local = ecp['lmax']
lmax = max(lmax_local,lmax)
nucl_charge_remove.append(ecp['zcore'])
klocmax = max(klocmax, len(ecp[str(lmax_local)]))
for l in range(lmax_local):
kmax = max(kmax,len(ecp[str(l)]))
lmax = lmax-1
ezfio.set_pseudo_pseudo_lmax(lmax)
ezfio.set_pseudo_nucl_charge_remove(nucl_charge_remove)
ezfio.set_pseudo_pseudo_klocmax(klocmax)
ezfio.set_pseudo_pseudo_kmax(kmax)
pseudo_n_k = [ [ 0 for _ in range(nucl_num) ] for _ in range(klocmax) ]
pseudo_v_k = [ [ 0. for _ in range(nucl_num) ] for _ in range(klocmax) ]
pseudo_dz_k = [ [ 0. for _ in range(nucl_num) ] for _ in range(klocmax) ]
pseudo_n_kl = [ [ [ 0 for _ in range(nucl_num) ] for _ in range(kmax) ] for _ in range(lmax+1) ]
pseudo_v_kl = [ [ [ 0. for _ in range(nucl_num) ] for _ in range(kmax) ] for _ in range(lmax+1) ]
pseudo_dz_kl = [ [ [ 0. for _ in range(nucl_num) ] for _ in range(kmax) ] for _ in range(lmax+1) ]
for ecp in res.pseudo:
lmax_local = ecp['lmax']
klocmax = len(ecp[str(lmax_local)])
atom = ecp['atom']-1
for kloc in range(klocmax):
try:
v, n, dz = ecp[str(lmax_local)][kloc]
pseudo_n_k[kloc][atom] = n-2
pseudo_v_k[kloc][atom] = v
pseudo_dz_k[kloc][atom] = dz
except:
pass
for l in range(lmax_local):
for k in range(kmax):
try:
v, n, dz = ecp[str(l)][k]
pseudo_n_kl[l][k][atom] = n-2
pseudo_v_kl[l][k][atom] = v
pseudo_dz_kl[l][k][atom] = dz
except:
pass
ezfio.set_pseudo_pseudo_n_k(pseudo_n_k)
ezfio.set_pseudo_pseudo_v_k(pseudo_v_k)
ezfio.set_pseudo_pseudo_dz_k(pseudo_dz_k)
ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
n_alpha = res.num_alpha
n_beta = res.num_beta
for i in range(nucl_num):
charge[i] -= nucl_charge_remove[i]
n_alpha -= nucl_charge_remove[i]/2
n_beta -= nucl_charge_remove[i]/2
ezfio.set_nuclei_nucl_charge(charge)
ezfio.set_electrons_elec_alpha_num(n_alpha)
ezfio.set_electrons_elec_beta_num(n_beta)
matrix, array_l_max_block, array_z_remove = parse_str(pseudo_str)
array_z_remove = map(float,array_z_remove)
except:
ezfio.set_pseudo_do_pseudo(False)
else:
ezfio.set_pseudo_do_pseudo(True)
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
# Z _ e f f , a l p h a / b e t a _ e l e c #
# ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ #
ezfio.set_pseudo_nucl_charge_remove(array_z_remove)
charge = ezfio.get_nuclei_nucl_charge()
charge = [ i - j for i, j in zip(charge, array_z_remove) ]
ezfio.set_nuclei_nucl_charge (charge)
import math
num_elec_diff = sum(array_z_remove)/2
nalpha = ezfio.get_electrons_elec_alpha_num() - num_elec_diff
nbeta = ezfio.get_electrons_elec_beta_num() - num_elec_diff
ezfio.set_electrons_elec_alpha_num(nalpha)
ezfio.set_electrons_elec_beta_num( nbeta )
# Change all the array 'cause EZFIO
# v_kl (v, l) => v_kl(l,v)
# v_kl => zip(*_v_kl)
# [[7.0, 79.74474797, -49.45159098], [1.0, 5.41040609, -4.60151975]]
# [(7.0, 1.0), (79.74474797, 5.41040609), (-49.45159098, -4.60151975)]
# ~#~#~#~#~ #
# L o c a l #
# ~#~#~#~#~ #
klocmax, m_coef, m_n, m_zeta = get_local_stuff(matrix)
ezfio.pseudo_pseudo_klocmax = klocmax
ezfio.pseudo_pseudo_v_k = zip(*m_coef)
ezfio.pseudo_pseudo_n_k = zip(*m_n)
ezfio.pseudo_pseudo_dz_k = zip(*m_zeta)
# ~#~#~#~#~#~#~#~#~ #
# N o n _ L o c a l #
# ~#~#~#~#~#~#~#~#~ #
l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc = get_non_local_stuff(
matrix)
ezfio.pseudo_pseudo_lmax = l_max_block - 1
ezfio.pseudo_pseudo_kmax = k_max
ezfio.pseudo_pseudo_v_kl = zip(*m_coef_noloc)
ezfio.pseudo_pseudo_n_kl = zip(*m_n_noloc)
ezfio.pseudo_pseudo_dz_kl = zip(*m_zeta_noloc)
def get_full_path(file_path):

23
scripts/qp_set_frozen_core.py
View File

@ -6,21 +6,26 @@ import sys
sys.path = [ os.environ["QP_ROOT"]+"/install/EZFIO/Python" ] + sys.path
from ezfio import ezfio
ezfio.set_filename(sys.argv[1])
filename = sys.argv[1]
if filename == '-q': filename = sys.argv[2]
ezfio.set_filename(filename)
nb = 0
for charge in ezfio.nuclei_nucl_charge:
if charge < 5:
pass
elif charge < 13:
nb += 1
else:
nb += 5
if not ezfio.pseudo_do_pseudo:
for charge in ezfio.nuclei_nucl_charge:
if charge < 5:
pass
elif charge < 13:
nb += 1
else:
nb += 5
mo_tot_num = ezfio.mo_basis_mo_tot_num
if len(sys.argv)>2:
if sys.argv[2] == '-q':
if '-q' in sys.argv:
print nb
sys.exit(0)

21
src/Bitmask/bitmasks.irp.f
View File

@ -21,14 +21,17 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int) ]
! Bitmask to include all possible MOs
END_DOC
integer :: i,j,n
n = mod(mo_tot_num-1,bit_kind_size)+1
full_ijkl_bitmask = 0_bit_kind
do i=1,N_int-1
full_ijkl_bitmask(i) = not(0_bit_kind)
enddo
do i=1,n
full_ijkl_bitmask(N_int) = ibset(full_ijkl_bitmask(N_int),i-1)
integer :: i,j,k
k=0
do j=1,N_int
full_ijkl_bitmask(j) = 0_bit_kind
do i=0,bit_kind_size-1
k=k+1
if (mo_class(k) /= 'Deleted') then
full_ijkl_bitmask(j) = ibset(full_ijkl_bitmask(j),i)
endif
if (k == mo_tot_num) exit
enddo
enddo
END_PROVIDER
@ -559,7 +562,7 @@ END_PROVIDER
n_core_inact_act_orb = 0
do i = 1, N_int
reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1))
reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1))
reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,2,1))
n_core_inact_act_orb +=popcnt(reunion_of_core_inact_act_bitmask(i,1))
enddo
END_PROVIDER

2
src/Davidson/davidson_parallel.irp.f
View File

@ -403,8 +403,8 @@ BEGIN_PROVIDER [ integer, nthreads_davidson ]
call getenv('NTHREADS_DAVIDSON',env)
if (trim(env) /= '') then
read(env,*) nthreads_davidson
call write_int(6,nthreads_davidson,'Number of threads for <Psi|H|Psi>')
endif
call write_int(6,nthreads_davidson,'Number of threads for Diagonalization')
END_PROVIDER

54
src/Davidson/diagonalization_hs2_dressed.irp.f
View File

@ -35,7 +35,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag)
integer, intent(in) :: dressing_state
double precision, allocatable :: H_jj(:), S2_jj(:)
double precision, allocatable :: H_jj(:)
double precision, external :: diag_H_mat_elem, diag_S_mat_elem
integer :: i,k
@ -44,7 +44,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
PROVIDE mo_bielec_integrals_in_map
allocate(H_jj(sze),S2_jj(sze))
allocate(H_jj(sze))
H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
!$OMP PARALLEL DEFAULT(NONE) &
@ -60,17 +60,32 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d
if (dressing_state > 0) then
do k=1,N_st
do i=1,sze
H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
H_jj(i) += u_in(i,k) * dressing_column_h(i,k)
enddo
enddo
endif
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
deallocate (H_jj,S2_jj)
integer :: N_st_diag_local
double precision, allocatable :: energies_local(:), s2_out_local(:), u_in_local(:,:)
logical :: converged
converged = .False.
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state,converged)
N_st_diag_local = N_st_diag
do while (.not.converged)
N_st_diag_local = 2 * N_st_diag_local
allocate (energies_local(N_st_diag_local), s2_out_local(N_st_diag_local), u_in_local(sze,N_st_diag_local))
u_in_local(1:sze,1:N_st_diag) = u_in(1:sze,1:N_st_diag)
call davidson_diag_hjj_sjj(dets_in,u_in_local,H_jj,s2_out_local,energies_local,dim_in,sze,N_st,N_st_diag_local,Nint,dressing_state,converged)
energies(1:N_st_diag) = energies_local(1:N_st_diag)
s2_out(1:N_st_diag) = s2_out_local(1:N_st_diag)
u_in(1:sze,1:N_st_diag) = u_in_local(1:sze,1:N_st_diag)
deallocate (energies_local, s2_out_local, u_in_local)
enddo
deallocate (H_jj)
end
subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state)
subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state,converged)
use bitmasks
implicit none
BEGIN_DOC
@ -105,13 +120,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
integer :: iter
integer :: i,j,k,l,m
logical :: converged
logical, intent(inout) :: converged
double precision, external :: u_dot_v, u_dot_u
integer :: k_pairs, kl
integer :: iter2
integer :: iter2, itertot
double precision, allocatable :: W(:,:), U(:,:), S(:,:), overlap(:,:)
double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:)
double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:)
@ -123,6 +138,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
integer :: shift, shift2, itermax, istate
double precision :: r1, r2
logical :: state_ok(N_st_diag*davidson_sze_max)
integer :: nproc_target
include 'constants.include.F'
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, h, lambda
@ -133,6 +149,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
endif
itermax = max(3,min(davidson_sze_max, sze/N_st_diag))
itertot = 0
PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse
@ -146,8 +163,22 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
call write_int(6,N_st,'Number of states')
call write_int(6,N_st_diag,'Number of states in diagonalization')
call write_int(6,sze,'Number of determinants')
nproc_target = nproc
r1 = 8.d0*(3.d0*dble(sze*N_st_diag*itermax+5.d0*(N_st_diag*itermax)**2 &
+ 4.d0*(N_st_diag*itermax)+nproc*(4.d0*N_det_alpha_unique+2.d0*N_st_diag*sze)))/(1024.d0**3)
+ 3.d0*(N_st_diag*itermax)+nproc*(4.d0*N_det_alpha_unique+2.d0*N_st_diag*sze)))/(1024.d0**3)
if (qp_max_mem > 0) then
do while (r1 > qp_max_mem)
nproc_target = nproc_target - 1
r1 = 8.d0*(3.d0*dble(sze*N_st_diag*itermax+5.d0*(N_st_diag*itermax)**2 &
+ 3.d0*(N_st_diag*itermax)+nproc_target*(4.d0*N_det_alpha_unique+2.d0*N_st_diag*sze)))/(1024.d0**3)
if (nproc_target == 0) then
nproc_target = 1
exit
endif
enddo
call omp_set_num_threads(nproc_target)
call write_int(6,nproc_target,'Number of threads for diagonalization')
endif
call write_double(6, r1, 'Memory(Gb)')
write(6,'(A)') ''
write_buffer = '====='
@ -220,6 +251,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
do while (.not.converged)
itertot = itertot+1
if (itertot == 8) then
exit
endif
do k=1,N_st_diag
do i=1,sze
@ -492,6 +527,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
y, s_, s_tmp, &
lambda &
)
call omp_set_num_threads(nproc)
end

1
src/Davidson/u0Hu0.irp.f
View File

@ -28,6 +28,7 @@ subroutine H_S2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze)
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det),v_t(N_st,N_det),s_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo

7
src/Determinants/EZFIO.cfg
View File

@ -40,6 +40,13 @@ doc: Force the wave function to be an eigenfunction of S^2
interface: ezfio,provider,ocaml
default: True
[used_weight]
type: integer
doc: 0: 1/(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
interface: ezfio,provider,ocaml
default: 0
[threshold_generators]
type: Threshold
doc: Thresholds on generators (fraction of the norm)

11
src/Determinants/Fock_diag.irp.f
View File

@ -21,10 +21,19 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint)
if (Ne(1) /= elec_alpha_num) then
print *, 'Error in build_fock_tmp (alpha)', Ne(1), Ne(2)
call debug_det(det_ref,N_int)
! print *, occ(:,1)
! print *, occ(:,2)
! do i=1,10000
! occ(i,1) = fock_diag_tmp(1,mo_tot_num+i) !traceback
! enddo
stop -1
endif
if (Ne(2) /= elec_beta_num) then
print *, 'Error in build_fock_tmp (beta)', Ne(1), Ne(2)
! print *, 'Error in build_fock_tmp (beta)', Ne(1), Ne(2)
! do i=1,10000
! occ(i,1) = fock_diag_tmp(1,mo_tot_num+i) !traceback
! enddo
stop -1
endif

68
src/Determinants/density_matrix.irp.f
View File

@ -347,7 +347,8 @@ subroutine set_natural_mos
double precision, allocatable :: tmp(:,:)
label = "Natural"
call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label)
call mo_as_svd_vectors_of_mo_matrix_eig(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,mo_occ,label)
soft_touch mo_occ
end
subroutine save_natural_mos
@ -361,6 +362,57 @@ subroutine save_natural_mos
end
BEGIN_PROVIDER [ double precision, l3_weight, (N_states) ]
implicit none
BEGIN_DOC
! Weight of the states in the selection : 1/(sum_i |c_i|^3)
END_DOC
integer :: i,k
double precision :: c
do i=1,N_states
l3_weight(i) = 1.d-31
do k=1,N_det
c = psi_coef(k,i)*psi_coef(k,i)
l3_weight(i) = l3_weight(i) + c*abs(psi_coef(k,i))
enddo
l3_weight(i) = min(1.d0/l3_weight(i), 100.d0)
enddo
if (mpi_master) then
print *, ''
print *, 'L3 weights'
print *, '----------'
print *, ''
print *, l3_weight(1:N_states)
print *, ''
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ]
implicit none
BEGIN_DOC
! Weight of the states in the selection : 1/c_0^2
END_DOC
integer :: i,k
double precision :: c
do i=1,N_states
c0_weight(i) = 1.d-31
c = maxval(psi_coef(:,i) * psi_coef(:,i))
c0_weight(i) = 1.d0/c
c0_weight(i) = min(c0_weight(i), 100.d0)
enddo
if (mpi_master) then
print *, ''
print *, 'c0 weights'
print *, '----------'
print *, ''
print *, c0_weight(1:N_states)
print *, ''
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
implicit none
BEGIN_DOC
@ -369,9 +421,17 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
logical :: exists
state_average_weight(:) = 1.d0
call ezfio_has_determinants_state_average_weight(exists)
if (exists) then
call ezfio_get_determinants_state_average_weight(state_average_weight)
if (used_weight == 0) then
state_average_weight(:) = c0_weight(:)
else if (used_weight == 1) then
state_average_weight(:) = 1./N_states
else if (used_weight == 3) then
state_average_weight(:) = l3_weight
else
call ezfio_has_determinants_state_average_weight(exists)
if (exists) then
call ezfio_get_determinants_state_average_weight(state_average_weight)
endif
endif
state_average_weight(:) = state_average_weight(:)+1.d-31
state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:)))

36
src/Determinants/determinants.irp.f
View File

@ -321,21 +321,24 @@ end subroutine
END_DOC
call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, &
psi_det_sorted_bit, psi_coef_sorted_bit)
psi_det_sorted_bit, psi_coef_sorted_bit, N_states)
END_PROVIDER
subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out)
subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out, N_st)
use bitmasks
implicit none
integer, intent(in) :: Ndet
integer, intent(in) :: Ndet, N_st
integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size)
double precision , intent(in) :: coef_in(psi_det_size,N_states)
double precision , intent(in) :: coef_in(psi_det_size,N_st)
integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size)
double precision , intent(out) :: coef_out(psi_det_size,N_states)
double precision , intent(out) :: coef_out(psi_det_size,N_st)
BEGIN_DOC
! Determinants are sorted are sorted according to their det_search_key.
! Useful to accelerate the search of a random determinant in the wave
! function.
!
! /!\ The first dimension of coef_out and coef_in need to be psi_det_size
!
END_DOC
integer :: i,j,k
integer, allocatable :: iorder(:)
@ -356,7 +359,7 @@ subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out)
det_out(j,1,i) = det_in(j,1,iorder(i))
det_out(j,2,i) = det_in(j,2,iorder(i))
enddo
do k=1,N_states
do k=1,N_st
coef_out(i,k) = coef_in(iorder(i),k)
enddo
enddo
@ -451,6 +454,26 @@ end
subroutine save_wavefunction_truncated(thr)
implicit none
double precision, intent(in) :: thr
use bitmasks
BEGIN_DOC
! Save the wave function into the EZFIO file
END_DOC
integer :: N_det_save,i
N_det_save = N_det
do i=1,N_det
if (psi_average_norm_contrib_sorted(i) < thr) then
N_det_save = i
exit
endif
enddo
if (mpi_master) then
call save_wavefunction_general(N_det_save,min(N_states,N_det_save),psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
endif
end
subroutine save_wavefunction
implicit none
use bitmasks
@ -718,6 +741,7 @@ subroutine apply_excitation(det, exc, res, ok, Nint)
case default
print *, degree
print *, "apply ex"
! print *, 1.d0/0.d0 ! For traceback
STOP
end select
! END INLINE

4
src/Determinants/psi_cas.irp.f
View File

@ -54,7 +54,7 @@ END_PROVIDER
! function.
END_DOC
call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, &
psi_cas_sorted_bit, psi_cas_coef_sorted_bit)
psi_cas_sorted_bit, psi_cas_coef_sorted_bit, N_states)
END_PROVIDER
@ -107,7 +107,7 @@ END_PROVIDER
! function.
END_DOC
call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, &
psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit)
psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit, N_states)
END_PROVIDER

9
src/Determinants/s2.irp.f
View File

@ -189,8 +189,8 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
allocate(vt(sze_8,N_st))
vt = 0.d0
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,1)
!$OMP DO SCHEDULE(static,1)
do sh2=sh,shortcut(0,1)
exa = 0
do ni=1,Nint
@ -227,11 +227,11 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
enddo
enddo
enddo
!$OMP END DO NOWAIT
enddo
!$OMP END DO NOWAIT
!$OMP DO SCHEDULE(dynamic)
do sh=1,shortcut(0,2)
!$OMP DO
do i=shortcut(sh,2),shortcut(sh+1,2)-1
org_i = sort_idx(i,2)
do j=shortcut(sh,2),i-1
@ -249,8 +249,9 @@ subroutine S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
end if
end do
end do
!$OMP END DO NOWAIT
enddo
!$OMP END DO NOWAIT
!$OMP BARRIER
do istate=1,N_st
do i=n,1,-1

16
src/Determinants/slater_rules.irp.f
View File

@ -639,11 +639,12 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
select case (degree)
case (2)
call get_double_excitation(key_i,key_j,exc,phase,Nint)
! Mono alpha, mono beta
if (exc(0,1,1) == 1) then
! Mono alpha, mono beta
if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then
s2 = -phase
endif
if(exc(1,1,1) == exc(1,2,2) )then
if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!!
hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
else if (exc(1,2,1) ==exc(1,1,2))then
hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
@ -654,8 +655,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
exc(1,2,1), &
exc(1,2,2) ,mo_integrals_map)
endif
! Double alpha
else if (exc(0,1,1) == 2) then
! Double alpha
hij = phase*(get_mo_bielec_integral( &
exc(1,1,1), &
exc(2,1,1), &
@ -666,8 +667,8 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
exc(2,1,1), &
exc(2,2,1), &
exc(1,2,1) ,mo_integrals_map) )
! Double beta
else if (exc(0,1,2) == 2) then
! Double beta
hij = phase*(get_mo_bielec_integral( &
exc(1,1,2), &
exc(2,1,2), &
@ -683,13 +684,13 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
! Mono alpha
if (exc(0,1,1) == 1) then
! Mono alpha
m = exc(1,1,1)
p = exc(1,2,1)
spin = 1
! Mono beta
else
! Mono beta
m = exc(1,1,2)
p = exc(1,2,2)
spin = 2
@ -697,7 +698,6 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
case (0)
print *," ZERO"
double precision, external :: diag_S_mat_elem
s2 = diag_S_mat_elem(key_i,Nint)
hij = diag_H_mat_elem(key_i,Nint)

13
src/Determinants/spindeterminants.irp.f
View File

@ -367,13 +367,13 @@ end
j = psi_bilinear_matrix_columns(k)
f = 0.d0
do l=1,N_states
f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l)
f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l) * state_average_weight(l)
enddo
det_alpha_norm(i) += f
det_beta_norm(j) += f
enddo
det_alpha_norm = det_alpha_norm / dble(N_states)
det_beta_norm = det_beta_norm / dble(N_states)
det_alpha_norm = det_alpha_norm
det_beta_norm = det_beta_norm
END_PROVIDER
@ -649,7 +649,7 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
do k=1,N_states
psi_coef_sorted_bit(idx,k) = psi_bilinear_matrix(i,j,k)
!$OMP ATOMIC
norm(k) += psi_bilinear_matrix(i,j,k)
norm(k) += psi_bilinear_matrix(i,j,k)*psi_bilinear_matrix(i,j,k)
enddo
endif
enddo
@ -671,7 +671,10 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate)
do i=1,N_det
norm(1) += psi_average_norm_contrib_sorted(i)
if (truncate) then
if (norm(1) >= 0.999999d0) then
if (norm(1) >= 1.d0) then
exit
endif
if (psi_average_norm_contrib_sorted(i) == 0.d0) then
exit
endif
endif

4
src/Determinants/zmq.irp.f
View File

@ -304,7 +304,7 @@ integer function zmq_get_psi_det(zmq_to_qp_run_socket, worker_id)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to broadcast zmq_get_psi_det'
endif
call broadcast_chunks_bit_kind(psi_det,N_det*N_int*2)
call broadcast_chunks_bit_kind(psi_det,size(psi_det))
IRP_ENDIF
end
@ -353,7 +353,7 @@ integer function zmq_get_psi_coef(zmq_to_qp_run_socket, worker_id)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to broadcast zmq_get_psi_coef'
endif
call broadcast_chunks_double(psi_coef,N_states*N_det)
call broadcast_chunks_double(psi_coef,size(psi_coef))
IRP_ENDIF
end

59
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -100,7 +100,7 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1)
call bielec_integrals_index(i(ii),j(ii),k(ii),l(ii),i2)
if (i1 /= i2) then
print *, i1, i2
print *, i(ii), j(jj), k(jj), l(jj)
print *, i(ii), j(ii), k(ii), l(ii)
stop 'bielec_integrals_index_reverse failed'
endif
endif
@ -421,8 +421,8 @@ double precision function mo_bielec_integral(i,j,k,l)
integer, intent(in) :: i,j,k,l
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
!DIR$ FORCEINLINE
PROVIDE mo_bielec_integrals_in_map
!DIR$ FORCEINLINE
mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
return
end
@ -438,24 +438,45 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE mo_bielec_integrals_in_map
double precision, external :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
integer :: ii, ii0
integer*8 :: ii_8, ii0_8
real(integral_kind) :: tmp
integer(key_kind) :: i1, idx
integer(key_kind) :: p,q,r,s,i2
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
ii0 = l-mo_integrals_cache_min
ii0 = ior(ii0, k-mo_integrals_cache_min)
ii0 = ior(ii0, j-mo_integrals_cache_min)
ii0_8 = int(l,8)-mo_integrals_cache_min_8
ii0_8 = ior( ishft(ii0_8,7), int(k,8)-mo_integrals_cache_min_8)
ii0_8 = ior( ishft(ii0_8,7), int(j,8)-mo_integrals_cache_min_8)
q = min(j,l)
s = max(j,l)
q = q+ishft(s*s-s,-1)
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash(i))
ii = ior(ii0, i-mo_integrals_cache_min)
if (iand(ii, -128) == 0) then
ii_8 = ior( ishft(ii0_8,7), int(i,8)-mo_integrals_cache_min_8)
out_val(i) = mo_integrals_cache(ii_8)
else
p = min(i,k)
r = max(i,k)
p = p+ishft(r*r-r,-1)
i1 = min(p,q)
i2 = max(p,q)
idx = i1+ishft(i2*i2-i2,-1)
!DIR$ FORCEINLINE
call map_get(map,idx,tmp)
out_val(i) = dble(tmp)
endif
enddo
if (key_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
subroutine get_mo_bielec_integrals_ij(k,l,sze,out_array,map)
@ -534,7 +555,7 @@ subroutine get_mo_bielec_integrals_coulomb_ii(k,l,sze,out_val,map)
call bielec_integrals_index(k,i,l,i,hash(i))
enddo
if (key_kind == 8) then
if (integral_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)
@ -567,7 +588,7 @@ subroutine get_mo_bielec_integrals_exch_ii(k,l,sze,out_val,map)
call bielec_integrals_index(k,i,i,l,hash(i))
enddo
if (key_kind == 8) then
if (integral_kind == 8) then
call map_get_many(map, hash, out_val, sze)
else
call map_get_many(map, hash, tmp_val, sze)

2
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -240,7 +240,7 @@ subroutine add_integrals_to_map(mask_ijkl)
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
character*(2048) :: output(1)
print*, 'i'
print *, 'i'
call bitstring_to_str( output(1), mask_ijkl(1,1), N_int )
print *, trim(output(1))
j = 0

26
src/Integrals_Bielec/qp_ao_ints.irp.f
View File

@ -1,10 +1,16 @@
program qp_ao_ints
use omp_lib
implicit none
IRP_IF MPI
include 'mpif.h'
IRP_ENDIF
integer :: ierr
BEGIN_DOC
! Increments a running calculation to compute AO integrals
END_DOC
integer :: i
PROVIDE zmq_context mpi_master zmq_state zmq_context
call switch_qp_run_to_master
@ -17,15 +23,27 @@ program qp_ao_ints
double precision :: integral, ao_bielec_integral
integral = ao_bielec_integral(1,1,1,1)
character*(64) :: state
call wait_for_state(zmq_state,state)
do while (state /= 'Stopped')
do
call wait_for_state('ao_integrals',zmq_state)
if (zmq_state(1:7) == 'Stopped') then
exit
endif
!$OMP PARALLEL DEFAULT(PRIVATE) PRIVATE(i)
i = omp_get_thread_num()
call ao_bielec_integrals_in_map_slave_tcp(i)
!$OMP END PARALLEL
call wait_for_state(zmq_state,state)
IRP_IF MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here, 'error in barrier'
endif
IRP_ENDIF
enddo
IRP_IF MPI
call MPI_finalize(i)
IRP_ENDIF
print *, 'Done'
end

4
src/MO_Basis/ao_ortho_canonical.irp.f
View File

@ -2,9 +2,7 @@
&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num ]
implicit none
BEGIN_DOC
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_cart_to_sphe_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
! ao_cart_to_sphe coefficients of the current basis set
END_DOC
integer :: i
integer, external :: ao_power_index

70
src/MO_Basis/utils.irp.f
View File

@ -112,6 +112,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
double precision, intent(in) :: matrix(lda,n)
integer :: i,j
double precision :: accu
double precision, allocatable :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
@ -137,12 +138,16 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
write (6,'(A)') 'Eigenvalues'
write (6,'(A)') '-----------'
write (6,'(A)') ''
write (6,'(A)') '======== ================'
write (6,'(A)') '======== ================ ================'
write (6,'(A)') ' MO Eigenvalue Cumulative '
write (6,'(A)') '======== ================ ================'
accu = 0.d0
do i=1,m
write (6,'(I8,1X,F16.10)') i,D(i)
accu = accu + D(i)
write (6,'(I8,1X,F16.10,1X,F16.10)') i,D(i), accu
enddo
write (6,'(A)') '======== ================'
write (6,'(A)') '======== ================ ================'
write (6,'(A)') ''
call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
@ -152,6 +157,65 @@ subroutine mo_as_svd_vectors_of_mo_matrix(matrix,lda,m,n,label)
mo_label = label
end
subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
implicit none
integer,intent(in) :: lda,m,n
character*(64), intent(in) :: label
double precision, intent(in) :: matrix(lda,n)
double precision, intent(out) :: eig(m)
integer :: i,j
double precision :: accu
double precision, allocatable :: mo_coef_new(:,:), U(:,:),D(:), A(:,:), Vt(:,:), work(:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, U, Vt, A
call write_time(6)
if (m /= mo_tot_num) then
print *, irp_here, ': Error : m/= mo_tot_num'
stop 1
endif
allocate(A(lda,n),U(lda,n),mo_coef_new(ao_num,m),D(m),Vt(lda,n))
do j=1,n
do i=1,m
A(i,j) = matrix(i,j)
enddo
enddo
mo_coef_new = mo_coef
call svd(A,lda,U,lda,D,Vt,lda,m,n)
write (6,'(A)') 'MOs are now **'//trim(label)//'**'
write (6,'(A)') ''
write (6,'(A)') 'Eigenvalues'
write (6,'(A)') '-----------'
write (6,'(A)') ''
write (6,'(A)') '======== ================ ================'
write (6,'(A)') ' MO Eigenvalue Cumulative '
write (6,'(A)') '======== ================ ================'
accu = 0.d0
do i=1,m
accu = accu + D(i)
write (6,'(I8,1X,F16.10,1X,F16.10)') i,D(i), accu
enddo
write (6,'(A)') '======== ================ ================'
write (6,'(A)') ''
call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
do i=1,m
eig(i) = D(i)
enddo
deallocate(A,mo_coef_new,U,Vt,D)
call write_time(6)
mo_label = label
end
subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
implicit none
integer,intent(in) :: n,m

2
src/Utils/fortran_mmap.c
View File

@ -33,7 +33,7 @@ void* mmap_fortran(char* filename, size_t bytes, int* file_descr, int read_only)
exit(EXIT_FAILURE);
}
result = lseek(fd, bytes, SEEK_SET);
result = lseek(fd, bytes+1, SEEK_SET);
if (result == -1) {
close(fd);
printf("%s:\n", filename);

Some files were not shown because too many files have changed in this diff Show More