diff --git a/Makefile b/Makefile index cdbe67c5..bce7735c 100644 --- a/Makefile +++ b/Makefile @@ -32,6 +32,7 @@ bin/irpf90: (echo Unable to download IRPF90 : $(WWW_SERVER)/$(IRPF90_TGZ) ; exit 1) tar -zxf $(IRPF90_TGZ) && rm $(IRPF90_TGZ) $(MAKE) -C irpf90 | tee install_irpf90.log + rm -rf -- $$PWD/bin/irpf90 $$PWD/bin/irpman ln -s $$PWD/irpf90/bin/irpf90 $$PWD/bin/irpf90 ln -s $$PWD/irpf90/bin/irpman $$PWD/bin/irpman @@ -53,7 +54,9 @@ bin/m4: QPACKAGE_ROOT=$$PWD ./scripts/install_m4.sh | tee install_m4.log -ocaml: ocaml/Qptypes.ml curl m4 +ocaml: curl m4 + - rm -f -- ocaml/Qptypes.ml + $(MAKE) ocaml/Qptypes.ml ocaml/Qptypes.ml: $(info $(BLUE)===== Installing ocaml =====$(BLACK)) diff --git a/README.md b/README.md new file mode 100644 index 00000000..3354ea1f --- /dev/null +++ b/README.md @@ -0,0 +1,8 @@ +Quantum package +=============== + +[![Gitter](https://badges.gitter.im/Join Chat.svg)](https://gitter.im/LCPQ/quantum_package?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) + +Set of quantum chemistry programs and libraries. + +For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>) diff --git a/README.rst b/README.rst deleted file mode 100644 index 5dfa5b65..00000000 --- a/README.rst +++ /dev/null @@ -1,9 +0,0 @@ -=============== -Quantum package -=============== - -Set of quantum chemistry programs and libraries. - -For more information, you can visit the -`wiki of the project `_ - diff --git a/data/basis/vdz b/data/basis/vdz new file mode 100644 index 00000000..8a1ac65b --- /dev/null +++ b/data/basis/vdz @@ -0,0 +1,1042 @@ +HYDROGEN +S 3 + 1 13.0107010 0.19682158E-01 + 2 1.9622572 0.13796524 + 3 0.44453796 0.47831935 +S 1 + 1 0.12194962 1.0000000 + +HELIUM +S 3 + 1 38.354936737 0.23814288905E-01 + 2 5.7689081479 0.15490906777 + 3 1.2399407035 0.46998096633 +S 1 + 1 0.29757815953 1.0000000 + +LITHIUM +S 5 + 1 266.27785516 0.64920150325E-02 + 2 40.069783447 0.47747863215E-01 + 3 9.0559944389 0.20268796111 + 4 2.4503009051 0.48606574817 + 5 0.72209571855 0.43626977955 +S 1 + 1 0.052810884721 1.0000000 +S 1 + 1 0.020960948798 1.0000000 + +BERYLLIUM +S 5 + 1 515.18616125 0.55615307983E-02 + 2 77.511037595 0.41190068062E-01 + 3 17.552481693 0.17913378108 + 4 4.8028940596 0.44736716455 + 5 1.4516214316 0.42009581920 +S 1 + 1 0.13281633620 1.0000000 +S 1 + 1 0.045837372213 1.0000000 + +BORON +S 5 + 1 839.31830086 -0.55929201074E-02 + 2 126.26464843 -0.41565520771E-01 + 3 28.620600763 -0.18299816983 + 4 7.8793722710 -0.46540391866 + 5 2.4088857172 -0.44173884791 +S 1 + 1 0.25105109036 1.0000000 +S 1 + 1 0.083648866069 1.0000000 +P 3 + 1 6.0332223619 -0.35603672456E-01 + 2 1.2499157866 -0.19895775769 + 3 0.33871676350 -0.50850202618 +P 1 + 1 0.096415632351 1.0000000 + +CARBON +S 5 + 1 1238.4016938 0.0054568832082 + 2 186.29004992 0.0406384092110 + 3 42.251176346 0.18025593888 + 4 11.676557932 0.46315121755 + 5 3.5930506482 0.44087173314 +S 1 + 1 0.40245147363 1.0000000 +S 1 + 1 0.13090182668 1.0000000 +P 3 + 1 9.4680970621 0.038387871728 + 2 2.0103545142 0.21117025112 + 3 0.54771004707 0.51328172114 +P 1 + 1 0.15268613795 1.0000000 + +NITROGEN +S 5 + 1 1712.8415853 -0.53934125305E-02 + 2 257.64812677 -0.40221581118E-01 + 3 58.458245853 -0.17931144990 + 4 16.198367905 -0.46376317823 + 5 5.0052600809 -0.44171422662 +S 1 + 1 0.58731856571 1.0000000 +S 1 + 1 0.18764592253 1.0000000 +P 3 + 1 13.571470233 -0.40072398852E-01 + 2 2.9257372874 -0.21807045028 + 3 0.79927750754 -0.51294466049 +P 1 + 1 0.21954348034 1.0000000 + +OXYGEN +S 5 + 1 2266.1767785 -0.53431809926E-02 + 2 340.87010191 -0.39890039230E-01 + 3 77.363135167 -0.17853911985 + 4 21.479644940 -0.46427684959 + 5 6.6589433124 -0.44309745172 +S 1 + 1 0.80975975668 1.0000000 +S 1 + 1 0.25530772234 1.0000000 +P 3 + 1 17.721504317 0.43394573193E-01 + 2 3.8635505440 0.23094120765 + 3 1.0480920883 0.51375311064 +P 1 + 1 0.27641544411 1.0000000 + +FLUORINE +S 5 + 1 2894.8325990 -0.53408255515E-02 + 2 435.41939120 -0.39904258866E-01 + 3 98.843328866 -0.17912768038 + 4 27.485198001 -0.46758090825 + 5 8.5405498171 -0.44653131020 +S 1 + 1 1.0654578038 1.0000000 +S 1 + 1 0.33247346748 1.0000000 +P 3 + 1 22.696633924 -0.45212874436E-01 + 2 4.9872339257 -0.23754317067 + 3 1.3491613954 -0.51287353587 +P 1 + 1 0.34829881977 1.0000000 + +NEON +S 5 + 1 3598.9736625 -0.53259297003E-02 + 2 541.32073112 -0.39817417969E-01 + 3 122.90450062 -0.17914358188 + 4 34.216617022 -0.46893582977 + 5 10.650584124 -0.44782537577 +S 1 + 1 1.3545953960 1.0000000 +S 1 + 1 0.41919362639 1.0000000 +P 3 + 1 28.424053785 -0.46031944795E-01 + 2 6.2822510953 -0.23993183041 + 3 1.6978715079 -0.50871724964 +P 1 + 1 0.43300700172 1.0000000 + +SODIUM +S 5 + 1 4098.2003908 -0.58535911879E-02 + 2 616.49374031 -0.43647161872E-01 + 3 139.96644001 -0.19431465884 + 4 39.073441051 -0.48685065731 + 5 11.929847205 -0.41881705137 +S 3 + 1 20.659966030 0.85949689854E-01 + 2 1.9838860978 -0.56359144041 + 3 0.64836323942 -0.51954009048 +S 1 + 1 0.52443967404E-01 1.0000000 +S 1 + 1 0.28048160742E-01 1.0000000 +P 5 + 1 75.401862017 0.154353625324E-01 + 2 17.274818978 0.997382931840E-01 + 3 5.1842347425 0.312095939659 + 4 1.6601211973 0.492956748074 + 5 0.51232528958 0.324203983180 + +MAGNESIUM +S 5 + 1 4953.8339196 -0.57778967498E-02 + 2 745.18044154 -0.43124761082E-01 + 3 169.21604972 -0.19268216987 + 4 47.300672019 -0.48641439116 + 5 14.461336973 -0.42550894077 +S 3 + 1 24.768174789 0.87956969984E-01 + 2 2.4940945349 -0.55165058128 + 3 0.87807584533 -0.53443294833 +S 1 + 1 0.87212782497E-01 1.0000000 +S 1 + 1 0.33599293780E-01 1.0000000 +P 5 + 1 98.053010494 -0.14480564601E-01 + 2 22.586932277 -0.95495750787E-01 + 3 6.8391509842 -0.30787672651 + 4 2.2332843818 -0.49936292886 + 5 0.71606599387 -0.31503476213 + +ALUMINUM +S 5 + 1 5887.5727030 0.13483347987E-02 + 2 885.61225996 0.10071576809E-01 + 3 201.13604899 0.45132454056E-01 + 4 56.284974674 0.11461268043 + 5 17.229551243 0.10159608943 +S 3 + 1 29.340249922 0.69347454208E-01 + 2 3.0439630420 -0.42528117679 + 3 1.1285539518 -0.41449832210 +S 1 + 1 0.14234175160 1.0000000 +S 1 + 1 0.54400192313E-01 1.0000000 +P 5 + 1 145.11918809 0.63963373134E-02 + 2 33.717894833 0.44189359965E-01 + 3 10.369863083 0.15581575993 + 4 3.5135616036 0.28635286951 + 5 1.1980050273 0.22921423248 +P 1 + 1 0.26583005913 1.0000000 +P 1 + 1 0.71003361994E-01 1.0000000 + +SILICON +S 5 + 1 6903.7118686 0.13373962995E-02 + 2 1038.4346419 0.99966546241E-02 + 3 235.87581480 0.44910165101E-01 + 4 66.069385169 0.11463638540 + 5 20.247945761 0.10280063858 +S 3 + 1 34.353481730 0.70837285010E-01 + 2 3.6370788192 -0.43028836252 + 3 1.4002048599 -0.41382774969 +S 1 + 1 0.20484414805 1.0000000 +S 1 + 1 0.77994095468E-01 1.0000000 +P 5 + 1 179.83907373 0.61916656462E-02 + 2 41.907258846 0.43399431982E-01 + 3 12.955294367 0.15632019351 + 4 4.4383267393 0.29419996982 + 5 1.5462247904 0.23536823814 +P 1 + 1 0.35607612302 1.0000000 +P 1 + 1 0.10008513762 1.0000000 + +PHOSPHOROUS +S 5 + 1 8002.4795106 0.57503489302E-02 + 2 1203.6813590 0.43007628959E-01 + 3 273.44227031 0.19363986252 + 4 76.655541517 0.49651693399 + 5 23.516927435 0.44983262479 +S 3 + 1 39.791683439 0.95188129789E-01 + 2 4.2770343323 -0.57649840368 + 3 1.6940256888 -0.54239583865 +S 1 + 1 0.27567674644 1.0000000 +S 1 + 1 0.10495590099 1.0000000 +P 5 + 1 219.50755823 0.92100565257E-02 + 2 51.274155030 0.65409765745E-01 + 3 15.921595892 0.24033730279 + 4 5.5069913481 0.46318321788 + 5 1.9537719426 0.37392563382 +P 1 + 1 0.47803397932 1.0000000 +P 1 + 1 0.13657952621 1.0000000 + +SULFUR +S 5 + 1 9184.9303010 -0.22294387756E-02 + 2 1381.5105503 -0.16683029937E-01 + 3 313.87147580 -0.75262436116E-01 + 4 88.053870623 -0.19376827038 + 5 27.039914905 -0.17718020803 +S 3 + 1 45.648731303 -0.10736062573 + 2 4.9664522326 0.65066293018 + 3 2.0116242047 0.59712155354 +S 1 + 1 0.35661077013 1.0000000 +S 1 + 1 0.13507221477 1.0000000 +P 5 + 1 261.98233439 -0.92729929822E-02 + 2 61.306894736 -0.66547669241E-01 + 3 19.103729887 -0.24828595903 + 4 6.6567720378 -0.48703847402 + 5 2.3959635161 -0.39337850312 +P 1 + 1 0.61776161679 1.0000000 +P 1 + 1 0.16993376871 1.0000000 + +CHLORINE +S 5 + 1 10449.8275660 0.19708362484E-02 + 2 1571.7365221 0.14754727977E-01 + 3 357.12065523 0.66679112875E-01 + 4 100.25185935 0.17228924084 + 5 30.812727554 0.15883786100 +S 3 + 1 51.923789434 -0.10009298909 + 2 5.7045760975 0.60841752753 + 3 2.3508376809 0.54352153355 +S 1 + 1 0.44605124672 1.0000000 +S 1 + 1 0.16848856190 1.0000000 +P 5 + 1 307.66790569 -0.87801484118E-02 + 2 72.102015515 -0.63563355471E-01 + 3 22.532680262 -0.24016428276 + 4 7.8991765444 -0.47798866557 + 5 2.8767268321 -0.38515850005 +P 1 + 1 0.77459363955 1.0000000 +P 1 + 1 0.21037699698 1.0000000 + +ARGON +S 5 + 1 11797.1197650 0.0020214479973 + 2 1774.3522753 0.0151398526880 + 3 403.18875733 0.0685254008730 + 4 113.24933999 0.17762929207 + 5 34.835298193 0.16496495001 +S 3 + 1 58.614775042 -0.10343394000 + 2 6.4922045078 0.63133365899 + 3 2.7117014403 0.54887572301 +S 1 + 1 0.54412974538 1.0000000 +S 1 + 1 0.20517411136 1.0000000 +P 5 + 1 356.28707256 -0.87321783271E-02 + 2 83.593132866 -0.63680317444E-01 + 3 26.186704029 -0.24311906886 + 4 9.2257275925 -0.48956069653 + 5 3.3922754954 -0.39229190040 +P 1 + 1 0.94740534010 1.0000000 +P 1 + 1 0.25659135062 1.0000000 + +POTASSIUM +S 6 + 1 31478.7467640 0.39838653994E-02 + 2 4726.8876066 0.30501759762E-01 + 3 1075.4345353 0.15073752622 + 4 303.39811023 0.51912939801 + 5 98.327112831 1.0366957005 + 6 33.636222177 0.76398963199 +S 3 + 1 65.639209962 -0.28242617106 + 2 7.3162592218 1.6914935860 + 3 2.8902580135 1.2965331953 +S 3 + 1 4.5459748965 -0.76343555273E-02 + 2 0.70404124062 0.25635718960E-01 + 3 0.28266888959 0.16606859208E-01 +S 1 + 1 0.29058164020E-01 1.0000000 +S 1 + 1 0.12111638157E-01 1.0000000 +P 5 + 1 361.22492154 0.20906479823E-01 + 2 84.670222166 0.15043641740 + 3 26.469088236 0.55440061077 + 4 9.2658077615 1.0409009991 + 5 3.3423388293 0.67825341194 +P 3 + 1 1.5100876104 0.75248191146 + 2 0.56568375163 1.3708585031 + 3 0.20817008495 0.66047633079 + +CALCIUM +S 6 + 1 35138.7139290 0.39482520740E-02 + 2 5276.4111348 0.30234243552E-01 + 3 1200.4692589 0.14952019681 + 4 338.71810542 0.51597345713 + 5 109.85385922 1.0339510296 + 6 37.608880299 0.76937933526 +S 3 + 1 73.107977555 -0.28268525011 + 2 8.2407705688 1.6796092142 + 3 3.2959812993 1.2803766016 +S 3 + 1 5.2341800914 -0.76868604561E-02 + 2 0.84187220515 0.25382375978E-01 + 3 0.36510294029 0.16512171511E-01 +S 1 + 1 0.51222402884E-01 1.0000000 +S 1 + 1 0.19825111408E-01 1.0000000 +P 5 + 1 413.11313893 0.20327135354E-01 + 2 96.935786224 0.14730276362 + 3 30.372154659 0.54887167322 + 4 10.684776830 1.0440659818 + 5 3.8821258350 0.68653490684 +P 3 + 1 1.7993016295 0.75410246871 + 2 0.69189056530 1.3409296599 + 3 0.26364024096 0.56391989435 + +SCANDIUM +S 6 + 1 38956.0818040 0.39293205858E-02 + 2 5849.5733637 0.30093221229E-01 + 3 1330.8813154 0.14890370950 + 4 375.55534165 0.51464282906 + 5 121.87261370 1.0337708070 + 6 41.760243729 0.77436851627 +S 3 + 1 81.060633953 -0.28318552316 + 2 9.2059823972 1.6770806984 + 3 3.7063215732 1.2594733678 +S 3 + 1 5.9888909988 -0.77821367493E-02 + 2 0.97363432378 0.25499692745E-01 + 3 0.42041019223 0.16191560820E-01 +S 1 + 1 0.059440551913 1.0000000 +S 1 + 1 0.022897806303 1.0000000 +P 5 + 1 466.31481262 0.019984300385 + 2 109.51217097 0.14561043072 + 3 34.375921827 0.54687466223 + 4 12.142096955 1.0479006012 + 5 4.4336767669 0.68894890327 +P 3 + 1 2.0971291866 0.75619214724 + 2 0.80977606956 1.3178212235 + 3 0.30834046588 0.54312268173 +D 4 + 1 19.240334928 0.27039082144E-01 + 2 5.1178995899 0.13803684743 + 3 1.6554278827 0.34869086403 + 4 0.54016355610 0.48594185717 +D 1 + 1 0.16211214518 1.0000000 + +TITANIUM +S 6 + 1 42961.5121850 0.39127635355E-02 + 2 6450.9759169 0.29969820489E-01 + 3 1467.7210915 0.14836352707 + 4 414.20997355 0.51347285324 + 5 134.48715840 1.0335365483 + 6 46.122209796 0.77854233930 +S 3 + 1 89.447762543 -0.28385401259 + 2 10.223346060 1.6772785333 + 3 4.1353774271 1.2411928456 +S 3 + 1 6.7896181452 -0.78399994518E-02 + 2 1.1106730691 0.25495493019E-01 + 3 0.47565975578 0.16061172892E-01 +S 1 + 1 0.65986956934E-01 1.0000000 +S 1 + 1 0.25210342250E-01 1.0000000 +P 5 + 1 522.03684782 0.19754179642E-01 + 2 122.68649489 0.14460677619 + 3 38.572903611 0.54669004165 + 4 13.672169319 1.0531647540 + 5 5.0118529359 0.69111213363 +P 3 + 1 2.4131928282 0.75803437136 + 2 0.93252270050 1.3036241399 + 3 0.35429058390 0.53638653300 +D 4 + 1 23.465125957 0.26536380115E-01 + 2 6.3332593832 0.13796453963 + 3 2.0766489946 0.35312644228 + 4 0.69027361954 0.48647124166 +D 1 + 1 0.21088738554 1.0000000 + +VANADIUM +S 6 + 1 47160.3760600 0.14498688908E-02 + 2 7081.4110871 0.11106435251E-01 + 3 1611.1621223 0.55005423585E-01 + 4 454.72940551 0.19060252591 + 5 147.71321208 0.38435022957 + 6 50.699538950 0.29095546792 +S 3 + 1 98.262492669 -0.10942337856 + 2 11.294293099 0.64539490399 + 3 4.5853360105 0.47117880777 +S 3 + 1 7.6359689588 -0.22454949054 + 2 1.2539836689 0.72594852764 + 3 0.53271935387 0.45560582748 +S 1 + 1 0.72246239566E-01 1.0000000 +S 1 + 1 0.27358087448E-01 1.0000000 +P 5 + 1 580.55044988 0.97315110917E-02 + 2 136.52341127 0.71531241137E-01 + 3 42.983958820 0.27197688414 + 4 15.282798763 0.52618988893 + 5 5.6202495154 0.34452533498 +P 3 + 1 2.7485386415 0.34040396496 + 2 1.0618550073 0.57983996120 + 3 0.40235518645 0.23911643083 +D 4 + 1 27.358434017 0.26641927050E-01 + 2 7.4540604253 0.13995311726 + 3 2.4633917847 0.35751066639 + 4 0.82480925277 0.48488354148 +D 1 + 1 0.25257904742 1.0000000 + +CHROMIUM +S 6 + 1 51528.0863490 0.14405823106E-02 + 2 7737.2103487 0.11036202287E-01 + 3 1760.3748470 0.54676651806E-01 + 4 496.87706544 0.18965038103 + 5 161.46520598 0.38295412850 + 6 55.466352268 0.29090050668 +S 3 + 1 107.54732999 -0.10932281100 + 2 12.408671897 0.64472599471 + 3 5.0423628826 0.46262712560 +S 3 + 1 8.5461640165 -0.22711013286 + 2 1.3900441221 0.73301527591 + 3 0.56066602876 0.44225565433 +S 1 + 1 0.71483705972E-01 1.0000000 +S 1 + 1 0.28250687604E-01 1.0000000 +P 5 + 1 640.48536096 0.96126715203E-02 + 2 150.69711194 0.70889834655E-01 + 3 47.503755296 0.27065258990 + 4 16.934120165 0.52437343414 + 5 6.2409680590 0.34107994714 +P 3 + 1 3.0885463206 0.33973986903 + 2 1.1791047769 0.57272062927 + 3 0.43369774432 0.24582728206 +D 4 + 1 27.559479426 0.30612488044D-01 + 2 7.4687020327 0.15593270944 + 3 2.4345903574 0.36984421276 + 4 0.78244754808 0.47071118077 +D 1 + 1 0.21995774311 1.0000000 + +MANGANESE +S 6 + 1 56137.0090370 0.14321304702E-02 + 2 8429.2063943 0.10972509162E-01 + 3 1917.8277233 0.54382468712E-01 + 4 541.36230198 0.18884335129 + 5 176.00069142 0.38198025054 + 6 60.500477010 0.29156772596 +S 3 + 1 117.17282882 -0.10933661328 + 2 13.596973368 0.64305039431 + 3 5.5483996341 0.45848970584 +S 3 + 1 9.4662853530 -0.22538977259 + 2 1.5595006070 0.72307758657 + 3 0.65230205868 0.45300721536 +S 1 + 1 0.84003734475E-01 1.0000000 +S 1 + 1 0.31256098581E-01 1.0000000 +P 5 + 1 706.00497535 0.95055518167E-02 + 2 166.19728820 0.70356271142E-01 + 3 52.452061906 0.27005556982 + 4 18.746932862 0.52574344602 + 5 6.9282991622 0.34254033223 +P 3 + 1 3.4772204938 0.33994073736 + 2 1.3406906449 0.57203836254 + 3 0.50498803038 0.23847605831 +D 4 + 1 35.423264935 0.26985304111E-01 + 2 9.7814221451 0.14383458648 + 3 3.2673488767 0.36418958377 + 4 1.1026472189 0.48152670661 +D 1 + 1 0.33743205934 1.0000000 + +IRON +S 6 + 1 60923.6406430 0.14302254466E-02 + 2 9147.8893982 0.10958790038E-01 + 3 2081.3505927 0.54332554248E-01 + 4 587.55977067 0.18884995009 + 5 191.09043990 0.38253069946 + 6 65.732730112 0.29308335984 +S 3 + 1 127.25891928 -0.10964564925 + 2 14.830913010 0.64387631332 + 3 6.0653307408 0.45472347323 +S 3 + 1 10.449943710 -0.22539639952 + 2 1.7245228003 0.72164398156 + 3 0.71772177325 0.44985492922 +S 1 + 1 0.91449828308E-01 1.0000000 +S 1 + 1 0.33706691021E-01 1.0000000 +P 5 + 1 773.43750995 0.94325735144E-02 + 2 182.15149714 0.70029620575E-01 + 3 57.547272758 0.26993651996 + 4 20.614988935 0.52700011047 + 5 7.6348557890 0.34284148028 +P 3 + 1 3.8719327990 0.33974402988 + 2 1.4924724132 0.56842594005 + 3 0.56061284958 0.23649365839 +D 4 + 1 38.968133419 0.27879664382E-01 + 2 10.800067078 0.14858319982 + 3 3.6136457999 0.36905479496 + 4 1.2129967888 0.47745100883 +D 1 + 1 0.36524393170 1.0000000 + +COBALT +S 6 + 1 65902.2082570 0.14284614936E-02 + 2 9895.3896027 0.10946072783E-01 + 3 2251.4305789 0.54285953890E-01 + 4 635.61097084 0.18885179079 + 5 206.78820681 0.38301634994 + 6 71.179242971 0.29443551266 +S 3 + 1 137.77268040 -0.10990221736 + 2 16.118079243 0.64455537395 + 3 6.6030327710 0.45116787924 +S 3 + 1 11.479915788 -0.22593846910 + 2 1.8956426324 0.72231409008 + 3 0.78466232067 0.44903812296 +S 1 + 1 0.98425774432E-01 1.0000000 +S 1 + 1 0.35945741932E-01 1.0000000 +P 5 + 1 843.64358575 0.93866097254E-02 + 2 198.76386994 0.69880208716E-01 + 3 62.854963098 0.27037070345 + 4 22.562842280 0.52904786880 + 5 8.3713209127 0.34357029579 +P 3 + 1 4.2858719800 0.34027999036 + 2 1.6508041817 0.56693392384 + 3 0.61834231096 0.23617979783 +D 4 + 1 42.927867612 0.28487788365E-01 + 2 11.942533053 0.15206951283 + 3 4.0046495664 0.37310913999 + 4 1.3413193804 0.47549837676 +D 1 + 1 0.40015009743 1.0000000 + +NICKEL +S 6 + 1 71074.8032110 0.14260386729E-02 + 2 10672.0209410 0.10928236994E-01 + 3 2428.1389007 0.54212626938E-01 + 4 685.53595148 0.18874768902 + 5 223.10072863 0.38324616985 + 6 76.842014042 0.29550637144 +S 3 + 1 148.71122016 -0.11014443059 + 2 17.459154987 0.64521426988 + 3 7.1625280665 0.44797838103 +S 3 + 1 12.556137125 -0.22645403224 + 2 2.0735740488 0.72320959286 + 3 0.85382640602 0.44868026476 +S 1 + 1 0.10536766271 1.0000000 +S 1 + 1 0.38134087688E-01 1.0000000 +P 5 + 1 916.73608662 0.93439635610E-02 + 2 216.06139913 0.69737374902E-01 + 3 68.383914817 0.27073495012 + 4 24.593843952 0.53078301549 + 5 9.1392960204 0.34410229438 +P 3 + 1 4.7193371746 0.34076082016 + 2 1.8161849234 0.56580169611 + 3 0.67840750720 0.23616717361 +D 4 + 1 47.093832108 0.28982316948E-01 + 2 13.146463975 0.15494995950 + 3 4.4170548925 0.37633115111 + 4 1.4771565078 0.47365096014 +D 1 + 1 0.43735921792 1.0000000 + +COPPER +S 6 + 1 76381.3480560 0.14336079896E-02 + 2 11468.7774990 0.10986749865E-01 + 3 2609.4246495 0.54513652465E-01 + 4 736.75033098 0.18990128258 + 5 239.82419958 0.38581959211 + 6 82.656829252 0.29790607498 +S 3 + 1 160.13544196 -0.11146778567 + 2 18.834177695 0.65349301031 + 3 7.7176595741 0.44770534421 +S 3 + 1 13.710846717 -0.22870911122 + 2 2.2349895670 0.73464423031 + 3 0.87818360069 0.43273070874 +S 1 + 1 0.87187458064E-01 1.0000000 +S 1 + 1 0.32969114665E-01 1.0000000 +P 5 + 1 991.24075782 0.93878498798E-02 + 2 233.69376116 0.70208282458E-01 + 3 74.020930927 0.27323522220 + 4 26.664967447 0.53580792728 + 5 9.9192087478 0.34575794906 +P 3 + 1 5.1519553926 0.34229108083 + 2 1.9638205828 0.56456592484 + 3 0.71560097037 0.24078584318 +D 4 + 1 47.335049590 0.32375547758E-01 + 2 13.161666077 0.16810218684 + 3 4.3693777244 0.38477707982 + 4 1.4132925109 0.46147880178 +D 1 + 1 0.38878001452 1.0000000 + +ZINC +S 6 + 1 82000.7116290 0.14210764000E-02 + 2 12312.4717770 0.10891499487E-01 + 3 2801.3944193 0.54057188059E-01 + 4 790.99424302 0.18847463904 + 5 257.56551079 0.38346549346 + 6 88.814933400 0.29723794197 +S 3 + 1 171.86353716 -0.11051849523 + 2 20.302534785 0.64607716984 + 3 8.3464123068 0.44220117322 +S 3 + 1 14.847536940 -0.22705309278 + 2 2.4495029507 0.72433217935 + 3 0.99845821824 0.44836495592 +S 1 + 1 0.11891307937 1.0000000 +S 1 + 1 0.42297428760E-01 1.0000000 +P 5 + 1 1071.5185372 0.92767797235E-02 + 2 252.69712152 0.69541149434E-01 + 3 80.100829126 0.27156772564 + 4 28.903393172 0.53401355573 + 5 10.768899879 0.34501323446 +P 3 + 1 5.6446212530 0.34129600164 + 2 2.1678291347 0.56390521973 + 3 0.80540898341 0.23676109735 +D 4 + 1 56.088939191 0.29588869140E-01 + 2 15.751908917 0.158725714045 + 3 5.3115812367 0.379762291591 + 4 1.7737904917 0.468989591719 +D 1 + 1 0.51975583665 1.0000000 + +GALLIUM +S 6 + 1 87842.1262960 0.14271588798E-02 + 2 13189.4964020 0.10938894097E-01 + 3 3000.9482141 0.54308200127E-01 + 4 847.38425966 0.18951055983 + 5 276.00980642 0.38615185918 + 6 95.216672071 0.30051082719 +S 3 + 1 184.00703968 -0.11124461963 + 2 21.827051433 0.64831322918 + 3 9.0026963761 0.44358383594 +S 3 + 1 16.033176938 -0.23014299555 + 2 2.6707724878 0.72946861516 + 3 1.1252834170 0.46214976951 +S 1 + 1 0.15967608829 1.0000000 +S 1 + 1 0.057643044753 1.0000000 +P 5 + 1 1167.2665844 0.90974021307E-02 + 2 275.38062789 0.68456273343E-01 + 3 87.375073070 0.26922293888 + 4 31.597254670 0.53507897882 + 5 11.824122798 0.35279062760 +P 3 + 1 6.2881845133 0.33938244530 + 2 2.5199853240 0.56871947608 + 3 1.0169726306 0.27717694939 +P 1 + 1 0.24748373604 1.0000000 +P 1 + 1 0.66685294163E-01 1.0000000 +D 4 + 1 65.354237674 0.27370146382E-01 + 2 18.504656747 0.15099463976 + 3 6.3179620632 0.37485328545 + 4 2.1641389426 0.47510606851 +D 1 + 1 0.66693092754 1.0000000 + +GERMANIUM +S 6 + 1 93889.8366420 0.14233976060E-02 + 2 14097.4975280 0.10910795654E-01 + 3 3207.5477309 0.54183705943E-01 + 4 905.76727269 0.18922820349 + 5 295.11014693 0.38612847001 + 6 101.84713141 0.30164050736 +S 3 + 1 196.56719662 -0.11118770940 + 2 23.405292522 0.64616007369 + 3 9.6839116702 0.44188904568 +S 3 + 1 17.269736544 -0.23027421375 + 2 2.8964622160 0.73017169398 + 3 1.2553621412 0.46197222255 +S 1 + 1 0.20213081492 1.0000000 +S 1 + 1 0.73867910982E-01 1.0000000 +P 5 + 1 1259.2085995 0.90115464252E-02 + 2 297.15626382 0.67986841689E-01 + 3 94.353387522 0.26853856488 + 4 34.176329677 0.53659649219 + 5 12.816139615 0.35633514961 +P 3 + 1 6.8471029784 0.33900693119 + 2 2.7717363939 0.56809365264 + 3 1.1458418175 0.27246539884 +P 1 + 1 0.30679631536 1.0000000 +P 1 + 1 0.89283644115E-01 1.0000000 +D 4 + 1 74.782168177 0.25755860205E-01 + 2 21.310849759 0.14536816132 + 3 7.3464792363 0.37134209859 + 4 2.5656271395 0.48002998436 +D 1 + 1 0.81981773070 1.0000000 + +ARSENIC +S 6 + 1 100146.5255400 0.14258349617E-02 + 2 15036.8617110 0.10930176963E-01 + 3 3421.2902833 0.54294174610E-01 + 4 966.16965717 0.18976078153 + 5 314.87394026 0.38775195453 + 6 108.70823790 0.30402812040 +S 3 + 1 209.54238950 -0.11162094204 + 2 25.038221139 0.64697607762 + 3 10.390964343 0.44223608673 +S 3 + 1 18.555090093 -0.22994190569 + 2 3.1281217449 0.73319107613 + 3 1.3884885073 0.45533653943 +S 1 + 1 0.24714362141 1.0000000 +S 1 + 1 0.91429428670E-01 1.0000000 +P 5 + 1 1355.6443507 0.89182507898E-02 + 2 319.99929270 0.67454750717E-01 + 3 101.67734092 0.26759772110 + 4 36.886323845 0.53776844520 + 5 13.861115909 0.35992570244 +P 3 + 1 7.4260666912 0.34036849637 + 2 3.0316247187 0.57030149334 + 3 1.2783078340 0.26606170238 +P 1 + 1 0.37568503356 1.0000000 +P 1 + 1 0.11394805454 1.0000000 +D 4 + 1 84.445514539 0.24518402724E-01 + 2 24.190416102 0.14107454677 + 3 8.4045015119 0.36875228915 + 4 2.9808970748 0.48409561362 +D 1 + 1 0.97909243359 1.0000000 + +SELENIUM +S 6 + 1 106612.2002700 0.14274889113E-02 + 2 16007.6047010 0.10943525114E-01 + 3 3642.1699707 0.54374171596E-01 + 4 1028.5912993 0.19018092947 + 5 335.30298888 0.38913021696 + 6 115.80129154 0.30620207088 +S 3 + 1 222.93325020 -0.11198808009 + 2 26.726257934 0.64752124207 + 3 11.124501923 0.44241976651 +S 3 + 1 19.888520061 -0.22857227762 + 2 3.3668473803 0.73591359951 + 3 1.5249277839 0.44330199577 +S 1 + 1 0.29630037912 1.0000000 +S 1 + 1 0.11009288974 1.0000000 +P 5 + 1 1455.9068120 0.88203597043E-02 + 2 343.75101831 0.66875851967E-01 + 3 109.29554964 0.26640578080 + 4 39.707711022 0.53834928422 + 5 14.950185232 0.36303281993 +P 3 + 1 8.0208962094 0.34153807025 + 2 3.2934649756 0.57257906583 + 3 1.4058602438 0.25549813222 +P 1 + 1 0.45076123226 1.0000000 +P 1 + 1 0.13353413325 1.0000000 +D 4 + 1 94.494024044 0.23490101098D-01 + 2 27.188185260 0.13747735767D+00 + 3 9.5091567352 0.36649929124D+00 + 4 3.4170516853 0.48750989884D+00 +D 1 + 1 1.1479590083 1.0000000 + +BROMINE +S 6 + 1 113286.3877600 0.14283037779E-02 + 2 17009.6263030 0.10950417496E-01 + 3 3870.1842567 0.54421006604E-01 + 4 1093.0357227 0.19047907695 + 5 356.39721797 0.39024642737 + 6 123.12539643 0.30814432514 +S 3 + 1 236.74084007 -0.11228065671 + 2 28.468661070 0.64775962312 + 3 11.883443722 0.44235575986 +S 3 + 1 21.269633312 -0.22642576323 + 2 3.6129226841 0.73823712008 + 3 1.6626648969 0.42683868694 +S 1 + 1 0.34823793232 1.0000000 +S 1 + 1 0.13019031394 1.0000000 +P 5 + 1 1560.2801881 0.87166669072E-02 + 2 368.47859205 0.66243637420E-01 + 3 117.22978849 0.26495610385 + 4 42.648909248 0.53839160587 + 5 16.087225096 0.36579387888 +P 3 + 1 8.6352810058 0.34248787366 + 2 3.5613665502 0.57500678213 + 3 1.5292626609 0.24330394172 +P 1 + 1 0.53064294848 1.0000000 +P 1 + 1 0.15702758965 1.0000000 +D 4 + 1 104.85518642 0.22650147581D-01 + 2 30.281143688 0.13455483230 + 3 10.651394267 0.36474454537 + 4 3.8699456233 0.49044587056 +D 1 + 1 1.3240876762 1.0000000 + +KRYPTON +S 6 + 1 120165.6487500 0.0014294234496 + 2 18042.5001690 0.0109595970200 + 3 4105.1800388 0.0544796474700 + 4 1159.4447248 0.19081239122 + 5 378.13810346 0.39140445949 + 6 130.67610045 0.31007346689 +S 3 + 1 250.96373353 -0.11258604011 + 2 30.265930380 0.64816322524 + 3 12.668110757 0.44240238994 +S 3 + 1 22.697940101 -0.22391964908 + 2 3.8673017398 0.74060953026 + 3 1.8013706671 0.40803658445 +S 1 + 1 0.40308639136 1.0000000 +S 1 + 1 0.15170527749 1.0000000 +P 5 + 1 1668.5736623 0.86215009803E-02 + 2 394.13862798 0.65665506884E-01 + 3 125.46644854 0.26366576697 + 4 45.704992729 0.53867658574 + 5 17.270133229 0.36865167502 +P 3 + 1 9.2699487770 0.34346189015 + 2 3.8364617719 0.57784776196 + 3 1.6467163312 0.23075615197 +P 1 + 1 0.61479686489 1.0000000 +P 1 + 1 0.18319467414 1.0000000 +D 4 + 1 115.55868603 0.21945485333D-01 + 2 33.477985864 0.13211212879 + 3 11.834450207 0.36333306964 + 4 4.3408254894 0.49300484446 +D 1 + 1 1.5079273157 1.0000000 + diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index 7fd94ac4..2aafe644 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -6,7 +6,7 @@ type t = { element : Element.t ; charge : Charge.t ; coord : Point3d.t ; -} +} with sexp (** Read xyz coordinates of the atom with unit u *) let of_string u s = @@ -23,7 +23,7 @@ let of_string u s = | [ name; x; y; z ] -> let e = Element.of_string name in { element = e ; - charge = Charge.of_int (Element.to_charge e); + charge = Element.to_charge e; coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") } | _ -> raise (AtomError s) @@ -33,6 +33,6 @@ let to_string u a = [ Element.to_string a.element ; Charge.to_string a.charge ; Point3d.to_string u a.coord ] - |> String.concat ?sep:(Some " ") + |> String.concat ~sep:" " ;; diff --git a/ocaml/Atom.mli b/ocaml/Atom.mli new file mode 100644 index 00000000..93f7d885 --- /dev/null +++ b/ocaml/Atom.mli @@ -0,0 +1,9 @@ +exception AtomError of string + +type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; } + +val t_of_sexp : Sexplib.Sexp.t -> t +val sexp_of_t : t -> Sexplib.Sexp.t + +val of_string : Units.units -> string -> t +val to_string : Units.units -> t -> string diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml index 844269ed..5b8d9d14 100644 --- a/ocaml/Basis.ml +++ b/ocaml/Basis.ml @@ -1,7 +1,7 @@ open Core.Std;; open Qptypes;; -type t = (Gto.t * Nucl_number.t) list;; +type t = (Gto.t * Nucl_number.t) list with sexp;; (** Read all the basis functions of an element *) let read in_channel at_number = @@ -36,8 +36,28 @@ let read_element in_channel at_number element = ;; let to_string b = - List.map ~f:(fun (g,n) -> - let n = Nucl_number.to_int n in - (Int.to_string n)^":"^(Gto.to_string g)) b + let new_nucleus n = + Printf.sprintf "Atom %d" n + in + + let rec do_work accu current_nucleus = function + | [] -> List.rev accu + | (g,n)::tail -> + let n = Nucl_number.to_int n + in + let accu = + if (n <> current_nucleus) then + (new_nucleus n)::""::accu + else + accu + in + do_work ((Gto.to_string g)::accu) n tail + in + do_work [new_nucleus 1] 1 b |> String.concat ~sep:"\n" ;; + +include To_md5;; +let to_md5 = to_md5 sexp_of_t +;; + diff --git a/ocaml/Basis.mli b/ocaml/Basis.mli index 276578f4..4da99266 100644 --- a/ocaml/Basis.mli +++ b/ocaml/Basis.mli @@ -15,3 +15,6 @@ val read_element : (** Convert the basis to a string *) val to_string : (Gto.t * Nucl_number.t) list -> string + +(** Convert the basis to an MD5 hash *) +val to_md5 : (Gto.t * Nucl_number.t) list -> MD5.t diff --git a/ocaml/Bit.ml b/ocaml/Bit.ml index c9c8abdc..28b8dac9 100644 --- a/ocaml/Bit.ml +++ b/ocaml/Bit.ml @@ -1,3 +1,4 @@ +open Core.Std;; (* Type for bits @@ -7,9 +8,10 @@ Zero | One *) -type bit = - | One - | Zero +type t = +| One +| Zero +with sexp let to_string = function | Zero -> "0" diff --git a/ocaml/Bit.mli b/ocaml/Bit.mli new file mode 100644 index 00000000..6dd5abdd --- /dev/null +++ b/ocaml/Bit.mli @@ -0,0 +1,10 @@ +type t = One | Zero with sexp + +(** String conversions for printing *) +val to_string : t -> string + +(** Logical operations *) +val and_operator : t -> t -> t +val or_operator : t -> t -> t +val xor_operator : t -> t -> t +val not_operator : t -> t diff --git a/ocaml/Bitlist.ml b/ocaml/Bitlist.ml index d197c42b..80efd2b2 100644 --- a/ocaml/Bitlist.ml +++ b/ocaml/Bitlist.ml @@ -1,4 +1,5 @@ open Qptypes;; +open Core.Std;; (* Type for bits strings @@ -7,7 +8,7 @@ Type for bits strings list of Bits *) -type bit_list = Bit.bit list +type t = Bit.t list (* String representation *) let to_string b = @@ -20,6 +21,13 @@ let to_string b = do_work "" b ;; +let of_string ?(zero='0') ?(one='1') s = + String.to_list s + |> List.rev_map ~f:( fun c -> + if (c = zero) then Bit.Zero + else if (c = one) then Bit.One + else (failwith ("Error in string "^s) ) ) +;; (* Create a bit list from an int64 *) let of_int64 i = @@ -27,7 +35,7 @@ let of_int64 i = | 0L -> [ Bit.Zero ] | 1L -> [ Bit.One ] | i -> let b = - match (Int64.logand i 1L ) with + match (Int64.bit_and i 1L ) with | 0L -> Bit.Zero | 1L -> Bit.One | _ -> raise (Failure "i land 1 not in (0,1)") @@ -51,15 +59,14 @@ let to_int64 l = let rec do_work accu = function | [] -> accu | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail - | Bit.One::tail -> do_work Int64.(logor one (shift_left accu 1)) tail + | Bit.One::tail -> do_work Int64.(bit_or one (shift_left accu 1)) tail in do_work Int64.zero (List.rev l) ;; (* Create a bit list from a list of int64 *) let of_int64_list l = - let list_of_lists = List.map of_int64 l in - let result = List.rev list_of_lists in - List.flatten result + List.map ~f:of_int64 l + |> List.concat ;; (* Compute n_int *) @@ -92,27 +99,28 @@ let to_int64_list l = in let l = do_work [] [] 1 l in - List.rev_map to_int64 l + List.rev_map ~f:to_int64 l ;; (* Create a bit list from a list of MO indices *) let of_mo_number_list n_int l = let n_int = N_int_number.to_int n_int in let length = n_int*64 in - let a = Array.make length (Bit.Zero) in - List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l; + let a = Array.create length (Bit.Zero) in + List.iter ~f:(fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l; Array.to_list a ;; let to_mo_number_list l = let a = Array.of_list l in + let mo_tot_num = MO_number.get_max () in let rec do_work accu = function | 0 -> accu | i -> begin let new_accu = match a.(i-1) with - | Bit.One -> (MO_number.of_int i)::accu + | Bit.One -> (MO_number.of_int ~max:mo_tot_num i)::accu | Bit.Zero -> accu in do_work new_accu (i-1) @@ -161,30 +169,4 @@ let popcnt b = in popcnt 0 b ;; -let test_module () = - let test = of_int64_list ([-1231L;255L]) in - print_string (to_string test); - print_newline (); - print_string (string_of_int (String.length (to_string test))); - print_newline (); - print_string ( Bit.to_string Bit.One ); - - let a = of_int64_list ([-1L;0L]) - and b = of_int64_list ([128L;127L]) - in begin - print_newline (); - print_newline (); - print_string (to_string a); - print_newline (); - print_string (to_string b); - print_newline (); - print_string (to_string (and_operator a b)); - print_newline (); - print_string (to_string (or_operator a b)); - print_newline (); - print_string (to_string (xor_operator a b)); - print_string (to_string a); - print_int (popcnt a); - end -;; diff --git a/ocaml/Bitlist.mli b/ocaml/Bitlist.mli new file mode 100644 index 00000000..c733712c --- /dev/null +++ b/ocaml/Bitlist.mli @@ -0,0 +1,35 @@ +type t = Bit.t list + +(** The zero bit list *) +val zero : Qptypes.N_int_number.t -> t + +(** Convert to a string for printing *) +val to_string : t -> string + +(** Convert to a string for printing *) +val of_string : ?zero:char -> ?one:char -> string -> t + +(** int64 conversion functions *) + +val of_int64 : int64 -> t +val to_int64 : t -> int64 + +val of_int64_list : int64 list -> t +val to_int64_list : t -> int64 list + +(** Get the number of needed int64 elements to encode the bit list *) +val n_int_of_mo_tot_num : int -> Qptypes.N_int_number.t + +(** Conversion to MO numbers *) +val to_mo_number_list : t -> Qptypes.MO_number.t list +val of_mo_number_list : + Qptypes.N_int_number.t -> Qptypes.MO_number.t list -> t + +(** Logical operators *) +val and_operator : t -> t -> t +val xor_operator : t -> t -> t +val or_operator : t -> t -> t +val not_operator : t -> t + +(** Count the number of bits set to one *) +val popcnt : t -> int diff --git a/ocaml/Charge.ml b/ocaml/Charge.ml index 6b58065a..b40469d2 100644 --- a/ocaml/Charge.ml +++ b/ocaml/Charge.ml @@ -1,11 +1,12 @@ open Core.Std;; -type t = float +type t = float with sexp;; let of_float x = x let of_int i = Float.of_int i let of_string s = Float.of_string s + let to_float x = x let to_int x = Float.to_int x let to_string x = diff --git a/ocaml/Charge.mli b/ocaml/Charge.mli index 027c2fad..07685531 100644 --- a/ocaml/Charge.mli +++ b/ocaml/Charge.mli @@ -1,8 +1,13 @@ -type t = float +type t = float with sexp +(** Float conversion functions *) val to_float : t -> float -val to_int : t -> int -val to_string: t -> string val of_float : float -> t + +(** Int conversion functions *) +val to_int : t -> int val of_int : int -> t + +(** String conversion functions *) +val to_string: t -> string val of_string: string -> t diff --git a/ocaml/Determinant.ml b/ocaml/Determinant.ml new file mode 100644 index 00000000..96291904 --- /dev/null +++ b/ocaml/Determinant.ml @@ -0,0 +1,73 @@ +open Core.Std;; +open Qptypes;; + +type t = int64 array with sexp + +let to_int64_array (x:t) = (x:int64 array) +;; + +let to_alpha_beta x = + let x = to_int64_array x in + let n_int = (Array.length x)/2 in + ( Array.init n_int ~f:(fun i -> x.(i)) , + Array.init n_int ~f:(fun i -> x.(i+n_int)) ) +;; + +let to_bitlist_couple x = + let (xa,xb) = to_alpha_beta x in + let xa = to_int64_array xa + |> Array.to_list + |> Bitlist.of_int64_list + and xb = to_int64_array xb + |> Array.to_list + |> Bitlist.of_int64_list + in (xa,xb) +;; + +let bitlist_to_string ~mo_tot_num x = + List.map x ~f:(fun i -> match i with + | Bit.Zero -> "-" + | Bit.One -> "+" ) + |> String.concat + |> String.sub ~pos:0 ~len:(MO_number.to_int mo_tot_num) +;; + + +let of_int64_array ~n_int ~alpha ~beta x = + assert ((Array.length x) = (N_int_number.to_int n_int)*2) ; + let (a,b) = to_bitlist_couple x + and alpha = Elec_alpha_number.to_int alpha + and beta = Elec_beta_number.to_int beta + in + if ( (Bitlist.popcnt a) <> alpha) then + begin + let mo_tot_num = MO_number.get_max () in + let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in + failwith (Printf.sprintf "Expected %d electrons in alpha determinant +%s" alpha (bitlist_to_string ~mo_tot_num:mo_tot_num a) ) + end; + if ( (Bitlist.popcnt b) <> beta ) then + begin + let mo_tot_num = MO_number.get_max () in + let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in + failwith (Printf.sprintf "Expected %d electrons in beta determinant +%s" beta (bitlist_to_string ~mo_tot_num:mo_tot_num b) ) + end; + x +;; + +let of_bitlist_couple ~alpha ~beta (xa,xb) = + let ba = Bitlist.to_int64_list xa in + let bb = Bitlist.to_int64_list xb in + let n_int = Bitlist.n_int_of_mo_tot_num (List.length xa) in + of_int64_array ~n_int:n_int ~alpha:alpha ~beta:beta (Array.of_list (ba@bb)) +;; + +let to_string ~mo_tot_num x = + let (xa,xb) = to_bitlist_couple x in + [ bitlist_to_string ~mo_tot_num:mo_tot_num xa ; + bitlist_to_string ~mo_tot_num:mo_tot_num xb ] + |> String.concat ~sep:"\n" +;; + + diff --git a/ocaml/Determinant.mli b/ocaml/Determinant.mli new file mode 100644 index 00000000..f01e49d9 --- /dev/null +++ b/ocaml/Determinant.mli @@ -0,0 +1,32 @@ +(** Determinants are stored as follows : + * <-------- N_int ----------> + * [| i1_alpha ; i2_alpha ; ... ; + * i1_beta ; i2_beta ; ... ; |] + * where each int64 is a list of 64 MOs. When the bit is set + * to 1, the MO is occupied. + *) +type t = int64 array with sexp + +(** Transform to an int64 array *) +val to_int64_array : t -> int64 array + +(** Create from an int64 array, checking the number of alpha + * and beta electrons *) +val of_int64_array : n_int:Qptypes.N_int_number.t -> + alpha:Qptypes.Elec_alpha_number.t -> + beta:Qptypes.Elec_beta_number.t -> + int64 array -> t + +(** Split into an alpha-only and a beta-only determinant *) +val to_alpha_beta : t -> (int64 array)*(int64 array) + +(** Transform to a bit list *) +val to_bitlist_couple : t -> Bitlist.t * Bitlist.t + +(** Create from a bit list *) +val of_bitlist_couple : alpha:Qptypes.Elec_alpha_number.t -> + beta:Qptypes.Elec_beta_number.t -> + Bitlist.t * Bitlist.t -> t + +(** String representation *) +val to_string : mo_tot_num:Qptypes.MO_number.t -> t -> string diff --git a/ocaml/Element.ml b/ocaml/Element.ml index c6ab549e..d282340b 100644 --- a/ocaml/Element.ml +++ b/ocaml/Element.ml @@ -8,7 +8,7 @@ type t = |Li|Be |B |C |N |O |F |Ne |Na|Mg |Al|Si|P |S |Cl|Ar |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr -;; +with sexp;; let of_string x = match (String.capitalize (String.lowercase x)) with @@ -132,47 +132,49 @@ let to_long_string = function | Kr -> "Krypton" ;; -let to_charge = function -| X -> 0 -| H -> 1 -| He -> 2 -| Li -> 3 -| Be -> 4 -| B -> 5 -| C -> 6 -| N -> 7 -| O -> 8 -| F -> 9 -| Ne -> 10 -| Na -> 11 -| Mg -> 12 -| Al -> 13 -| Si -> 14 -| P -> 15 -| S -> 16 -| Cl -> 17 -| Ar -> 18 -| K -> 19 -| Ca -> 20 -| Sc -> 21 -| Ti -> 22 -| V -> 23 -| Cr -> 24 -| Mn -> 25 -| Fe -> 26 -| Co -> 27 -| Ni -> 28 -| Cu -> 29 -| Zn -> 30 -| Ga -> 31 -| Ge -> 32 -| As -> 33 -| Se -> 34 -| Br -> 35 -| Kr -> 36 +let to_charge c = + let result = match c with + | X -> 0 + | H -> 1 + | He -> 2 + | Li -> 3 + | Be -> 4 + | B -> 5 + | C -> 6 + | N -> 7 + | O -> 8 + | F -> 9 + | Ne -> 10 + | Na -> 11 + | Mg -> 12 + | Al -> 13 + | Si -> 14 + | P -> 15 + | S -> 16 + | Cl -> 17 + | Ar -> 18 + | K -> 19 + | Ca -> 20 + | Sc -> 21 + | Ti -> 22 + | V -> 23 + | Cr -> 24 + | Mn -> 25 + | Fe -> 26 + | Co -> 27 + | Ni -> 28 + | Cu -> 29 + | Zn -> 30 + | Ga -> 31 + | Ge -> 32 + | As -> 33 + | Se -> 34 + | Br -> 35 + | Kr -> 36 + in Charge.of_int result ;; -let of_charge = function +let of_charge c = match (Charge.to_int c) with | 0 -> X | 1 -> H | 2 -> He diff --git a/ocaml/Element.mli b/ocaml/Element.mli new file mode 100644 index 00000000..48bd3c61 --- /dev/null +++ b/ocaml/Element.mli @@ -0,0 +1,18 @@ +exception ElementError of string + +type t = +|X +|H |He +|Li|Be |B |C |N |O |F |Ne +|Na|Mg |Al|Si|P |S |Cl|Ar +|K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr +with sexp + +(** String conversion functions *) +val of_string : string -> t +val to_string : t -> string +val to_long_string : t -> string + +(** get the positive charge *) +val to_charge : t -> Charge.t +val of_charge : Charge.t -> t diff --git a/ocaml/Excitation.ml b/ocaml/Excitation.ml index 89d76f1f..d58c3093 100644 --- a/ocaml/Excitation.ml +++ b/ocaml/Excitation.ml @@ -1,22 +1,14 @@ open Core.Std;; open Qptypes;; -module Hole : sig - type t - val to_mo_class : t -> MO_class.t - val of_mo_class : MO_class.t -> t -end = struct - type t = MO_class.t +module Hole = struct + type t = MO_class.t with sexp let of_mo_class x = x let to_mo_class x = x end -module Particle : sig - type t - val to_mo_class : t -> MO_class.t - val of_mo_class : MO_class.t -> t -end = struct - type t = MO_class.t +module Particle = struct + type t = MO_class.t with sexp let of_mo_class x = x let to_mo_class x = x end @@ -24,10 +16,7 @@ end type t = | Single of Hole.t*Particle.t | Double of Hole.t*Particle.t*Hole.t*Particle.t -;; - -let failwith s = raise (Failure s) -;; +with sexp;; let create_single ~hole ~particle = MO_class.( diff --git a/ocaml/Excitation.mli b/ocaml/Excitation.mli new file mode 100644 index 00000000..982cfd0e --- /dev/null +++ b/ocaml/Excitation.mli @@ -0,0 +1,30 @@ +module Hole : + sig + type t + val to_mo_class : t -> MO_class.t + val of_mo_class : MO_class.t -> t + val t_of_sexp : Sexplib.Sexp.t -> t + val sexp_of_t : t -> Sexplib.Sexp.t + end +module Particle : + sig + type t + val to_mo_class : t -> MO_class.t + val of_mo_class : MO_class.t -> t + val t_of_sexp : Sexplib.Sexp.t -> t + val sexp_of_t : t -> Sexplib.Sexp.t + end + +type t = + | Single of Hole.t * Particle.t + | Double of Hole.t * Particle.t * Hole.t * Particle.t +with sexp + +val create_single : hole:MO_class.t -> particle:MO_class.t -> t + +val double_of_singles : t -> t -> t + +val create_double : hole1:MO_class.t -> particle1:MO_class.t -> + hole2:MO_class.t -> particle2:MO_class.t -> t + +val to_string : t -> string diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index fc1953b2..de9de237 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -7,7 +7,7 @@ exception End_Of_Basis type t = { sym : Symmetry.t ; lc : ((Primitive.t * AO_coef.t) list) -} +} with sexp ;; let of_prim_coef_list pc = @@ -52,24 +52,37 @@ let read_one in_channel = match buffer with | [ j ; expo ; coef ] -> begin - let p = { Primitive.sym = sym ; - Primitive.expo = AO_expo.of_float - (Float.of_string expo) - } + let p = + Primitive.of_sym_expo sym + (AO_expo.of_float (Float.of_string expo) ) and c = AO_coef.of_float (Float.of_string coef) in read_lines ( (p,c)::result) (i-1) end | _ -> raise (GTO_Read_Failure line_buffer) end in read_lines [] n + |> List.rev |> of_prim_coef_list ;; (** Transform the gto to a string *) let to_string { sym = sym ; lc = lc } = - let f (p,c) = Printf.sprintf "( %s, %f )" (Primitive.to_string p) (AO_coef.to_float c) + let result = + Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc) in - Printf.sprintf "[ %s : %s ]" (Symmetry.to_string sym) - (String.concat (List.map ~f:f lc) ~sep:", ") + let rec do_work accu i = function + | [] -> List.rev accu + | (p,c)::tail -> + let p = AO_expo.to_float p.Primitive.expo + and c = AO_coef.to_float c + in + let result = + Printf.sprintf "%3d %16f %16f" i p c + in + do_work (result::accu) (i+1) tail + in + (do_work [result] 1 lc) + |> String.concat ~sep:"\n" ;; + diff --git a/ocaml/Gto.mli b/ocaml/Gto.mli new file mode 100644 index 00000000..fad133a3 --- /dev/null +++ b/ocaml/Gto.mli @@ -0,0 +1,16 @@ +exception GTO_Read_Failure of string +exception End_Of_Basis +type t = + { sym : Symmetry.t ; + lc : (Primitive.t * Qptypes.AO_coef.t) list; + } with sexp + +(** Create from a list of Primitive.t * Qptypes.AO_coef.t *) +val of_prim_coef_list : + (Primitive.t * Qptypes.AO_coef.t) list -> t + +(** Read from a file *) +val read_one : in_channel -> t + +(** Convert to string for printing *) +val to_string : t -> string diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index 276e4b8d..0e61dc54 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -12,10 +12,12 @@ module Ao_basis : sig ao_power : Symmetry.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; - } + } with sexp ;; val read : unit -> t val to_string : t -> string + val to_md5 : t -> MD5.t + val to_rst : t -> Rst_string.t end = struct type t = { ao_basis : string ; @@ -26,7 +28,7 @@ end = struct ao_power : Symmetry.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "ao_basis";; @@ -44,28 +46,29 @@ end = struct ;; let read_ao_prim_num () = - (Ezfio.get_ao_basis_ao_prim_num () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio |> Array.map ~f:AO_prim_number.of_int ;; let read_ao_prim_num_max () = - (Ezfio.get_ao_basis_ao_prim_num () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_ao_basis_ao_prim_num () + |> Ezfio.flattened_ezfio |> Array.fold ~f:(fun x y -> if x>y then x else y) ~init:0 |> AO_prim_number.of_int ;; let read_ao_nucl () = - (Ezfio.get_ao_basis_ao_nucl () ).Ezfio.data - |> Ezfio.flattened_ezfio_data - |> Array.map ~f:Nucl_number.of_int + let nmax = Nucl_number.get_max () in + Ezfio.get_ao_basis_ao_nucl () + |> Ezfio.flattened_ezfio + |> Array.map ~f:(fun x-> Nucl_number.of_int ~max:nmax x) ;; let read_ao_power () = let x = Ezfio.get_ao_basis_ao_power () in let dim = x.Ezfio.dim.(0) in - let data = Ezfio.flattened_ezfio_data x.Ezfio.data in + let data = Ezfio.flattened_ezfio x in let result = Array.init dim ~f:(fun x -> "") in for i=1 to dim do @@ -80,14 +83,14 @@ end = struct ;; let read_ao_coef () = - (Ezfio.get_ao_basis_ao_coef () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_ao_basis_ao_coef () + |> Ezfio.flattened_ezfio |> Array.map ~f:AO_coef.of_float ;; let read_ao_expo () = - (Ezfio.get_ao_basis_ao_expo () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_ao_basis_ao_expo () + |> Ezfio.flattened_ezfio |> Array.map ~f:AO_expo.of_float ;; @@ -102,10 +105,105 @@ end = struct ao_expo = read_ao_expo () ; } ;; + + let to_long_basis b = + let ao_num = AO_number.to_int b.ao_num in + let gto_array = Array.init (AO_number.to_int b.ao_num) + ~f:(fun i -> + let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in + let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in + let prims = List.init ao_prim_num ~f:(fun j -> + let prim = { Primitive.sym = s ; + Primitive.expo = b.ao_expo.(ao_num*j+i) + } + in + let coef = b.ao_coef.(ao_num*j+i) in + (prim,coef) + ) in + Gto.of_prim_coef_list prims + ) + in + let rec do_work accu sym gto nucl = + match (sym, gto, nucl) with + | (s::srest, g::grest, n::nrest) -> + do_work ((s,g,n)::accu) srest grest nrest + | ([],[],[]) -> List.rev accu + | _ -> assert false + in + do_work [] + (Array.to_list b.ao_power) + (Array.to_list gto_array) + (Array.to_list b.ao_nucl) + ;; + let to_basis b = + to_long_basis b + |> Long_basis.to_basis + ;; + + let to_rst b = + let print_sym = + let l = List.init (Array.length b.ao_power) ~f:( + fun i -> ( (i+1),b.ao_nucl.(i),b.ao_power.(i) ) ) in + let rec do_work = function + | [] -> [] + | (i,n,x)::tail -> + (Printf.sprintf " %5d %6d %-8s\n" i (Nucl_number.to_int n) (Symmetry.Xyz.to_string x)):: + (do_work tail) + in do_work l + |> String.concat + in + + let short_basis = to_basis b in + Printf.sprintf " +Name of the AO basis :: + + ao_basis = %s + +Basis set :: + +%s + + +======= ========= =========== + Basis Nucleus Symmetries +======= ========= =========== +%s +======= ========= =========== + +" b.ao_basis + (Basis.to_string short_basis + |> String.split ~on:'\n' + |> List.map ~f:(fun x-> " "^x) + |> String.concat ~sep:"\n" + ) print_sym + |> Rst_string.of_string + ;; + + let read_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + in + let rec extract_basis = function + | [] -> failwith "Error in basis set" + | line :: tail -> + let line = String.strip line in + if line = "Basis set ::" then + String.concat tail ~sep:"\n" + else + extract_basis tail + in + extract_basis s + ;; + + let to_md5 b = + let short_basis = to_basis b in + Basis.to_md5 short_basis + ;; + let to_string b = Printf.sprintf " -ao_basis = \"%s\" +ao_basis = %s ao_num = %s ao_prim_num = %s ao_prim_num_max = %s @@ -113,6 +211,7 @@ ao_nucl = %s ao_power = %s ao_coef = %s ao_expo = %s +md5 = %s " b.ao_basis (AO_number.to_string b.ao_num) @@ -127,7 +226,8 @@ ao_expo = %s |> String.concat ~sep:", ") (b.ao_expo |> Array.to_list |> List.map ~f:AO_expo.to_string |> String.concat ~sep:", ") + (to_md5 b |> MD5.to_string ) + ;; end - diff --git a/ocaml/Input_bi_integrals.ml b/ocaml/Input_bi_integrals.ml index 3322f4c2..8cc8203c 100644 --- a/ocaml/Input_bi_integrals.ml +++ b/ocaml/Input_bi_integrals.ml @@ -11,10 +11,13 @@ module Bielec_integrals : sig threshold_ao : Threshold.t; threshold_mo : Threshold.t; direct : bool; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t -> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = { read_ao_integrals : bool; @@ -24,7 +27,7 @@ end = struct threshold_ao : Threshold.t; threshold_mo : Threshold.t; direct : bool; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "bielec_integrals";; @@ -38,6 +41,11 @@ end = struct Ezfio.get_bielec_integrals_read_ao_integrals () ;; + let write_read_ao_integrals = + Ezfio.set_bielec_integrals_read_ao_integrals + ;; + + let read_read_mo_integrals () = if not (Ezfio.has_bielec_integrals_read_mo_integrals ()) then get_default "read_mo_integrals" @@ -47,6 +55,11 @@ end = struct Ezfio.get_bielec_integrals_read_mo_integrals () ;; + let write_read_mo_integrals = + Ezfio.set_bielec_integrals_read_mo_integrals + ;; + + let read_write_ao_integrals () = if not (Ezfio.has_bielec_integrals_write_ao_integrals ()) then get_default "write_ao_integrals" @@ -56,6 +69,11 @@ end = struct Ezfio.get_bielec_integrals_write_ao_integrals () ;; + let write_write_ao_integrals = + Ezfio.set_bielec_integrals_write_ao_integrals + ;; + + let read_write_mo_integrals () = if not (Ezfio.has_bielec_integrals_write_mo_integrals ()) then get_default "write_mo_integrals" @@ -65,6 +83,11 @@ end = struct Ezfio.get_bielec_integrals_write_mo_integrals () ;; + let write_write_mo_integrals = + Ezfio.set_bielec_integrals_write_mo_integrals + ;; + + let read_direct () = if not (Ezfio.has_bielec_integrals_direct ()) then get_default "direct" @@ -74,6 +97,11 @@ end = struct Ezfio.get_bielec_integrals_direct () ;; + let write_direct = + Ezfio.set_bielec_integrals_direct + ;; + + let read_threshold_ao () = if not (Ezfio.has_bielec_integrals_threshold_ao ()) then get_default "threshold_ao" @@ -84,6 +112,12 @@ end = struct |> Threshold.of_float ;; + let write_threshold_ao t = + Threshold.to_float t + |> Ezfio.set_bielec_integrals_threshold_ao + ;; + + let read_threshold_mo () = if not (Ezfio.has_bielec_integrals_threshold_mo ()) then get_default "threshold_mo" @@ -94,6 +128,12 @@ end = struct |> Threshold.of_float ;; + let write_threshold_mo t = + Threshold.to_float t + |> Ezfio.set_bielec_integrals_threshold_mo + ;; + + let read ()= let result = { read_ao_integrals = read_read_ao_integrals(); @@ -113,15 +153,31 @@ end = struct result ;; + let write b = + if (b.read_ao_integrals && + b.write_ao_integrals) then + failwith "Read and Write AO integrals are both true."; + if (b.read_mo_integrals && + b.write_mo_integrals) then + failwith "Read and Write MO integrals are both true."; + write_read_ao_integrals b.read_ao_integrals; + write_read_mo_integrals b.read_mo_integrals; + write_write_ao_integrals b.write_ao_integrals ; + write_write_mo_integrals b.write_mo_integrals ; + write_threshold_ao b.threshold_ao; + write_threshold_mo b.threshold_mo; + write_direct b.direct; + ;; + let to_string b = Printf.sprintf " -read_ao_integrals = %s -read_mo_integrals = %s +read_ao_integrals = %s +read_mo_integrals = %s write_ao_integrals = %s write_mo_integrals = %s -threshold_ao = %s -threshold_mo = %s -direct = %s +threshold_ao = %s +threshold_mo = %s +direct = %s " (Bool.to_string b.read_ao_integrals) (Bool.to_string b.read_mo_integrals) @@ -130,6 +186,55 @@ direct = %s (Threshold.to_string b.threshold_ao) (Threshold.to_string b.threshold_mo) (Bool.to_string b.direct) + ;; + + let to_rst b = + Printf.sprintf " +Read AO/MO integrals from disk :: + + read_ao_integrals = %s + read_mo_integrals = %s + +Write AO/MO integrals to disk :: + + write_ao_integrals = %s + write_mo_integrals = %s + +Thresholds on integrals :: + + threshold_ao = %s + threshold_mo = %s + +Direct calculation of integrals :: + + direct = %s + +" + (Bool.to_string b.read_ao_integrals) + (Bool.to_string b.read_mo_integrals) + (Bool.to_string b.write_ao_integrals) + (Bool.to_string b.write_mo_integrals) + (Threshold.to_string b.threshold_ao) + (Threshold.to_string b.threshold_mo) + (Bool.to_string b.direct) + |> Rst_string.of_string + ;; + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + end diff --git a/ocaml/Input_bitmasks.ml b/ocaml/Input_bitmasks.ml index be82f603..4b02a735 100644 --- a/ocaml/Input_bitmasks.ml +++ b/ocaml/Input_bitmasks.ml @@ -8,7 +8,7 @@ module Bitmasks : sig bit_kind : Bit_kind.t; n_mask_gen : Bitmask_number.t; generators : int64 array; - } + } with sexp ;; val read : unit -> t val to_string : t -> string @@ -18,7 +18,7 @@ end = struct bit_kind : Bit_kind.t; n_mask_gen : Bitmask_number.t; generators : int64 array; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "bitmasks";; @@ -72,8 +72,8 @@ end = struct in Ezfio.set_bitmasks_generators generators end; - (Ezfio.get_bitmasks_generators ()).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_bitmasks_generators () + |> Ezfio.flattened_ezfio ;; let read () = diff --git a/ocaml/Input_cis.ml b/ocaml/Input_cis.ml index 6b521872..79cf03e8 100644 --- a/ocaml/Input_cis.ml +++ b/ocaml/Input_cis.ml @@ -10,10 +10,11 @@ module Cis_dressed : sig mp2_dressing : bool; standard_doubles : bool; en_2_2 : bool; - } + } with sexp ;; val read : unit -> t val to_string : t -> string + val to_rst : t -> Rst_string.t end = struct type t = { n_state_cis : States_number.t; @@ -22,7 +23,7 @@ end = struct mp2_dressing : bool; standard_doubles : bool; en_2_2 : bool; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "cis_dressed";; @@ -95,9 +96,9 @@ end = struct let to_string b = Printf.sprintf " -n_state_cis = %s -n_core_cis = %s -n_act_cis = %s +n_state_cis = %s +n_core_cis = %s +n_act_cis = %s mp2_dressing = %s standard_doubles = %s en_2_2 = %s @@ -108,6 +109,42 @@ en_2_2 = %s (Bool.to_string b.mp2_dressing) (Bool.to_string b.standard_doubles) (Bool.to_string b.en_2_2) + ;; + + + let to_rst b = + Printf.sprintf " +Number of states :: + + n_state_cis = %s + +Core and active MOs :: + + n_core_cis = %s + n_act_cis = %s + +Dress with MP2 perturbation :: + + mp2_dressing = %s + +Use standard double-excitations :: + + standard_doubles = %s + +Epstein-Nesbet 2x2 diagonalization :: + + en_2_2 = %s + +" + (States_number.to_string b.n_state_cis) + (Positive_int.to_string b.n_core_cis) + (Positive_int.to_string b.n_act_cis) + (Bool.to_string b.mp2_dressing) + (Bool.to_string b.standard_doubles) + (Bool.to_string b.en_2_2) + + |> Rst_string.of_string + ;; end diff --git a/ocaml/Input_cisd_sc2.ml b/ocaml/Input_cisd_sc2.ml index 4a5461f4..47dca5e5 100644 --- a/ocaml/Input_cisd_sc2.ml +++ b/ocaml/Input_cisd_sc2.ml @@ -4,19 +4,22 @@ open Core.Std;; module Cisd_sc2 : sig type t = - { n_det_max_cisd_sc2 : Det_number.t; + { n_det_max_cisd_sc2 : Det_number_max.t; pt2_max : PT2_energy.t; do_pt2_end : bool; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t -> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = - { n_det_max_cisd_sc2 : Det_number.t; + { n_det_max_cisd_sc2 : Det_number_max.t; pt2_max : PT2_energy.t; do_pt2_end : bool; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "cisd_sc2_selected";; @@ -28,7 +31,12 @@ end = struct |> Ezfio.set_cisd_sc2_selected_n_det_max_cisd_sc2 ; Ezfio.get_cisd_sc2_selected_n_det_max_cisd_sc2 () - |> Det_number.of_int + |> Det_number_max.of_int + ;; + + let write_n_det_max_cisd_sc2 n = + Det_number_max.to_int n + |> Ezfio.set_cisd_sc2_selected_n_det_max_cisd_sc2 ;; @@ -42,6 +50,11 @@ end = struct |> PT2_energy.of_float ;; + let write_pt2_max p = + PT2_energy.to_float p + |> Ezfio.set_cisd_sc2_selected_pt2_max + ;; + let read_do_pt2_end () = if not (Ezfio.has_cisd_sc2_selected_do_pt2_end ()) then @@ -52,6 +65,10 @@ end = struct Ezfio.get_cisd_sc2_selected_do_pt2_end () ;; + let write_do_pt2_end = + Ezfio.set_cisd_sc2_selected_do_pt2_end + ;; + let read () = { n_det_max_cisd_sc2 = read_n_det_max_cisd_sc2 (); @@ -60,15 +77,64 @@ end = struct } ;; + let write { n_det_max_cisd_sc2 ; + pt2_max ; + do_pt2_end ; + } = + write_n_det_max_cisd_sc2 n_det_max_cisd_sc2; + write_pt2_max pt2_max; + write_do_pt2_end do_pt2_end; + ;; + let to_string b = Printf.sprintf " n_det_max_cisd_sc2 = %s -pt2_max = %s -do_pt2_end = %s +pt2_max = %s +do_pt2_end = %s " - (Det_number.to_string b.n_det_max_cisd_sc2) + (Det_number_max.to_string b.n_det_max_cisd_sc2) (PT2_energy.to_string b.pt2_max) (Bool.to_string b.do_pt2_end) + ;; + + let to_rst b = + Printf.sprintf " +Stop when the `n_det` > `n_det_max_cisd_sc2` :: + + n_det_max_cisd_sc2 = %s + +Stop when -E(PT2) < `pt2_max` :: + + pt2_max = %s + +Compute E(PT2) at the end :: + + do_pt2_end = %s + +" + (Det_number_max.to_string b.n_det_max_cisd_sc2) + (PT2_energy.to_string b.pt2_max) + (Bool.to_string b.do_pt2_end) + |> Rst_string.of_string + ;; + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + + + end diff --git a/ocaml/Input_determinants.ml b/ocaml/Input_determinants.ml index fb2f599f..f4eeee6a 100644 --- a/ocaml/Input_determinants.ml +++ b/ocaml/Input_determinants.ml @@ -6,11 +6,11 @@ module Determinants : sig type t = { n_int : N_int_number.t; bit_kind : Bit_kind.t; - mo_label : Non_empty_string.t; + mo_label : MO_label.t; n_det : Det_number.t; n_states : States_number.t; n_states_diag : States_number.t; - n_det_max_jacobi : Det_number.t; + n_det_max_jacobi : Strictly_positive_int.t; threshold_generators : Threshold.t; threshold_selectors : Threshold.t; read_wf : bool; @@ -18,19 +18,21 @@ module Determinants : sig s2_eig : bool; psi_coef : Det_coef.t array; psi_det : Determinant.t array; - } - ;; - val read : unit -> t + } with sexp + val read : unit -> t + val write : t -> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = { n_int : N_int_number.t; bit_kind : Bit_kind.t; - mo_label : Non_empty_string.t; + mo_label : MO_label.t; n_det : Det_number.t; n_states : States_number.t; n_states_diag : States_number.t; - n_det_max_jacobi : Det_number.t; + n_det_max_jacobi : Strictly_positive_int.t; threshold_generators : Threshold.t; threshold_selectors : Threshold.t; read_wf : bool; @@ -38,7 +40,7 @@ end = struct s2_eig : bool; psi_coef : Det_coef.t array; psi_det : Determinant.t array; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "determinants";; @@ -54,6 +56,12 @@ end = struct |> N_int_number.of_int ;; + let write_n_int n = + N_int_number.to_int n + |> Ezfio.set_determinants_n_int + ;; + + let read_bit_kind () = if not (Ezfio.has_determinants_bit_kind ()) then Lazy.force Qpackage.bit_kind @@ -64,15 +72,27 @@ end = struct |> Bit_kind.of_int ;; + let write_bit_kind b = + Bit_kind.to_int b + |> Ezfio.set_determinants_bit_kind + ;; + + let read_mo_label () = if (not (Ezfio.has_determinants_mo_label ())) then Ezfio.get_mo_basis_mo_label () |> Ezfio.set_determinants_mo_label ; Ezfio.get_determinants_mo_label () - |> Non_empty_string.of_string + |> MO_label.of_string ;; + let write_mo_label l = + MO_label.to_string l + |> Ezfio.set_determinants_mo_label + ;; + + let read_n_det () = if not (Ezfio.has_determinants_n_det ()) then Ezfio.set_determinants_n_det 1 @@ -81,6 +101,12 @@ end = struct |> Det_number.of_int ;; + let write_n_det n = + Det_number.to_int n + |> Ezfio.set_determinants_n_det + ;; + + let read_n_states () = if not (Ezfio.has_determinants_n_states ()) then Ezfio.set_determinants_n_states 1 @@ -89,6 +115,12 @@ end = struct |> States_number.of_int ;; + let write_n_states n = + States_number.to_int n + |> Ezfio.set_determinants_n_states + ;; + + let read_n_states_diag () = if not (Ezfio.has_determinants_n_states_diag ()) then read_n_states () @@ -99,6 +131,12 @@ end = struct |> States_number.of_int ;; + let write_n_states_diag n = + States_number.to_int n + |> Ezfio.set_determinants_n_states_diag + ;; + + let read_n_det_max_jacobi () = if not (Ezfio.has_determinants_n_det_max_jacobi ()) then get_default "n_det_max_jacobi" @@ -106,9 +144,15 @@ end = struct |> Ezfio.set_determinants_n_det_max_jacobi ; Ezfio.get_determinants_n_det_max_jacobi () - |> Det_number.of_int + |> Strictly_positive_int.of_int ;; + let write_n_det_max_jacobi n = + Strictly_positive_int.to_int n + |> Ezfio.set_determinants_n_det_max_jacobi + ;; + + let read_threshold_generators () = if not (Ezfio.has_determinants_threshold_generators ()) then get_default "threshold_generators" @@ -119,6 +163,12 @@ end = struct |> Threshold.of_float ;; + let write_threshold_generators t = + Threshold.to_float t + |> Ezfio.set_determinants_threshold_generators + ;; + + let read_threshold_selectors () = if not (Ezfio.has_determinants_threshold_selectors ()) then get_default "threshold_selectors" @@ -129,6 +179,12 @@ end = struct |> Threshold.of_float ;; + let write_threshold_selectors t = + Threshold.to_float t + |> Ezfio.set_determinants_threshold_selectors + ;; + + let read_read_wf () = if not (Ezfio.has_determinants_read_wf ()) then get_default "read_wf" @@ -138,6 +194,9 @@ end = struct Ezfio.get_determinants_read_wf () ;; + let write_read_wf = Ezfio.set_determinants_read_wf ;; + + let read_expected_s2 () = if not (Ezfio.has_determinants_expected_s2 ()) then begin @@ -153,6 +212,12 @@ end = struct |> Positive_float.of_float ;; + let write_expected_s2 s2 = + Positive_float.to_float s2 + |> Ezfio.set_determinants_expected_s2 + ;; + + let read_s2_eig () = if not (Ezfio.has_determinants_s2_eig ()) then get_default "s2_eig" @@ -162,59 +227,80 @@ end = struct Ezfio.get_determinants_s2_eig () ;; + let write_s2_eig = Ezfio.set_determinants_s2_eig ;; + + let read_psi_coef () = if not (Ezfio.has_determinants_psi_coef ()) then Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| 1 |] ~data:[1.] |> Ezfio.set_determinants_psi_coef ; - (Ezfio.get_determinants_psi_coef ()).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_determinants_psi_coef () + |> Ezfio.flattened_ezfio |> Array.map ~f:Det_coef.of_float ;; + let write_psi_coef ~n_det c = + let n_det = Det_number.to_int n_det + and c = Array.to_list c + |> List.map ~f:Det_coef.to_float + in + Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_det |] ~data:c + |> Ezfio.set_determinants_psi_coef + ;; + + let read_psi_det () = - let n_int = read_n_int () in + let n_int = read_n_int () + and n_alpha = Ezfio.get_electrons_elec_alpha_num () + |> Elec_alpha_number.of_int + and n_beta = Ezfio.get_electrons_elec_beta_num () + |> Elec_beta_number.of_int + in if not (Ezfio.has_determinants_psi_det ()) then begin + let mo_tot_num = MO_number.get_max () in let rec build_data accu = function | 0 -> accu - | n -> build_data ((MO_number.of_int n)::accu) (n-1) + | n -> build_data ((MO_number.of_int ~max:mo_tot_num n)::accu) (n-1) in - let det_a = build_data [] (Ezfio.get_electrons_elec_alpha_num ()) + let det_a = build_data [] (Elec_alpha_number.to_int n_alpha) |> Bitlist.of_mo_number_list n_int - and det_b = build_data [] (Ezfio.get_electrons_elec_beta_num ()) + and det_b = build_data [] (Elec_beta_number.to_int n_beta) |> Bitlist.of_mo_number_list n_int in let data = ( (Bitlist.to_int64_list det_a) @ (Bitlist.to_int64_list det_b) ) in Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; 1 |] ~data:data - |> Ezfio.set_determinants_psi_det + |> Ezfio.set_determinants_psi_det ; end ; let n_int = N_int_number.to_int n_int in - let rec transform accu1 accu2 n_rest = function - | [] -> - let accu1 = List.rev accu1 - |> Array.of_list - |> Determinant.of_int64_array - in - List.rev (accu1::accu2) |> Array.of_list - | i::rest -> - if (n_rest > 0) then - transform (i::accu1) accu2 (n_rest-1) rest - else - let accu1 = List.rev accu1 - |> Array.of_list - |> Determinant.of_int64_array - in - transform [] (accu1::accu2) (2*n_int) rest + let psi_det_array = Ezfio.get_determinants_psi_det () in + let dim = psi_det_array.Ezfio.dim + and data = Ezfio.flattened_ezfio psi_det_array in - (Ezfio.get_determinants_psi_det ()).Ezfio.data - |> Ezfio.flattened_ezfio_data - |> Array.to_list - |> transform [] [] (2*n_int) + assert (n_int = dim.(0)); + assert (dim.(1) = 2); + assert (dim.(2) = (Det_number.to_int (read_n_det ()))); + List.init dim.(2) ~f:(fun i -> + Array.sub ~pos:(2*n_int*i) ~len:(2*n_int) data) + |> List.map ~f:(Determinant.of_int64_array + ~n_int:(N_int_number.of_int n_int) + ~alpha:n_alpha ~beta:n_beta ) + |> Array.of_list ;; + let write_psi_det ~n_int ~n_det d = + let data = Array.to_list d + |> Array.concat + |> Array.to_list + in + Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data + |> Ezfio.set_determinants_psi_det + ;; + + let read () = { n_int = read_n_int () ; bit_kind = read_bit_kind () ; @@ -233,7 +319,108 @@ end = struct } ;; + let write { n_int ; + bit_kind ; + mo_label ; + n_det ; + n_states ; + n_states_diag ; + n_det_max_jacobi ; + threshold_generators ; + threshold_selectors ; + read_wf ; + expected_s2 ; + s2_eig ; + psi_coef ; + psi_det ; + } = + write_n_int n_int ; + write_bit_kind bit_kind; + write_mo_label mo_label; + write_n_det n_det; + write_n_states n_states; + write_n_states_diag n_states_diag; + write_n_det_max_jacobi n_det_max_jacobi; + write_threshold_generators threshold_generators; + write_threshold_selectors threshold_selectors; + write_read_wf read_wf; + write_expected_s2 expected_s2; + write_s2_eig s2_eig; + write_psi_coef ~n_det:n_det psi_coef; + write_psi_det ~n_int:n_int ~n_det:n_det psi_det; + ;; + + + let to_rst b = + let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in + let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in + let det_text = + List.map2_exn ~f:(fun coef det -> + Printf.sprintf " %F\n%s\n" + (Det_coef.to_float coef) + (Determinant.to_string ~mo_tot_num:mo_tot_num det + |> String.split ~on:'\n' + |> List.map ~f:(fun x -> " "^x) + |> String.concat ~sep:"\n" + ) + ) (Array.to_list b.psi_coef) (Array.to_list b.psi_det) + |> String.concat ~sep:"\n" + in + Printf.sprintf " +Read the current wave function :: + + read_wf = %s + +Label of the MOs on which the determinants were computed :: + + mo_label = %s + +Force the selected wave function to be an eigenfunction of S^2. +If true, input the expected value of S^2 :: + + s2_eig = %s + expected_s2 = %s + +Thresholds on generators and selectors (fraction of the norm) :: + + threshold_generators = %s + threshold_selectors = %s + +Number of requested states, and number of states used for the +Davidson diagonalization :: + + n_states = %s + n_states_diag = %s + +Maximum size of the Hamiltonian matrix that will be fully diagonalized :: + + n_det_max_jacobi = %s + +Number of determinants :: + + n_det = %s + +Determinants :: + +%s +" + (b.read_wf |> Bool.to_string) + (b.mo_label |> MO_label.to_string) + (b.s2_eig |> Bool.to_string) + (b.expected_s2 |> Positive_float.to_string) + (b.threshold_generators |> Threshold.to_string) + (b.threshold_selectors |> Threshold.to_string) + (b.n_states |> States_number.to_string) + (b.n_states_diag |> States_number.to_string) + (b.n_det_max_jacobi |> Strictly_positive_int.to_string) + (b.n_det |> Det_number.to_string) + det_text + |> Rst_string.of_string + ;; + let to_string b = + let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in + let mo_tot_num = MO_number.of_int mo_tot_num ~max:mo_tot_num in Printf.sprintf " n_int = %s bit_kind = %s @@ -252,11 +439,11 @@ psi_det = %s " (b.n_int |> N_int_number.to_string) (b.bit_kind |> Bit_kind.to_string) - (b.mo_label |> Non_empty_string.to_string) + (b.mo_label |> MO_label.to_string) (b.n_det |> Det_number.to_string) (b.n_states |> States_number.to_string) (b.n_states_diag |> States_number.to_string) - (b.n_det_max_jacobi |> Det_number.to_string) + (b.n_det_max_jacobi |> Strictly_positive_int.to_string) (b.threshold_generators |> Threshold.to_string) (b.threshold_selectors |> Threshold.to_string) (b.read_wf |> Bool.to_string) @@ -264,13 +451,96 @@ psi_det = %s (b.s2_eig |> Bool.to_string) (b.psi_coef |> Array.to_list |> List.map ~f:Det_coef.to_string |> String.concat ~sep:", ") - (b.psi_det |> Array.map ~f:(fun x -> Determinant.to_int64_array x - |> Array.map ~f:(fun x-> - Int64.to_string x )|> Array.to_list |> - String.concat ~sep:", ") |> Array.to_list - |> String.concat ~sep:" | ") - ; -;; + (b.psi_det |> Array.to_list |> List.map ~f:(Determinant.to_string + ~mo_tot_num:mo_tot_num) |> String.concat ~sep:"\n\n") + ;; + + let of_rst r = + let r = Rst_string.to_string r + in + + (* Split into header and determinants data *) + let idx = String.substr_index_exn r ~pos:0 ~pattern:"\nDeterminants" + in + let (header, dets) = + (String.prefix r idx, String.suffix r ((String.length r)-idx) ) + in + + (* Handle header *) + let header = r + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + if (line = "") then + false + else + ( (String.contains line '=') && (line.[0] = ' ') ) + ) + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + |> String.strip + )^")" ) + |> String.concat + in + + (* Handle determinant coefs *) + let dets = match ( dets + |> String.split ~on:'\n' + |> List.map ~f:(String.strip) + ) with + | _::lines -> lines + | _ -> failwith "Error in determinants" + in + + let psi_coef = + let rec read_coefs accu = function + | [] -> List.rev accu + | ""::c::tail -> + read_coefs (c::accu) tail + | _::tail -> read_coefs accu tail + in + let a = read_coefs [] dets + |> String.concat ~sep:" " + in + "(psi_coef ("^a^"))" + in + + (* Handle determinants *) + let psi_det = + let n_alpha = Ezfio.get_electrons_elec_alpha_num () + |> Elec_alpha_number.of_int + and n_beta = Ezfio.get_electrons_elec_beta_num () + |> Elec_beta_number.of_int + in + let rec read_dets accu = function + | [] -> List.rev accu + | ""::c::alpha::beta::tail -> + begin + let alpha = String.rev alpha |> Bitlist.of_string ~zero:'-' ~one:'+' + and beta = String.rev beta |> Bitlist.of_string ~zero:'-' ~one:'+' + in + let newdet = Determinant.of_bitlist_couple + ~alpha:n_alpha ~beta:n_beta (alpha,beta) + |> Determinant.sexp_of_t |> Sexplib.Sexp.to_string + in + read_dets (newdet::accu) tail + end + | _::tail -> read_dets accu tail + in + let a = read_dets [] dets + |> String.concat + in + "(psi_det ("^a^"))" + in + + let bitkind = Printf.sprintf "(bit_kind %d)" (Lazy.force Qpackage.bit_kind + |> Bit_kind.to_int) + and n_int = Printf.sprintf "(n_int %d)" (N_int_number.get_max ()) in + let s = String.concat [ header ; bitkind ; n_int ; psi_coef ; psi_det] + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; end diff --git a/ocaml/Input_electrons.ml b/ocaml/Input_electrons.ml index 75785aea..e4bc3ee2 100644 --- a/ocaml/Input_electrons.ml +++ b/ocaml/Input_electrons.ml @@ -6,17 +6,19 @@ module Electrons : sig type t = { elec_alpha_num : Elec_alpha_number.t; elec_beta_num : Elec_beta_number.t; - elec_num : Elec_number.t; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t -> unit + val read_elec_num : unit -> Elec_number.t val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = { elec_alpha_num : Elec_alpha_number.t; elec_beta_num : Elec_beta_number.t; - elec_num : Elec_number.t; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "electrons";; @@ -26,11 +28,22 @@ end = struct |> Elec_alpha_number.of_int ;; + let write_elec_alpha_num n = + Elec_alpha_number.to_int n + |> Ezfio.set_electrons_elec_alpha_num + ;; + + let read_elec_beta_num() = Ezfio.get_electrons_elec_beta_num () |> Elec_beta_number.of_int ;; + let write_elec_beta_num n = + Elec_beta_number.to_int n + |> Ezfio.set_electrons_elec_beta_num + ;; + let read_elec_num () = let na = Ezfio.get_electrons_elec_alpha_num () and nb = Ezfio.get_electrons_elec_beta_num () @@ -42,19 +55,57 @@ end = struct let read () = { elec_alpha_num = read_elec_alpha_num (); elec_beta_num = read_elec_beta_num (); - elec_num = read_elec_num (); } ;; - let to_string b = + let write { elec_alpha_num ; elec_beta_num } = + write_elec_alpha_num elec_alpha_num; + write_elec_beta_num elec_beta_num; + ;; + + + let to_rst b = Printf.sprintf " -elec_alpha_num = %s +Spin multiplicity is %s. + +Number of alpha and beta electrons :: + + elec_alpha_num = %s + elec_beta_num = %s + +" + (Multiplicity.of_alpha_beta b.elec_alpha_num b.elec_beta_num + |> Multiplicity.to_string) + (Elec_alpha_number.to_string b.elec_alpha_num) + (Elec_beta_number.to_string b.elec_beta_num) + |> Rst_string.of_string + ;; + + let to_string b = + Printf.sprintf "elec_alpha_num = %s elec_beta_num = %s elec_num = %s " (Elec_alpha_number.to_string b.elec_alpha_num) (Elec_beta_number.to_string b.elec_beta_num) - (Elec_number.to_string b.elec_num) + (Elec_number.to_string (read_elec_num ())) + ;; + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + end diff --git a/ocaml/Input_full_ci.ml b/ocaml/Input_full_ci.ml index bcff6681..08c305c8 100644 --- a/ocaml/Input_full_ci.ml +++ b/ocaml/Input_full_ci.ml @@ -4,19 +4,22 @@ open Core.Std;; module Full_ci : sig type t = - { n_det_max_fci : Det_number.t; + { n_det_max_fci : Det_number_max.t; pt2_max : PT2_energy.t; do_pt2_end : bool; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t-> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = - { n_det_max_fci : Det_number.t; + { n_det_max_fci : Det_number_max.t; pt2_max : PT2_energy.t; do_pt2_end : bool; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "full_ci";; @@ -28,7 +31,12 @@ end = struct |> Ezfio.set_full_ci_n_det_max_fci ; Ezfio.get_full_ci_n_det_max_fci () - |> Det_number.of_int + |> Det_number_max.of_int + ;; + + let write_n_det_max_fci ndet = + Det_number_max.to_int ndet + |> Ezfio.set_full_ci_n_det_max_fci ;; let read_pt2_max () = @@ -41,6 +49,11 @@ end = struct |> PT2_energy.of_float ;; + let write_pt2_max pt2_max = + PT2_energy.to_float pt2_max + |> Ezfio.set_full_ci_pt2_max + ;; + let read_do_pt2_end () = if not (Ezfio.has_full_ci_do_pt2_end ()) then get_default "do_pt2_end" @@ -50,6 +63,10 @@ end = struct Ezfio.get_full_ci_do_pt2_end () ;; + let write_do_pt2_end = + Ezfio.set_full_ci_do_pt2_end + ;; + let read () = { n_det_max_fci = read_n_det_max_fci (); @@ -58,15 +75,64 @@ end = struct } ;; + + let write { n_det_max_fci ; + pt2_max ; + do_pt2_end ; + } = + write_n_det_max_fci n_det_max_fci; + write_pt2_max pt2_max; + write_do_pt2_end do_pt2_end; + ;; + let to_string b = Printf.sprintf " -n_det_max_fci = %s -pt2_max = %s -do_pt2_end = %s +n_det_max_fci = %s +pt2_max = %s +do_pt2_end = %s " - (Det_number.to_string b.n_det_max_fci) + (Det_number_max.to_string b.n_det_max_fci) (PT2_energy.to_string b.pt2_max) (Bool.to_string b.do_pt2_end) + ;; + + let to_rst b = + Printf.sprintf " +Stop when the `n_det` > `n_det_max_fci` :: + + n_det_max_fci = %s + +Stop when -E(PT2) < `pt2_max` :: + + pt2_max = %s + +Compute E(PT2) at the end :: + + do_pt2_end = %s + +" + (Det_number_max.to_string b.n_det_max_fci) + (PT2_energy.to_string b.pt2_max) + (Bool.to_string b.do_pt2_end) + |> Rst_string.of_string + ;; + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + + end diff --git a/ocaml/Input_hartree_fock.ml b/ocaml/Input_hartree_fock.ml index a612ae74..56b18b31 100644 --- a/ocaml/Input_hartree_fock.ml +++ b/ocaml/Input_hartree_fock.ml @@ -6,15 +6,18 @@ module Hartree_fock : sig type t = { n_it_scf_max : Strictly_positive_int.t; thresh_scf : Threshold.t; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t -> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = { n_it_scf_max : Strictly_positive_int.t; thresh_scf : Threshold.t; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "hartree_fock";; @@ -29,14 +32,25 @@ end = struct |> Strictly_positive_int.of_int ;; - let read_thresh_scf() = + let write_n_it_scf_max n_it_scf_max = + Strictly_positive_int.to_int n_it_scf_max + |> Ezfio.set_hartree_fock_n_it_scf_max + ;; + + let read_thresh_scf () = if not (Ezfio.has_hartree_fock_thresh_scf()) then get_default "thresh_scf" |> Float.of_string |> Ezfio.set_hartree_fock_thresh_scf ; Ezfio.get_hartree_fock_thresh_scf () - |> Threshold.of_float ;; + |> Threshold.of_float + ;; + + let write_thresh_scf thresh_scf = + Threshold.to_float thresh_scf + |> Ezfio.set_hartree_fock_thresh_scf + ;; let read () = @@ -45,13 +59,56 @@ end = struct } ;; + + let write { n_it_scf_max ; + thresh_scf ; + } = + write_n_it_scf_max n_it_scf_max; + write_thresh_scf thresh_scf + ;; + + let to_string b = Printf.sprintf " -n_it_scf_max = %s -thresh_scf = %s +n_it_scf_max = %s +thresh_scf = %s +" + (Strictly_positive_int.to_string b.n_it_scf_max) + (Threshold.to_string b.thresh_scf) + ;; + + let to_rst b = + Printf.sprintf " +Max number of SCF iterations :: + + n_it_scf_max = %s + +SCF convergence criterion (on energy) :: + + thresh_scf = %s + " (Strictly_positive_int.to_string b.n_it_scf_max) (Threshold.to_string b.thresh_scf) + |> Rst_string.of_string + ;; + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + + end diff --git a/ocaml/Input_mo_basis.ml b/ocaml/Input_mo_basis.ml index fab7810a..d53e78b7 100644 --- a/ocaml/Input_mo_basis.ml +++ b/ocaml/Input_mo_basis.ml @@ -5,30 +5,31 @@ open Core.Std;; module Mo_basis : sig type t = { mo_tot_num : MO_number.t ; - mo_label : Non_empty_string.t; - mo_occ : Positive_float.t array; - mo_coef : MO_coef.t array; - } + mo_label : MO_label.t; + mo_occ : MO_occ.t array; + mo_coef : (MO_coef.t array) array; + } with sexp ;; val read : unit -> t val to_string : t -> string + val to_rst : t -> Rst_string.t end = struct type t = { mo_tot_num : MO_number.t ; - mo_label : Non_empty_string.t; - mo_occ : Positive_float.t array; - mo_coef : MO_coef.t array; - } + mo_label : MO_label.t; + mo_occ : MO_occ.t array; + mo_coef : (MO_coef.t array) array; + } with sexp ;; let get_default = Qpackage.get_ezfio_default "mo_basis";; let read_mo_label () = if not (Ezfio.has_mo_basis_mo_label ()) then - Ezfio.set_mo_basis_mo_label "Unknown" + Ezfio.set_mo_basis_mo_label "None" ; Ezfio.get_mo_basis_mo_label () - |> Non_empty_string.of_string + |> MO_label.of_string ;; let read_mo_tot_num () = @@ -50,15 +51,20 @@ end = struct ~dim:[| mo_tot_num |] ~data:data |> Ezfio.set_mo_basis_mo_occ end; - (Ezfio.get_mo_basis_mo_occ () ).Ezfio.data - |> Ezfio.flattened_ezfio_data - |> Array.map ~f:Positive_float.of_float + Ezfio.flattened_ezfio (Ezfio.get_mo_basis_mo_occ () ) + |> Array.map ~f:MO_occ.of_float ;; let read_mo_coef () = - (Ezfio.get_mo_basis_mo_coef () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + let a = Ezfio.get_mo_basis_mo_coef () + |> Ezfio.flattened_ezfio |> Array.map ~f:MO_coef.of_float + in + let mo_tot_num = read_mo_tot_num () |> MO_number.to_int in + let ao_num = (Array.length a)/mo_tot_num in + Array.init mo_tot_num ~f:(fun j -> + Array.sub ~pos:(j*ao_num) ~len:(ao_num) a + ) ;; let read () = @@ -69,6 +75,91 @@ end = struct } ;; + let mo_coef_to_string mo_coef = + let ao_num = Array.length mo_coef.(0) + and mo_tot_num = Array.length mo_coef in + let rec print_five imin imax = + match (imax-imin+1) with + | 1 -> + let header = [ Printf.sprintf " #%15d" (imin+1) ; ] in + let new_lines = + List.init ao_num ~f:(fun i -> + Printf.sprintf " %3d %15.10f " (i+1) + (MO_coef.to_float mo_coef.(imin ).(i)) ) + in header @ new_lines + | 2 -> + let header = [ Printf.sprintf " #%15d %15d" (imin+1) (imin+2) ; ] in + let new_lines = + List.init ao_num ~f:(fun i -> + Printf.sprintf " %3d %15.10f %15.10f" (i+1) + (MO_coef.to_float mo_coef.(imin ).(i)) + (MO_coef.to_float mo_coef.(imin+1).(i)) ) + in header @ new_lines + | 3 -> + let header = [ Printf.sprintf " #%15d %15d %15d" + (imin+1) (imin+2) (imin+3); ] in + let new_lines = + List.init ao_num ~f:(fun i -> + Printf.sprintf " %3d %15.10f %15.10f %15.10f" (i+1) + (MO_coef.to_float mo_coef.(imin ).(i)) + (MO_coef.to_float mo_coef.(imin+1).(i)) + (MO_coef.to_float mo_coef.(imin+2).(i)) ) + in header @ new_lines + | 4 -> + let header = [ Printf.sprintf " #%15d %15d %15d %15d" + (imin+1) (imin+2) (imin+3) (imin+4) ; ] in + let new_lines = + List.init ao_num ~f:(fun i -> + Printf.sprintf " %3d %15.10f %15.10f %15.10f %15.10f" (i+1) + (MO_coef.to_float mo_coef.(imin ).(i)) + (MO_coef.to_float mo_coef.(imin+1).(i)) + (MO_coef.to_float mo_coef.(imin+2).(i)) + (MO_coef.to_float mo_coef.(imin+3).(i)) ) + in header @ new_lines + | 5 -> + let header = [ Printf.sprintf " #%15d %15d %15d %15d %15d" + (imin+1) (imin+2) (imin+3) (imin+4) (imin+5) ; ] in + let new_lines = + List.init ao_num ~f:(fun i -> + Printf.sprintf " %3d %15.10f %15.10f %15.10f %15.10f %15.10f" (i+1) + (MO_coef.to_float mo_coef.(imin ).(i)) + (MO_coef.to_float mo_coef.(imin+1).(i)) + (MO_coef.to_float mo_coef.(imin+2).(i)) + (MO_coef.to_float mo_coef.(imin+3).(i)) + (MO_coef.to_float mo_coef.(imin+4).(i)) ) + in header @ new_lines + | _ -> assert false + in + let rec create_list accu i = + if (i+4 < mo_tot_num) then + create_list ( (print_five i (i+3) |> String.concat ~sep:"\n")::accu ) (i+4) + else + (print_five i (mo_tot_num-1) |> String.concat ~sep:"\n")::accu |> List.rev + in + create_list [] 0 |> String.concat ~sep:"\n\n" + ;; + + let to_rst b = + Printf.sprintf " +Label of the molecular orbitals :: + + mo_label = %s + +Total number of MOs :: + + mo_tot_num = %s + +MO coefficients :: + +%s +" + (MO_label.to_string b.mo_label) + (MO_number.to_string b.mo_tot_num) + (mo_coef_to_string b.mo_coef) + |> Rst_string.of_string + + ;; + let to_string b = Printf.sprintf " mo_label = %s @@ -76,12 +167,15 @@ mo_tot_num = \"%s\" mo_occ = %s mo_coef = %s " - (Non_empty_string.to_string b.mo_label) + (MO_label.to_string b.mo_label) (MO_number.to_string b.mo_tot_num) (b.mo_occ |> Array.to_list |> List.map - ~f:(Positive_float.to_string) |> String.concat ~sep:", " ) - (b.mo_coef |> Array.to_list |> List.map - ~f:(MO_coef.to_string) |> String.concat ~sep:", " ) + ~f:(MO_occ.to_string) |> String.concat ~sep:", " ) + (b.mo_coef |> Array.map + ~f:(fun x-> Array.map ~f:MO_coef.to_string x |> String.concat_array + ~sep:"," ) |> + String.concat_array ~sep:"\n" ) + ;; end diff --git a/ocaml/Input_nuclei.ml b/ocaml/Input_nuclei.ml index c4139edb..11eeb7d4 100644 --- a/ocaml/Input_nuclei.ml +++ b/ocaml/Input_nuclei.ml @@ -8,43 +8,80 @@ module Nuclei : sig nucl_label : Element.t array; nucl_charge : Charge.t array; nucl_coord : Point3d.t array; - } + } with sexp ;; - val read : unit -> t + val read : unit -> t + val write : t -> unit val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t end = struct type t = { nucl_num : Nucl_number.t ; nucl_label : Element.t array; nucl_charge : Charge.t array; nucl_coord : Point3d.t array; - } + } with sexp ;; let get_default = Qpackage.get_ezfio_default "nuclei";; let read_nucl_num () = - Ezfio.get_nuclei_nucl_num () - |> Nucl_number.of_int + let nmax = Nucl_number.get_max () in + Nucl_number.of_int ~max:nmax nmax ;; + let write_nucl_num n = + Nucl_number.to_int n + |> Ezfio.set_nuclei_nucl_num + ;; + + let read_nucl_label () = - (Ezfio.get_nuclei_nucl_label ()).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_nuclei_nucl_label () + |> Ezfio.flattened_ezfio |> Array.map ~f:Element.of_string ;; + let write_nucl_label ~nucl_num labels = + let nucl_num = + Nucl_number.to_int nucl_num + in + let labels = + Array.to_list labels + |> List.map ~f:Element.to_string + in + Ezfio.ezfio_array_of_list ~rank:1 + ~dim:[| nucl_num |] ~data:labels + |> Ezfio.set_nuclei_nucl_label + ;; + + let read_nucl_charge () = - (Ezfio.get_nuclei_nucl_charge () ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_nuclei_nucl_charge () + |> Ezfio.flattened_ezfio |> Array.map ~f:Charge.of_float ;; + let write_nucl_charge ~nucl_num charges = + let nucl_num = + Nucl_number.to_int nucl_num + in + let charges = + Array.to_list charges + |> List.map ~f:Charge.to_float + in + Ezfio.ezfio_array_of_list ~rank:1 + ~dim:[| nucl_num |] ~data:charges + |> Ezfio.set_nuclei_nucl_charge + ;; + + let read_nucl_coord () = let nucl_num = Nucl_number.to_int (read_nucl_num ()) in let raw_data = - (Ezfio.get_nuclei_nucl_coord() ).Ezfio.data - |> Ezfio.flattened_ezfio_data + Ezfio.get_nuclei_nucl_coord() + |> Ezfio.flattened_ezfio in let zero = Point3d.of_string Units.Bohr "0. 0. 0." in let result = Array.create nucl_num zero in @@ -57,6 +94,22 @@ end = struct result ;; + let write_nucl_coord ~nucl_num coord = + let nucl_num = + Nucl_number.to_int nucl_num + in + let coord = Array.to_list coord in + let coord = + (List.map ~f:(fun x-> x.Point3d.x) coord) @ + (List.map ~f:(fun x-> x.Point3d.y) coord) @ + (List.map ~f:(fun x-> x.Point3d.z) coord) + in + Ezfio.ezfio_array_of_list ~rank:2 + ~dim:[| nucl_num ; 3 |] ~data:coord + |> Ezfio.set_nuclei_nucl_coord + ;; + + let read () = { nucl_num = read_nucl_num (); nucl_label = read_nucl_label () ; @@ -65,6 +118,18 @@ end = struct } ;; + let write { nucl_num ; + nucl_label ; + nucl_charge ; + nucl_coord ; + } = + write_nucl_num nucl_num ; + write_nucl_label ~nucl_num:nucl_num nucl_label; + write_nucl_charge ~nucl_num:nucl_num nucl_charge; + write_nucl_coord ~nucl_num:nucl_num nucl_coord; + ;; + + let to_string b = Printf.sprintf " nucl_num = %s @@ -79,6 +144,75 @@ nucl_coord = %s ~f:(Charge.to_string) |> String.concat ~sep:", " ) (b.nucl_coord |> Array.to_list |> List.map ~f:(Point3d.to_string Units.Bohr) |> String.concat ~sep:"\n" ) + ;; + + + let to_rst b = + let nucl_num = Nucl_number.to_int b.nucl_num in + let text = + ( Printf.sprintf " %d\n " + nucl_num + ) :: ( + List.init nucl_num ~f:(fun i-> + Printf.sprintf " %-3s %d %s" + (b.nucl_label.(i) |> Element.to_string) + (b.nucl_charge.(i) |> Charge.to_int ) + (b.nucl_coord.(i) |> Point3d.to_string Units.Angstrom) ) + ) |> String.concat ~sep:"\n" + in + Printf.sprintf " +Nuclear coordinates in xyz format (Angstroms) :: + +%s + +" text + |> Rst_string.of_string + ;; + + let of_rst s = + let l = Rst_string.to_string s + |> String.split ~on:'\n' + in + (* Find lines containing the xyz data *) + let rec extract_begin = function + | [] -> raise Not_found + | line::tail -> + let line = String.strip line in + if (String.length line > 3) && + (String.sub line ~pos:((String.length line)-2) + ~len:2 = "::") then + tail + else + extract_begin tail + in + (* Create a list of Atom.t *) + let nmax = Nucl_number.get_max () in + let atom_list = + match (extract_begin l) with + | _ :: nucl_num :: title :: lines -> + begin + let nucl_num = nucl_num + |> String.strip + |> Int.of_string + |> Nucl_number.of_int ~max:nmax + and lines = Array.of_list lines + in + List.init (Nucl_number.to_int nucl_num) ~f:(fun i -> + Atom.of_string Units.Angstrom lines.(i)) + end + | _ -> failwith "Error in xyz format" + in + (* Create the Nuclei.t data structure *) + { nucl_num = List.length atom_list + |> Nucl_number.of_int ~max:nmax; + nucl_label = List.map atom_list ~f:(fun x -> + x.Atom.element) |> Array.of_list ; + nucl_charge = List.map atom_list ~f:(fun x -> + x.Atom.charge ) |> Array.of_list ; + nucl_coord = List.map atom_list ~f:(fun x -> + x.Atom.coord ) |> Array.of_list ; + } + ;; end diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index c59a628a..38091079 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -1,7 +1,7 @@ open Core.Std;; open Qptypes;; -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list;; +type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list with sexp let of_basis b = let rec do_work accu = function @@ -19,11 +19,36 @@ let of_basis b = |> List.rev ;; +let to_basis b = + let rec do_work accu = function + | [] -> List.rev accu + | (s,g,n)::tail -> + let first_sym = + Symmetry.Xyz.of_symmetry g.Gto.sym + |> List.hd_exn + in + let new_accu = + if ( s = first_sym ) then + (g,n)::accu + else + accu + in + do_work new_accu tail + in + do_work [] b +;; let to_string b = - List.map ~f:(fun (x,y,z) -> - (Int.to_string (Nucl_number.to_int z))^":"^ - (Symmetry.Xyz.to_string x)^" "^(Gto.to_string y) + let middle = List.map ~f:(fun (x,y,z) -> + "( "^((Int.to_string (Nucl_number.to_int z)))^", "^ + (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) + ^" )" ) b - |> String.concat ~sep:"\n" + |> String.concat ~sep:",\n" + in "("^middle^")" ;; + +include To_md5;; +let to_md5 = to_md5 sexp_of_t +;; + diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli index 6c9c8db3..7e69ecce 100644 --- a/ocaml/Long_basis.mli +++ b/ocaml/Long_basis.mli @@ -5,12 +5,16 @@ open Qptypes;; * all the D orbitals are converted to xx, xy, xz, yy, yx * etc *) -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list +type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list with sexp (** Transform a basis to a long basis *) val of_basis : (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list +(** Transform a long basis to a basis *) +val to_basis : + (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list + (** Convert the basis into its string representation *) val to_string : (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string diff --git a/ocaml/MO_class.ml b/ocaml/MO_class.ml index 93b8007b..4fc03da2 100644 --- a/ocaml/MO_class.ml +++ b/ocaml/MO_class.ml @@ -7,7 +7,7 @@ type t = | Active of MO_number.t list | Virtual of MO_number.t list | Deleted of MO_number.t list -;; +with sexp let to_string x = diff --git a/ocaml/MO_class.mli b/ocaml/MO_class.mli new file mode 100644 index 00000000..057d4b20 --- /dev/null +++ b/ocaml/MO_class.mli @@ -0,0 +1,21 @@ +type t = + | Core of Qptypes.MO_number.t list + | Inactive of Qptypes.MO_number.t list + | Active of Qptypes.MO_number.t list + | Virtual of Qptypes.MO_number.t list + | Deleted of Qptypes.MO_number.t list +with sexp + + +(** Create different excitation classes *) +val create_core : string -> t +val create_inactive : string -> t +val create_active : string -> t +val create_virtual : string -> t +val create_deleted : string -> t + +(** Convert to a Bitlist.t *) +val to_bitlist : Qptypes.N_int_number.t -> t -> Bitlist.t + +(** Convert to string for printing *) +val to_string : t -> string diff --git a/ocaml/MO_label.ml b/ocaml/MO_label.ml new file mode 100644 index 00000000..f8de6210 --- /dev/null +++ b/ocaml/MO_label.ml @@ -0,0 +1,32 @@ +open Core.Std;; + +type t = +| Guess +| Canonical +| Natural +| Localized +| Orthonormalized +| None +with sexp +;; + +let to_string = function + | Guess -> "Guess" + | Canonical -> "Canonical" + | Orthonormalized -> "Orthonormalized" + | Natural -> "Natural" + | Localized -> "Localized" + | None -> "None" +;; + +let of_string s = + match String.lowercase s with + | "guess" -> Guess + | "canonical" -> Canonical + | "natural" -> Natural + | "localized" -> Localized + | "orthonormalized" -> Orthonormalized + | "none" -> None + | _ -> failwith "MO_label should be one of: + Guess | Orthonormalized | Canonical | Natural | Localized | None." +;; diff --git a/ocaml/MO_label.mli b/ocaml/MO_label.mli new file mode 100644 index 00000000..d5061095 --- /dev/null +++ b/ocaml/MO_label.mli @@ -0,0 +1,15 @@ +type t = + | Guess + | Canonical + | Natural + | Localized + | Orthonormalized + | None +with sexp + +(** String representation *) +val to_string : t -> string + +(** Build from string representation *) +val of_string : string -> t + diff --git a/ocaml/Makefile b/ocaml/Makefile index 82838cee..61e531a3 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -1,5 +1,3 @@ -#TODO : Opam auto-installer in makefile - # Check if QPACKAGE_ROOT is defined ifndef QPACKAGE_ROOT @@ -12,8 +10,8 @@ endif LIBS= PKGS= -OCAMLCFLAGS=-g -OCAMLBUILD=ocamlbuild -j 0 -cflags $(OCAMLCFLAGS) -lflags -g +OCAMLCFLAGS="-g" +OCAMLBUILD=ocamlbuild -j 0 -syntax camlp4o -cflags $(OCAMLCFLAGS) -lflags $(OCAMLCFLAGS) MLFILES=$(wildcard *.ml) ezfio.ml Qptypes.ml MLIFILES=$(wildcard *.mli) ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml)) @@ -21,11 +19,21 @@ ALL_EXE=$(patsubst %.ml,%.native,$(wildcard qp_*.ml)) .PHONY: executables default + default: $(ALL_TESTS) $(ALL_EXE) executables: $(MAKE) -C $(QPACKAGE_ROOT)/data executables +external_libs: + opam install cryptokit core + +qpackage.odocl: $(MLIFILES) + ls $(MLIFILES) | sed "s/\.mli//" > qpackage.odocl + +doc: qpackage.odocl + $(OCAMLBUILD) qpackage.docdir/index.html -use-ocamlfind $(PKGS) + %.inferred.mli: $(MLFILES) $(OCAMLBUILD) $*.inferred.mli -use-ocamlfind $(PKGS) mv _build/$*.inferred.mli . diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 7592b85a..f295420b 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -7,7 +7,7 @@ type t = { nuclei : Atom.t list ; elec_alpha : Elec_alpha_number.t ; elec_beta : Elec_beta_number.t ; -} +} with sexp let get_charge { nuclei ; elec_alpha ; elec_beta } = let result = (Elec_alpha_number.to_int elec_alpha) + @@ -25,7 +25,8 @@ let get_multiplicity m = ;; let get_nucl_num m = - Nucl_number.of_int (List.length m.nuclei) + let nmax = (List.length m.nuclei) in + Nucl_number.of_int nmax ~max:nmax ;; let name m = @@ -78,7 +79,7 @@ let to_string m = ;; let of_xyz_string - ?(charge=0) ?(multiplicity=(Multiplicity.of_int 1)) + ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) ?(units=Units.Angstrom) s = let l = String.split s ~on:'\n' @@ -90,7 +91,7 @@ let of_xyz_string elec_alpha=(Elec_alpha_number.of_int 1) ; elec_beta=(Elec_beta_number.of_int 0) } |> Charge.to_int - ) + 1 - charge + ) + 1 - (Charge.to_int charge) |> Elec_number.of_int in let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in @@ -112,10 +113,16 @@ let of_xyz_string let of_xyz_file - ?(charge=0) ?(multiplicity=(Multiplicity.of_int 1)) + ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) + ?(units=Units.Angstrom) filename = let (_,buffer) = In_channel.read_all filename |> String.lsplit2_exn ~on:'\n' in let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in - of_xyz_string ~charge:charge ~multiplicity:multiplicity buffer + of_xyz_string ~charge:charge ~multiplicity:multiplicity + ~units:units buffer +;; + +include To_md5;; +let to_md5 = to_md5 sexp_of_t ;; diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli new file mode 100644 index 00000000..176441d4 --- /dev/null +++ b/ocaml/Molecule.mli @@ -0,0 +1,38 @@ +exception MultiplicityError of string + +type t = { + nuclei : Atom.t list; + elec_alpha : Qptypes.Elec_alpha_number.t; + elec_beta : Qptypes.Elec_beta_number.t; +} with sexp + +(** Returns the charge of the molecule *) +val get_charge : t -> Charge.t + +(** Returns the multiplicity of the molecule *) +val get_multiplicity : t -> Multiplicity.t + +(** Returns the number of nuclei *) +val get_nucl_num : t -> Qptypes.Nucl_number.t + +(** The name of the molecule *) +val name : t -> string + +(** Conversion for printing *) +val to_string : t -> string + + +(** Creates a molecule from an xyz file *) +val of_xyz_file : + ?charge:Charge.t -> + ?multiplicity:Multiplicity.t -> + ?units:Units.units -> string -> t + +(** Creates a molecule from an xyz file in a string *) +val of_xyz_string : + ?charge:Charge.t -> + ?multiplicity:Multiplicity.t -> + ?units:Units.units -> string -> t + +(** Computes the MD5 hash *) +val to_md5 : t -> Qptypes.MD5.t diff --git a/ocaml/Multiplicity.ml b/ocaml/Multiplicity.ml index a66233a8..e56341c6 100644 --- a/ocaml/Multiplicity.ml +++ b/ocaml/Multiplicity.ml @@ -1,7 +1,8 @@ open Core.Std;; open Qptypes ;; -type t = Strictly_positive_int.t;; +type t = Strictly_positive_int.t with sexp + let of_int = Strictly_positive_int.of_int ;; let to_int = Strictly_positive_int.to_int ;; diff --git a/ocaml/Multiplicity.mli b/ocaml/Multiplicity.mli new file mode 100644 index 00000000..c6f8c6bf --- /dev/null +++ b/ocaml/Multiplicity.mli @@ -0,0 +1,19 @@ +type t = Qptypes.Strictly_positive_int.t with sexp + +(** Conversion from int *) +val of_int : int -> t +val to_int : t -> int + +(** Computation from the number of alpha and beta electrons *) +val of_alpha_beta : + Qptypes.Elec_alpha_number.t -> + Qptypes.Elec_beta_number.t -> t + +(** Generation of the number of alpha and beta electrons *) +val to_alpha_beta : + Qptypes.Elec_number.t -> t -> + Qptypes.Elec_alpha_number.t * Qptypes.Elec_beta_number.t + +(** Conversion to string for printing *) +val to_string : t-> string + diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index daafd139..18f091f1 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -1,10 +1,11 @@ open Core.Std;; +open Qptypes;; type t = { x : float ; y : float ; z : float ; -} +} with sexp (** Read x y z coordinates in string s with units u *) let of_string u s = @@ -28,9 +29,10 @@ let distance2 p1 p2 = let { x=x1 ; y=y1 ; z=z1 } = p1 and { x=x2 ; y=y2 ; z=z2 } = p2 in (x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1) + |> Positive_float.of_float ;; -let distance p1 p2 = sqrt (distance2 p1 p2) +let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2)) ;; let to_string u p = @@ -39,6 +41,6 @@ let to_string u p = | Units.Angstrom -> Units.bohr_to_angstrom in let { x=x ; y=y ; z=z } = p in - Printf.sprintf "%f %f %f" (x*.f) (y*.f) (z*.f) + Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f) ;; diff --git a/ocaml/Point3d.mli b/ocaml/Point3d.mli new file mode 100644 index 00000000..066a4514 --- /dev/null +++ b/ocaml/Point3d.mli @@ -0,0 +1,17 @@ +type t = +{ x : float; + y : float; + z : float; +} with sexp + +(** Create from an xyz string *) +val of_string : Units.units -> string -> t + +(** Convert to a string for printing *) +val to_string : Units.units -> t -> string + +(** Computes the squared distance between 2 points *) +val distance2 : t -> t -> Qptypes.Positive_float.t + +(** Computes the distance between 2 points *) +val distance : t -> t -> float diff --git a/ocaml/Primitive.ml b/ocaml/Primitive.ml index 82a94fa4..a6377d6f 100644 --- a/ocaml/Primitive.ml +++ b/ocaml/Primitive.ml @@ -4,7 +4,7 @@ open Core.Std;; type t = { sym : Symmetry.t ; expo : AO_expo.t ; -} +} with sexp let to_string p = let { sym = s ; expo = e } = p in @@ -13,3 +13,6 @@ let to_string p = (AO_expo.to_float e) ;; +let of_sym_expo s e = + { sym=s ; expo=e} +;; diff --git a/ocaml/Primitive.mli b/ocaml/Primitive.mli new file mode 100644 index 00000000..77cb633a --- /dev/null +++ b/ocaml/Primitive.mli @@ -0,0 +1,11 @@ +type t = +{ sym : Symmetry.t; + expo : Qptypes.AO_expo.t; +} with sexp + +(** Conversion to string for printing *) +val to_string : t -> string + +(** Creation *) +val of_sym_expo : Symmetry.t -> Qptypes.AO_expo.t -> t + diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index fd40547f..ed112de3 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -1,3 +1,6 @@ +open Core.Std;; + +(* let rec transpose = function | [] -> [] | []::tail -> transpose tail @@ -7,5 +10,20 @@ let rec transpose = function in new_head @ new_tail ;; +*) +let input_to_sexp s = + let result = + String.split_lines s + |> List.filter ~f:(fun x-> + (String.strip x) <> "") + |> List.map ~f:(fun x-> + "("^(String.tr '=' ' ' x)^")") + |> String.concat + in + print_endline ("("^result^")"); + "("^result^")" + |> Sexp.of_string +;; + diff --git a/ocaml/README.rst b/ocaml/README.rst index 414428dd..2f33b268 100644 --- a/ocaml/README.rst +++ b/ocaml/README.rst @@ -4,3 +4,9 @@ Ocaml scripts This directory contains all the scripts that control the input/output with the user. + +All executables start with `qp_` and all tests start with `test_`. Modules +file names start with a capital letter. + +Info on how to extend the `qp_edit` tool is given in `README_qp_edit.rst`. + diff --git a/ocaml/README_qp_edit.rst b/ocaml/README_qp_edit.rst new file mode 100644 index 00000000..4525d21a --- /dev/null +++ b/ocaml/README_qp_edit.rst @@ -0,0 +1,220 @@ +Adding a new block +================== + +In this section, we assume we will add the `New_keyword` keyword. + +Create the `Input_new_keyword.ml` file +-------------------------------------- + +Copy for example the `Input_full_ci.ml` file as a starting point. + +The template is the following, where `r_x`, `r_y`, ..., `last_r` are the records +of the block. + +.. code-block:: ocaml + + module New_keyword : sig + type t = + { r_x : Type_of_x.t + r_y : Y_type.t + ... + last_r : bool + } with sexp + ;; + val read : unit -> t + val write : t -> unit + val to_string : t -> string + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t + end = struct + type t = + { r_x : Type_of_x.t + r_y : Y_type.t + ... + last_r : bool + } with sexp + ;; + + let get_default = Qpackage.get_ezfio_default "new_keyword";; + + + ... + + end + +The following functions need to be defined:: + + val read : unit -> t + val write : t -> unit + val to_rst : t -> Rst_string.t + val of_rst : Rst_string.t -> t + + +The type `t` has to be defined in a same way in the `sig` and the `struct`. + +For each record of the type `t`, use types defined in the `Qptypes.ml` file as +much as possible. + +The `get_default` function will fetch the default values in the `ezfio_defaults` file +in the `new_keyword` block. + +For each record `r_x` of the type `t`, create a `read_r_x ()` function +and a `write_r_x r_x` function that performs the I/O in the EZFIO. +To set a default value in the `read_r_x` function, use the following template +(assuming that the `Type_of_x` is built from a `double precision` value in +the EZFIO file). + +.. code-block:: ocaml + + let read_r_x () = + if not (Ezfio.has_new_keyword_r_x ()) then + get_default "r_x" + |> Float.of_string + |> Ezfio.set_new_keyword_r_x + ; + Ezfio.get_new_keyword_r_x () + |> Type_of_x.of_float + ;; + + let write_r_x r_x = + Type_of_x.to_float r_x + |> Ezfio.set_new_keyword_r_x + ;; + + +Then, create a `read` and a `write` function as + +.. code-block:: ocaml + + let read () = + { r_x = read_r_x () ; + r_y = read_r_y () ; + ... + last_r = read_last_r () ; + } + ;; + + let write { r_x ; + r_y + ... + last_r ; + } = + write_r_x r_x; + write_r_y r_y; + ... + write_last_r last_r; + ;; + +Finally, create the functions to write an RST string as + +.. code-block:: ocaml + + let to_rst b = + Printf.sprintf " + You can put here some documentation as long as there is no equal sign. + The record entries should be indented on the right with a blank line + before and a blank line after, as they would be in a rst file. + + Here is the text for r_x + + r_x = %s + + And here is the text for r_y + + r_y = %s + + ... + Finally, the text for last_r + + last_r = %s + " + (Type_of_x.to_string b.r_x) + (Y_type.to_string b.r_y) + ... + (Bool.to_string b.last_r) + ;; + + +and you can use this function to read it back: + +.. code-block:: ocaml + + let of_rst s = + let s = Rst_string.to_string s + |> String.split ~on:'\n' + |> List.filter ~f:(fun line -> + String.contains line '=') + |> List.map ~f:(fun line -> + "("^( + String.tr line ~target:'=' ~replacement:' ' + )^")" ) + |> String.concat + in + Sexp.of_string ("("^s^")") + |> t_of_sexp + ;; + + + + + +Add module to `Input.ml` file +----------------------------- + +Append module to the `Input.ml` file. Use the name of the `Input_new_keyword.ml` without the +`.ml` suffix. + +.. code-block:: ocaml + + include Input_new_keyword;; + + +In the `qp_edit.ml` file +------------------------ + +vim search strings are given in brackets. + +1. (`/type keyword`) : Add a new entry to the keyword type corresponding to the block to add: + +.. code-block:: ocaml + + type keyword = + ... + | New_keyword + ;; + + + +2. (`/keyword_to_string`) : Add a new entry to the `keyword_to_string` function for the title of the block + +.. code-block:: ocaml + + let keyword_to_string = function + ... + | New_keyword -> "My new keyword" + ;; + + +3. (`/let get s`) : Add a new call to the to_rst function of the `Input.New_keyword` module + +.. code-block:: ocaml + + let get s = + let header = (make_header s) + and rst = match s with + ... + | New_keyword -> + Input.New_keyword.(to_rst (read ())) + ... + + +4. (`/let set s`) : Add a new call to the of_rst function of the `Input.New_keyword` module + +.. code-block:: ocaml + + match s with + ... + | New_keyword -> + Input.New_keyword.(write (of_rst str)) + ;; + diff --git a/ocaml/Range.ml b/ocaml/Range.ml index d93be879..9607c6a3 100644 --- a/ocaml/Range.ml +++ b/ocaml/Range.ml @@ -12,7 +12,7 @@ open Core.Std;; *) -type t = int list ;; +type t = int list with sexp let expand_range r = match String.lsplit2 ~on:'-' r with diff --git a/ocaml/Range.mli b/ocaml/Range.mli new file mode 100644 index 00000000..2d56a0fa --- /dev/null +++ b/ocaml/Range.mli @@ -0,0 +1,10 @@ +type t = int list with sexp + +(** A range is a sorted list of ints in an interval. + It is created using a string : + "[a-b]" : range between a and b (included) + "[a]" : the list with only one integer a + "a" : equivalent to "[a]" + *) +val of_string : string -> t +val to_string : t -> string diff --git a/ocaml/Symmetry.ml b/ocaml/Symmetry.ml index 828d0dfa..9e130606 100644 --- a/ocaml/Symmetry.ml +++ b/ocaml/Symmetry.ml @@ -1,7 +1,7 @@ open Qptypes;; open Core.Std;; -type t = S|P|D|F|G|H|I|J|K|L +type t = S|P|D|F|G|H|I|J|K|L with sexp let to_string = function | S -> "S" @@ -53,22 +53,31 @@ let to_l = function | J -> Positive_int.of_int 7 | K -> Positive_int.of_int 8 | L -> Positive_int.of_int 9 +;; + +let of_l i = + let i = Positive_int.to_int i in + match i with + | 0 -> S + | 1 -> P + | 2 -> D + | 3 -> F + | 4 -> G + | 5 -> H + | 6 -> I + | 7 -> J + | 8 -> K + | 9 -> L + | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L")) +;; type st = t ;; -module Xyz : sig +module Xyz = struct type t = { x: Positive_int.t ; y: Positive_int.t ; - z: Positive_int.t } - val of_string : string -> t - val to_string : t -> string - val get_l : t -> Positive_int.t - val of_symmetry : st -> t list -end = struct - type t = { x: Positive_int.t ; - y: Positive_int.t ; - z: Positive_int.t } + z: Positive_int.t } with sexp type state_type = Null | X | Y | Z (** Builds an XYZ triplet from a string. @@ -127,7 +136,9 @@ end = struct | 1 -> "z" | i -> Printf.sprintf "z%d" i in - x^y^z + let result = (x^y^z) in + if (result = "") then "s" + else result ;; (** Returns the l quantum number from a XYZ powers triplet *) @@ -167,7 +178,10 @@ end = struct in create_x [] { x=(to_l sym) ; y=Positive_int.of_int 0 ; z=Positive_int.of_int 0 } - ;; + ;; + (** Returns the symmetry corresponding to the XYZ triplet *) + let to_symmetry sym = of_l (get_l sym) + ;; end diff --git a/ocaml/Symmetry.mli b/ocaml/Symmetry.mli new file mode 100644 index 00000000..baffeb2e --- /dev/null +++ b/ocaml/Symmetry.mli @@ -0,0 +1,36 @@ +type t = S | P | D | F | G | H | I | J | K | L with sexp + +(** Creatio from strings *) +val to_string : t -> string +val of_string : string -> t +val of_char : char -> t + +(** Connexion with l quantum number *) +val to_l : t -> Qptypes.Positive_int.t +val of_l : Qptypes.Positive_int.t -> t + +type st = t +module Xyz : + sig + type t = { + x : Qptypes.Positive_int.t; + y : Qptypes.Positive_int.t; + z : Qptypes.Positive_int.t; + } with sexp + + (** The string format contains the powers of x,y and z in a + format like "x2z3" *) + + val of_string : string -> t + val to_string : t -> string + + (** Returns the quantum number l *) + val get_l : t -> Qptypes.Positive_int.t + + (** Returns a list of XYZ powers for a given symmetry *) + val of_symmetry : st -> t list + + (** Returns the symmetry corresponding to the XYZ powers *) + val to_symmetry : t -> st + + end diff --git a/ocaml/To_md5.ml b/ocaml/To_md5.ml new file mode 100644 index 00000000..33015d75 --- /dev/null +++ b/ocaml/To_md5.ml @@ -0,0 +1,11 @@ +open Core.Std;; +open Qptypes;; + +let to_md5 sexp_of_t t = + sexp_of_t t + |> Sexp.to_string + |> Cryptokit.hash_string (Cryptokit.Hash.md5 ()) + |> Cryptokit.transform_string (Cryptokit.Hexa.encode ()) + |> MD5.of_string +;; + diff --git a/ocaml/Units.mli b/ocaml/Units.mli new file mode 100644 index 00000000..eff3fb23 --- /dev/null +++ b/ocaml/Units.mli @@ -0,0 +1,7 @@ +type units = Bohr | Angstrom + +(** Conversion functions *) +val angstrom_to_bohr : float +val bohr_to_angstrom : float + + diff --git a/ocaml/_tags b/ocaml/_tags index 55cb5085..519d558f 100644 --- a/ocaml/_tags +++ b/ocaml/_tags @@ -1,3 +1,2 @@ -true: package(core) +true: package(core,sexplib.syntax,cryptokit) true: thread - diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index d35cb42b..441519ea 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -26,8 +26,8 @@ let run ?o b c m xyz_file = (* Read molecule *) let molecule = - Molecule.of_xyz_file xyz_file ~charge:c - ~multiplicity:(Multiplicity.of_int m) + (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c) + ~multiplicity:(Multiplicity.of_int m) ) in (* Build EZFIO File name *) @@ -74,9 +74,10 @@ let run ?o b c m xyz_file = (* Write Basis set *) let basis = + let nmax = Nucl_number.get_max () in let rec do_work (accu:(Atom.t*Nucl_number.t) list) (n:int) = function | [] -> accu - | e::tail -> let new_accu = (e,(Nucl_number.of_int n))::accu in + | e::tail -> let new_accu = (e,(Nucl_number.of_int ~max:nmax n))::accu in do_work new_accu (n+1) tail in do_work [] 1 nuclei diff --git a/ocaml/qp_edit.ml b/ocaml/qp_edit.ml new file mode 100644 index 00000000..e8c93b9f --- /dev/null +++ b/ocaml/qp_edit.ml @@ -0,0 +1,209 @@ +open Qputils;; +open Qptypes;; +open Core.Std;; + +let file_header filename = Printf.sprintf +" +================================================================== + Quantum Package +================================================================== + +Editing file `%s` + +" filename + +type keyword = +| Ao_basis +| Bielec_integrals +| Cisd_sc2 +| Determinants +| Electrons +| Full_ci +| Hartree_fock +| Mo_basis +| Nuclei +;; + +let keyword_to_string = function +| Ao_basis -> "AO basis" +| Bielec_integrals -> "Two electron integrals" +| Cisd_sc2 -> "CISD (SC)^2" +| Determinants -> "Determinants" +| Electrons -> "Electrons" +| Full_ci -> "Selected Full-CI" +| Hartree_fock -> "Hartree-Fock" +| Mo_basis -> "MO basis" +| Nuclei -> "Molecule" +;; + +let make_header kw = + let s = keyword_to_string kw in + let l = String.length s in + "\n\n"^s^"\n"^(String.init l ~f:(fun _ -> '='))^"\n\n" +;; + +let get s = + let header = (make_header s) + and rst = match s with + | Full_ci -> + Input.Full_ci.(to_rst (read ())) + | Hartree_fock -> + Input.Hartree_fock.(to_rst (read ())) + | Mo_basis -> + Input.Mo_basis.(to_rst (read ())) + | Electrons -> + Input.Electrons.(to_rst (read ())) + | Determinants -> + Input.Determinants.(to_rst (read ())) + | Cisd_sc2 -> + Input.Cisd_sc2.(to_rst (read ())) + | Nuclei -> + Input.Nuclei.(to_rst (read ())) + | Ao_basis -> + Input.Ao_basis.(to_rst (read ())) + | Bielec_integrals -> + Input.Bielec_integrals.(to_rst (read ())) + + in header^(Rst_string.to_string rst) +;; + +let set str s = + let header = (make_header s) in + let index_begin = String.substr_index_exn ~pos:0 ~pattern:header str in + let index_begin = index_begin + (String.length header) in + let index_end = + match ( String.substr_index ~pos:(index_begin+(String.length header)+1) + ~pattern:"==" str) with + | Some i -> i + | None -> String.length str + in + let l = index_end - index_begin in + let str = String.sub ~pos:index_begin ~len:l str + |> Rst_string.of_string + in + match s with + (* + | Mo_basis -> + *) + | Hartree_fock -> + Input.Hartree_fock.(write (of_rst str )) + | Full_ci -> + Input.Full_ci.(write (of_rst str)) + | Electrons -> + Input.Electrons.(write (of_rst str)) + | Determinants -> + Input.Determinants.(write (of_rst str)) + | Cisd_sc2 -> + Input.Cisd_sc2.(write (of_rst str)) + | Nuclei -> + Input.Nuclei.(write (of_rst str)) + | Bielec_integrals -> + Input.Bielec_integrals.(write (of_rst str)) + (* + | Ao_basis -> + *) + +;; + + +let create_temp_file ezfio_filename fields = + let temp_filename = Filename.temp_file "qp_edit_" ".rst" in + Out_channel.with_file temp_filename ~f:(fun out_channel -> + (file_header ezfio_filename) :: (List.map ~f:get fields) + |> String.concat ~sep:"\n" + |> Out_channel.output_string out_channel + ); + temp_filename +;; + +let run ezfio_filename = + + (* Open EZFIO *) + if (not (Sys.file_exists_exn ezfio_filename)) then + failwith (ezfio_filename^" does not exists"); + + Ezfio.set_file ezfio_filename; + + (* + let output = (file_header ezfio_filename) :: ( + List.map ~f:get [ + Ao_basis ; + Mo_basis ; + ]) + in + String.concat output + |> print_string + *) + + let tasks = [ + Nuclei ; + Electrons ; + Bielec_integrals ; + Hartree_fock ; + Cisd_sc2 ; + Full_ci ; + Determinants ; + ] + in + + (* Create the temp file *) + let temp_filename = create_temp_file ezfio_filename tasks in + + (* Open the temp file with external editor *) + let editor = + match Sys.getenv "EDITOR" with + | Some editor -> editor + | None -> "vi" + in + let command = Printf.sprintf "%s %s" editor temp_filename in + Sys.command_exn command; + + (* Re-read the temp file *) + let temp_string = + In_channel.with_file temp_filename ~f:(fun in_channel -> + In_channel.input_all in_channel) + in + List.iter ~f:(fun x -> set temp_string x) tasks; + + (* Remove temp_file *) + Sys.remove temp_filename; +;; + + + +let spec = + let open Command.Spec in + empty +(* + +> flag "i" (optional string) + ~doc:"Prints input data" + +> flag "o" (optional string) + ~doc:"Prints output data" +*) + +> anon ("ezfio_file" %: string) +;; + +let command = + Command.basic + ~summary: "Quantum Package command" + ~readme:(fun () -> + " +Edit input data + ") + spec +(* (fun i o ezfio_file () -> *) + (*fun ezfio_file () -> + try + run ezfio_file + with + | _ msg -> print_string ("\n\nError\n\n"^msg^"\n\n") + *) + (fun ezfio_file () -> run ezfio_file) +;; + +let () = + Command.run command +;; + + + diff --git a/ocaml/qp_print.ml b/ocaml/qp_print.ml index 6734b00a..6405d0ee 100644 --- a/ocaml/qp_print.ml +++ b/ocaml/qp_print.ml @@ -52,13 +52,14 @@ let run_i ~action ezfio_filename = let compute_charge () = let input = Input.Electrons.read () in - let nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data) - |> Ezfio.flattened_ezfio_data |> Array.map ~f:(Float.to_int) + let nucl_charge = Ezfio.get_nuclei_nucl_charge () + |> Ezfio.flattened_ezfio |> Array.map ~f:(Float.to_int) and n_alpha = input.Input.Electrons.elec_alpha_num |> Elec_alpha_number.to_int and n_beta = input.Input.Electrons.elec_beta_num |> Elec_beta_number.to_int in Array.fold ~init:(-n_alpha-n_beta) ~f:(fun x y -> x+y) nucl_charge + |> Charge.of_int in let compute_multiplicity () = @@ -70,17 +71,17 @@ let run_i ~action ezfio_filename = let create_molecule () = let nucl_num = Ezfio.get_nuclei_nucl_num () - and nucl_charge = Ezfio.((get_nuclei_nucl_charge ()).data) - |> Ezfio.flattened_ezfio_data - and nucl_coord = Ezfio.((get_nuclei_nucl_coord ()).data ) - |> Ezfio.flattened_ezfio_data + and nucl_charge = Ezfio.get_nuclei_nucl_charge () + |> Ezfio.flattened_ezfio + and nucl_coord = Ezfio.get_nuclei_nucl_coord () + |> Ezfio.flattened_ezfio in let nucl_label = if (Ezfio.has_nuclei_nucl_label ()) then - Ezfio.((get_nuclei_nucl_label ()).data) |> Ezfio.flattened_ezfio_data + Ezfio.get_nuclei_nucl_label () |> Ezfio.flattened_ezfio else Array.map ~f:(fun x-> x - |> Float.to_int + |> Charge.of_float |> Element.of_charge |> Element.to_string ) nucl_charge in diff --git a/ocaml/qp_set_ddci.ml b/ocaml/qp_set_ddci.ml new file mode 100644 index 00000000..48a56153 --- /dev/null +++ b/ocaml/qp_set_ddci.ml @@ -0,0 +1,303 @@ +open Qputils;; +open Qptypes;; +open Core.Std;; + +(* + * Command-line arguments + * ---------------------- + *) + +let build_mask from upto n_int = + let from = MO_number.to_int from + and upto = MO_number.to_int upto + and n_int = N_int_number.to_int n_int + in + let rec build_mask bit = function + | 0 -> [] + | i -> + if ( i = upto ) then + Bit.One::(build_mask Bit.One (i-1)) + else if ( i = from ) then + Bit.One::(build_mask Bit.Zero (i-1)) + else + bit::(build_mask bit (i-1)) + in + let starting_bit = + if ( (upto >= n_int*64) || (upto < 0) ) then Bit.One + else Bit.Zero + in + build_mask starting_bit (n_int*64) + |> List.rev +;; + + +let failure s = raise (Failure s) +;; + +type t = + | Core + | Inactive + | Active + | Virtual + | Deleted + | None +;; + +let t_to_string = function + | Core -> "core" + | Inactive -> "inactive" + | Active -> "active" + | Virtual -> "virtual" + | Deleted -> "deleted" + | None -> assert false +;; + +let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename = + + Ezfio.set_file ezfio_filename ; + if not (Ezfio.has_mo_basis_mo_tot_num ()) then + failure "mo_basis/mo_tot_num not found" ; + + let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in + let n_int = + try N_int_number.of_int (Ezfio.get_determinants_n_int ()) + with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num + in + + + let mo_class = Array.init mo_tot_num ~f:(fun i -> None) in + + (* Check input data *) + let apply_class l = + let rec apply_class t = function + | [] -> () + | k::tail -> let i = MO_number.to_int k in + begin + match mo_class.(i-1) with + | None -> mo_class.(i-1) <- t ; + apply_class t tail; + | x -> failure + (Printf.sprintf "Orbital %d is defined both in the %s and %s spaces" + i (t_to_string x) (t_to_string t)) + end + in + match l with + | MO_class.Core x -> apply_class Core x + | MO_class.Inactive x -> apply_class Inactive x + | MO_class.Active x -> apply_class Active x + | MO_class.Virtual x -> apply_class Virtual x + | MO_class.Deleted x -> apply_class Deleted x + in + + let core_input = core in + let core = MO_class.create_core core in + let inact = MO_class.create_inactive inact in + let act = MO_class.create_active act in + let virt = MO_class.create_virtual virt in + let del = MO_class.create_deleted del in + + apply_class core ; + apply_class inact ; + apply_class act ; + apply_class virt ; + apply_class del ; + + for i=1 to (Array.length mo_class) + do + if (mo_class.(i-1) = None) then + failure (Printf.sprintf "Orbital %d is not specified (mo_tot_num = %d)" i mo_tot_num) + done; + + + (* Debug output *) + MO_class.to_string core |> print_endline ; + MO_class.to_string inact |> print_endline ; + MO_class.to_string act |> print_endline ; + MO_class.to_string virt |> print_endline ; + MO_class.to_string del |> print_endline ; + + (* Create masks *) + let ia = Excitation.create_single inact act + and aa = Excitation.create_single act act + and av = Excitation.create_single act virt + and iv = Excitation.create_single inact virt + in + let single_excitations = [| ia ; aa ; av ; iv |] + |> Array.map ~f:Excitation.(fun x -> + match x with + | Single (x,y) -> + ( MO_class.to_bitlist n_int (Hole.to_mo_class x), + MO_class.to_bitlist n_int (Particle.to_mo_class y) ) + | Double _ -> assert false + ) + + and double_excitations = [| + Excitation.double_of_singles ia ia ; + Excitation.double_of_singles ia aa ; + Excitation.double_of_singles ia iv ; + Excitation.double_of_singles ia av ; + + Excitation.double_of_singles aa aa ; + Excitation.double_of_singles aa iv ; + Excitation.double_of_singles aa av ; + + Excitation.double_of_singles iv aa ; + Excitation.double_of_singles iv av ; + +(* Excitation.double_of_singles iv iv ; *) + |] + + |> Array.map ~f:Excitation.(fun x -> + match x with + | Single _ -> assert false + | Double (x,y,z,t) -> + ( MO_class.to_bitlist n_int (Hole.to_mo_class x), + MO_class.to_bitlist n_int (Particle.to_mo_class y) , + MO_class.to_bitlist n_int (Hole.to_mo_class z), + MO_class.to_bitlist n_int (Particle.to_mo_class t) ) + ) + in + + let extract_hole (h,_) = h + and extract_particle (_,p) = p + and extract_hole1 (h,_,_,_) = h + and extract_particle1 (_,p,_,_) = p + and extract_hole2 (_,_,h,_) = h + and extract_particle2 (_,_,_,p) = p + in + let result_ref = + let core = MO_class.create_inactive core_input in + let cv = Excitation.create_single core virt in + let cv = match cv with + | Excitation.Single (x,y) -> + ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x), + MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) + | Excitation.Double _ -> assert false + in + let iv = match iv with + | Excitation.Single (x,y) -> + ( MO_class.to_bitlist n_int (Excitation.Hole.to_mo_class x), + MO_class.to_bitlist n_int (Excitation.Particle.to_mo_class y) ) + | Excitation.Double _ -> assert false + in + [ Bitlist.or_operator (extract_hole iv) (extract_hole cv); + extract_particle iv ] + in + + let n_single = Array.length single_excitations in + let n_mask = Array.length double_excitations in + let zero = List.init (N_int_number.to_int n_int) ~f:(fun i -> 0L) + |> Bitlist.of_int64_list in + let result_gen = (List.init n_single ~f:(fun i-> [ + extract_hole single_excitations.(i) ; + extract_particle single_excitations.(i) ; + extract_hole1 double_excitations.(i) ; + extract_particle1 double_excitations.(i) ; + extract_hole2 double_excitations.(i) ; + extract_particle2 double_excitations.(i) ; ]) + )@(List.init (n_mask-n_single) ~f:(fun i-> [ + zero ; zero ; + extract_hole1 double_excitations.(n_single+i) ; + extract_particle1 double_excitations.(n_single+i) ; + extract_hole2 double_excitations.(n_single+i) ; + extract_particle2 double_excitations.(n_single+i) ; ]) + ) + |> List.concat + in + + (* Print bitmasks *) + print_endline "Reference Bitmasks:"; + List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_ref; + + print_endline "Generators Bitmasks:"; + List.iter ~f:(fun x-> print_endline (Bitlist.to_string x)) result_gen; + + (* Transform to int64 *) + let result_gen = List.map ~f:(fun x -> + let y = Bitlist.to_int64_list x in y@y ) + result_gen + |> List.concat + in + let result_ref = List.map ~f:(fun x -> + let y = Bitlist.to_int64_list x in y@y ) + result_ref + |> List.concat + in + + (* Write generators masks *) + Ezfio.set_bitmasks_n_int (N_int_number.to_int n_int); + Ezfio.set_bitmasks_bit_kind 8; + Ezfio.set_bitmasks_n_mask_gen n_mask; + Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 6; n_mask|]) ~data:result_gen + |> Ezfio.set_bitmasks_generators ; + + (* Write reference masks *) + Ezfio.set_bitmasks_n_mask_ref 1; + Ezfio.ezfio_array_of_list ~rank:4 ~dim:([| (N_int_number.to_int n_int) ; 2; 2; 1|]) ~data:result_ref + |> Ezfio.set_bitmasks_reference ; + + +;; + +let ezfio_file = + let failure filename = + eprintf "'%s' is not an EZFIO file.\n%!" filename; + exit 1 + in + Command.Spec.Arg_type.create + (fun filename -> + match Sys.is_directory filename with + | `Yes -> + begin + match Sys.is_file (filename ^ "/.version") with + | `Yes -> filename + | _ -> failure filename + end + | _ -> failure filename + ) +;; + +let default range = + let failure filename = + eprintf "'%s' is not a regular file.\n%!" filename; + exit 1 + in + Command.Spec.Arg_type.create + (fun filename -> + match Sys.is_directory filename with + | `Yes -> + begin + match Sys.is_file (filename^"/.version") with + | `Yes -> filename + | _ -> failure filename + end + | _ -> failure filename + ) +;; + +let spec = + let open Command.Spec in + empty + +> flag "core" (optional string) ~doc:"range Range of core orbitals" + +> flag "inact" (optional string) ~doc:"range Range of inactive orbitals" + +> flag "act" (optional string) ~doc:"range Range of active orbitals" + +> flag "virt" (optional string) ~doc:"range Range of virtual orbitals" + +> flag "del" (optional string) ~doc:"range Range of deleted orbitals" + +> anon ("ezfio_filename" %: ezfio_file) +;; + +let command = + Command.basic + ~summary: "Quantum Package command" + ~readme:(fun () -> + "Set the orbital classes in an EZFIO directory + The range of MOs has the form : \"[36-53,72-107,126-131]\" + ") + spec + (fun core inact act virt del ezfio_filename () -> run ?core ?inact ?act ?virt ?del ezfio_filename ) +;; + +let () = + Command.run command + + diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index fa3bbdd5..63a53775 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -36,40 +36,11 @@ let input_data = " * Non_empty_string : string assert (x <> \"\") ; -* MO_number : int - assert (x > 0) ; - if (x > 1000) then - warning \"More than 1000 MOs\"; - if (Ezfio.has_mo_basis_mo_tot_num ()) then - assert (x <= (Ezfio.get_mo_basis_mo_tot_num ())); -* AO_number : int - assert (x > 0) ; - if (x > 1000) then - warning \"More than 1000 AOs\"; - if (Ezfio.has_ao_basis_ao_num ()) then - assert (x <= (Ezfio.get_ao_basis_ao_num ())); - -* Nucl_number : int - assert (x > 0) ; - if (x > 1000) then - warning \"More than 1000 atoms\"; - if (Ezfio.has_nuclei_nucl_num ()) then - assert (x <= (Ezfio.get_nuclei_nucl_num ())); - -* N_int_number : int - assert (x > 0) ; - if (x > 100) then - warning \"N_int > 100\"; - if (Ezfio.has_determinants_n_int ()) then - assert (x = (Ezfio.get_determinants_n_int ())); - -* Det_number : int +* Det_number_max : int assert (x > 0) ; if (x > 100000000) then warning \"More than 100 million determinants\"; - if (Ezfio.has_determinants_det_num ()) then - assert (x <= (Ezfio.get_determinants_det_num ())); * States_number : int assert (x > 0) ; @@ -97,6 +68,9 @@ let input_data = " * MO_coef : float +* MO_occ : float + assert (x >= 0.); + * AO_coef : float * AO_expo : float @@ -121,42 +95,56 @@ let input_data = " * Elec_number : int assert (x > 0) ; +* MD5 : string + assert ((String.length x) = 32); + +* Rst_string : string + +" +;; + +let input_ezfio = " +* MO_number : int + mo_basis_mo_tot_num + 1 : 10000 + More than 10000 MOs + +* AO_number : int + ao_basis_ao_num + 1 : 10000 + More than 10000 AOs + +* Nucl_number : int + nuclei_nucl_num + 1 : 10000 + More than 10000 nuclei + +* N_int_number : int + determinants_n_int + 1 : 30 + N_int > 30 + +* Det_number : int + determinants_n_det + 1 : 100000000 + More than 100 million determinants + " ;; let untouched = " -module Determinant : sig - type t - val to_int64_array : t -> int64 array - val of_int64_array : int64 array -> t - val to_string : t -> string -end = struct - type t = int64 array - let to_int64_array x = x - let of_int64_array x = - if (Ezfio.has_determinants_n_int ()) then - begin - let n_int = Ezfio.get_determinants_n_int () in - assert ((Array.length x) = n_int*2) - end - ; x - let to_string x = Array.to_list x - |> List.map ~f:Int64.to_string - |> String.concat ~sep:\", \" -end - " let template = format_of_string " module %s : sig - type t + type t with sexp val to_%s : t -> %s - val of_%s : %s -> t + val of_%s : %s %s -> t val to_string : t -> string end = struct - type t = %s + type t = %s with sexp let to_%s x = x - let of_%s x = ( %s x ) + let of_%s %s x = ( %s x ) let to_string x = %s.to_string x end @@ -169,13 +157,18 @@ let parse_input input= | [] -> result | ( "" , "" )::tail -> parse result tail | ( t , text )::tail -> - let name , typ = String.lsplit2_exn ~on:':' t + let name,typ,params,params_val = + match String.split ~on:':' t with + | [name;typ] -> (name,typ,"","") + | name::typ::params::params_val -> (name,typ,params, + (String.concat params_val ~sep:":") ) + | _ -> assert false in let typ = String.strip typ and name = String.strip name in let typ_cap = String.capitalize typ in - let newstring = Printf.sprintf template name typ typ typ typ typ typ typ - ( String.strip text ) typ_cap + let newstring = Printf.sprintf template name typ typ typ params_val typ typ + typ typ params ( String.strip text ) typ_cap in List.rev (parse (newstring::result) tail ) in @@ -186,9 +179,76 @@ let parse_input input= |> print_string ;; -let () = - parse_input input_data ; - print_endline untouched + +let ezfio_template = format_of_string " +module %s : sig + type t with sexp + val to_%s : t -> %s + val get_max : unit -> %s + val of_%s : ?min:%s -> ?max:%s -> %s -> t + val to_string : t -> string +end = struct + type t = %s with sexp + let to_string x = %s.to_string x + let get_max () = + if (Ezfio.has_%s ()) then + Ezfio.get_%s () + else + %s + let get_min () = + %s + let to_%s x = x + let of_%s ?(min=get_min ()) ?(max=get_max ()) x = + begin + assert (x >= min) ; + if (x > %s) then + warning \"%s\"; + begin + match max with + | %s -> () + | i -> assert ( x <= i ) + end ; + x + end +end +" ;; +let parse_input_ezfio input= + let parse s = + match ( + String.split s ~on:'\n' + |> List.filter ~f:(fun x -> (String.strip x) <> "") + ) with + | [] -> "" + | a :: b :: c :: d :: [] -> + begin + let (name,typ) = String.lsplit2_exn ~on:':' a + and ezfio_func = b + and (min, max) = String.lsplit2_exn ~on:':' c + and msg = d + in + let name :: typ :: ezfio_func :: min :: max :: msg :: [] = + match (name :: typ :: ezfio_func :: min :: max :: msg :: []) with + | l -> List.map ~f:String.strip l + | _ -> assert false + in + Printf.sprintf ezfio_template + name typ typ typ typ typ typ typ typ (String.capitalize typ) + ezfio_func ezfio_func max min typ typ max msg min + end + | _ -> failwith "Error in input_ezfio" + in + String.split ~on:'*' input + |> List.map ~f:parse + |> String.concat + |> print_string +;; + +let () = + parse_input input_data ; + parse_input_ezfio input_ezfio; + print_endline untouched; + + diff --git a/ocaml/test_basis.ml b/ocaml/test_basis.ml index 1f029bf6..79492101 100644 --- a/ocaml/test_basis.ml +++ b/ocaml/test_basis.ml @@ -21,9 +21,21 @@ let test_module () = (Basis.read_element basis_channel (Nucl_number.of_int 2) Element.F) in - Long_basis.of_basis basis - |> Long_basis.to_string - |> print_endline + print_string "Long basis\n==========\n"; + let long_basis = + Long_basis.of_basis basis + in + print_endline (Long_basis.to_string long_basis); + + let short_basis = + Long_basis.to_basis long_basis + in + if (short_basis <> basis) then + print_endline "(short_basis <> basis)" + ; + print_string "Short basis\n===========\n"; + print_endline (Basis.to_string basis); + print_endline ("MD5: "^(Basis.to_md5 basis |> MD5.to_string)); ;; test_module (); diff --git a/ocaml/test_determinants.ml b/ocaml/test_determinants.ml new file mode 100644 index 00000000..b2649828 --- /dev/null +++ b/ocaml/test_determinants.ml @@ -0,0 +1,15 @@ +open Qptypes;; + +let test_module () = + let mo_tot_num = MO_number.of_int 10 in + let det = + [| 15L ; 7L |] + |> Determinant.of_int64_array + ~n_int:(N_int_number.of_int 1) + ~alpha:(Elec_alpha_number.of_int 4) + ~beta:(Elec_beta_number.of_int 3) + in + Printf.printf "%s\n" (Determinant.to_string (~mo_tot_num:mo_tot_num) det) +;; + +test_module ();; diff --git a/ocaml/test_elements.ml b/ocaml/test_elements.ml index ef650e1e..119be84b 100644 --- a/ocaml/test_elements.ml +++ b/ocaml/test_elements.ml @@ -1,6 +1,6 @@ let test_module () = let atom = Element.of_string "Cobalt" in - Printf.printf "%s %d\n" (Element.to_string atom) (Element.to_charge atom) + Printf.printf "%s %d\n" (Element.to_string atom) (Charge.to_int (Element.to_charge atom)) ;; test_module ();; diff --git a/ocaml/test_gto.ml b/ocaml/test_gto.ml index 3675f501..9b834cd2 100644 --- a/ocaml/test_gto.ml +++ b/ocaml/test_gto.ml @@ -13,11 +13,20 @@ let test_gto_1 () = let in_channel = open_in "/home/scemama/quantum_package/data/basis/cc-pvdz" in ignore (input_line in_channel); let gto = Gto.read_one in_channel in - print_string (Gto.to_string gto); - let gto = Gto.read_one in_channel in - print_string (Gto.to_string gto); - let gto = Gto.read_one in_channel in - print_string (Gto.to_string gto); + print_endline (Gto.to_string gto); + In_channel.seek in_channel 0L; + ignore (input_line in_channel); + let gto2 = Gto.read_one in_channel in + print_endline (Gto.to_string gto2); + let gto3 = Gto.read_one in_channel in + print_endline (Gto.to_string gto3); + if (gto2 = gto) then + print_endline "gto2 = gto"; + if (gto3 = gto) then + print_endline "gto3 = gto"; + if (gto3 = gto3) then + print_endline "gto3 = gto3"; + ;; let test_gto_2 () = @@ -34,7 +43,7 @@ let test_gto () = ;; let test_module () = - test_gto() + test_gto_1() ;; test_module ();; diff --git a/ocaml/test_input.ml b/ocaml/test_input.ml index 917ced2a..613f3b1a 100644 --- a/ocaml/test_input.ml +++ b/ocaml/test_input.ml @@ -1,15 +1,26 @@ +open Qptypes;; + let test_ao () = Ezfio.set_file "F2.ezfio" ; let b = Input.Ao_basis.read () in print_endline (Input.Ao_basis.to_string b); + print_endline (Input.Ao_basis.to_rst b |> Rst_string.to_string); ;; let test_bielec_intergals () = Ezfio.set_file "F2.ezfio" ; let b = Input.Bielec_integrals.read () in - print_endline (Input.Bielec_integrals.to_string b); + let output = Input.Bielec_integrals.to_string b + in + print_endline output; + let rst = Input.Bielec_integrals.to_rst b in + let b2 = Input.Bielec_integrals.of_rst rst in + if (b = b2) then + print_endline "OK" + else + print_endline "rst failed"; ;; let test_bitmasks () = @@ -30,7 +41,14 @@ let test_dets () = Ezfio.set_file "F2.ezfio" ; let b = Input.Determinants.read () in - print_endline (Input.Determinants.to_string b); + print_endline (Input.Determinants.to_rst b |> Rst_string.to_string ) ; + print_endline (Input.Determinants.sexp_of_t b |> Sexplib.Sexp.to_string ) ; + let r = Input.Determinants.to_rst b in + let b2 = Input.Determinants.of_rst r in + if (b2 = b) then + print_endline "OK" + else + print_endline "Failed" ;; let test_cisd_sc2 () = @@ -38,6 +56,13 @@ let test_cisd_sc2 () = let b = Input.Cisd_sc2.read () in print_endline (Input.Cisd_sc2.to_string b); + let rst = Input.Cisd_sc2.to_rst b in + let b2 = Input.Cisd_sc2.of_rst rst in + if (b = b2) then + print_endline "OK" + else + print_endline "rst failed"; + ;; let test_electrons () = @@ -45,6 +70,12 @@ let test_electrons () = let b = Input.Electrons.read () in print_endline (Input.Electrons.to_string b); + let rst = Input.Electrons.to_rst b in + let new_b = Input.Electrons.of_rst rst in + if (b = new_b) then + print_endline "OK" + else + print_endline "Failed in rst" ;; let test_fci () = @@ -52,6 +83,13 @@ let test_fci () = let b = Input.Full_ci.read () in print_endline (Input.Full_ci.to_string b); + let rst = Input.Full_ci.to_rst b in + let new_b = Input.Full_ci.of_rst rst in + print_endline (Input.Full_ci.to_string b); + if (b = new_b) then + print_endline "OK" + else + print_endline "Failed in rst" ;; let test_hf () = @@ -59,6 +97,13 @@ let test_hf () = let b = Input.Hartree_fock.read () in print_endline (Input.Hartree_fock.to_string b); + let rst = Input.Hartree_fock.to_rst b in + let new_b = Input.Hartree_fock.of_rst rst in + print_endline (Input.Hartree_fock.to_string b); + if (b = new_b) then + print_endline "OK" + else + print_endline "Failed in rst" ;; let test_mo () = @@ -70,19 +115,27 @@ let test_mo () = let test_nucl () = Ezfio.set_file "F2.ezfio" ; - let b = Input.Nuclei.read () - in + let b = Input.Nuclei.read () in + let rst = Input.Nuclei.to_rst b in + let new_b = Input.Nuclei.of_rst rst in print_endline (Input.Nuclei.to_string b); + if (b = new_b) then + print_endline "OK" + else + print_endline "Failed in rst" ;; (* -test_hf ();; test_ao ();; -test_bielec_intergals ();; test_bitmasks (); test_cis (); -test_dets (); test_cisd_sc2 (); +test_dets (); +test_hf ();; test_mo ();; -*) test_nucl (); +test_bielec_intergals ();; +test_electrons(); +*) +test_dets (); + diff --git a/ocaml/test_mo_label.ml b/ocaml/test_mo_label.ml new file mode 100644 index 00000000..714ca6e3 --- /dev/null +++ b/ocaml/test_mo_label.ml @@ -0,0 +1,5 @@ +let () = + let m = MO_label.of_string "canonical" in + let s = MO_label.to_string m in + print_string s +;; diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml index 7cd7e838..3abd7e9a 100644 --- a/ocaml/test_molecule.ml +++ b/ocaml/test_molecule.ml @@ -21,7 +21,7 @@ H 1.0 0.54386314 0.00000000 -0.92559535 print_string "---\n"; let m = Molecule.of_xyz_string xyz in print_endline (Molecule.name m) ; - let m = Molecule.of_xyz_string xyz ~charge:1 ~multiplicity:(Multiplicity.of_int 2) + let m = Molecule.of_xyz_string xyz ~charge:(Charge.of_int 1) ~multiplicity:(Multiplicity.of_int 2) in print_endline (Molecule.name m) ; let xyz = @@ -31,7 +31,7 @@ O 1.65102147 0.00000000 -2.35602344 H 0.54386314 0.00000000 -0.92559535 " in - let m = Molecule.of_xyz_string xyz ~charge:(-2) + let m = Molecule.of_xyz_string xyz ~charge:(Charge.of_int (-2)) in print_endline (Molecule.name m) ; print_endline (Molecule.to_string m); print_string "---------\n"; diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 2114f103..25b086e0 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -222,6 +222,9 @@ class H_apply(object): self.data["size_max"] = str(1024*128) self.data["copy_buffer"] = """ call copy_H_apply_buffer_to_wf + if (s2_eig) then + call make_s2_eigenfunction + endif SOFT_TOUCH psi_det psi_coef N_det selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0 selection_criterion = selection_criterion_min diff --git a/scripts/install_ocaml.sh b/scripts/install_ocaml.sh index 19dbc55b..fe9dde2f 100755 --- a/scripts/install_ocaml.sh +++ b/scripts/install_ocaml.sh @@ -4,16 +4,22 @@ # Thu Oct 23 21:58:40 CEST 2014 QPACKAGE_ROOT=${PWD} +PACKAGES="core cryptokit" +if [[ -f quantum_package.rc ]] +then + source quantum_package.rc +fi make -C ocaml Qptypes.ml &> /dev/null if [[ $? -ne 0 ]] then + rm -rf -- ${HOME}/ocamlbrew scripts/fetch_from_web.py "https://raw.github.com/hcarty/ocamlbrew/master/ocamlbrew-install" ocamlbrew-install.sh cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash | tee ocamlbrew_install.log grep "source " ocamlbrew_install.log | grep "etc/ocamlbrew.bashrc" >> quantum_package.rc source quantum_package.rc - echo Y | opam install core + echo Y | opam install ${PACKAGES} fi make -C ocaml Qptypes.ml diff --git a/scripts/upgrade_ezfio.sh b/scripts/upgrade_ezfio.sh new file mode 100755 index 00000000..daa51405 --- /dev/null +++ b/scripts/upgrade_ezfio.sh @@ -0,0 +1,26 @@ +#!/bin/bash +# +# Upgrades the EZFIO library from the web. +# Tue Nov 4 00:53:13 CET 2014 + +if [[ -z ${QPACKAGE_ROOT} ]] +then + print "The QPACKAGE_ROOT environment variable is not set." + print "Please reload the quantum_package.rc file." +fi + +cd -- ${QPACKAGE_ROOT} +mv -- ${QPACKAGE_ROOT}/EZFIO ${QPACKAGE_ROOT}/EZFIO.old + +make EZFIO + +if [[ $? -eq 0 ]] +then + rm -rf -- ${QPACKAGE_ROOT}/EZFIO.old + echo "Successfully updated EZFIO" +else + rm -rf -- ${QPACKAGE_ROOT}/EZFIO + mv -- ${QPACKAGE_ROOT}/EZFIO.old ${QPACKAGE_ROOT}/EZFIO + echo "Failed to update EZFIO" +fi + diff --git a/scripts/upgrade_irpf90.sh b/scripts/upgrade_irpf90.sh new file mode 100755 index 00000000..a3aeac48 --- /dev/null +++ b/scripts/upgrade_irpf90.sh @@ -0,0 +1,25 @@ +#!/bin/bash +# +# Upgrades IRPF90 from the web. +# Tue Nov 4 00:53:13 CET 2014 + +if [[ -z ${QPACKAGE_ROOT} ]] +then + print "The QPACKAGE_ROOT environment variable is not set." + print "Please reload the quantum_package.rc file." +fi + +cd -- ${QPACKAGE_ROOT} +mv -- ${QPACKAGE_ROOT}/irpf90 ${QPACKAGE_ROOT}/irpf90.old + +make irpf90 + +if [[ $? -eq 0 ]] +then + rm -rf -- ${QPACKAGE_ROOT}/irpf90.old + echo "Successfully updated IRPF90" +else + rm -rf -- ${QPACKAGE_ROOT}/irpf90 + mv -- ${QPACKAGE_ROOT}/irpf90.old ${QPACKAGE_ROOT}/irpf90 + echo "Failed to update IRPF90" +fi diff --git a/setup_environment.sh b/setup_environment.sh index 028b4f2a..bc300076 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -53,7 +53,7 @@ echo $RED " To complete the installation, add the following line to your ~/.bashrc: -source quantum_package.rc +source ${QPACKAGE_ROOT}/quantum_package.rc ======================================================= " $BLACK diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index b6e0295d..ff0a83f6 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -91,6 +91,17 @@ subroutine bielec_integrals_index_reverse(i,j,k,l,i1) endif enddo enddo + do ii=1,8 + if (i(ii) /= 0) then + 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) + stop 'bielec_integrals_index_reverse failed' + endif + endif + enddo + end diff --git a/src/Bitmask/bitmasks.ezfio_config b/src/Bitmask/bitmasks.ezfio_config index 79f1ae6c..988510fb 100644 --- a/src/Bitmask/bitmasks.ezfio_config +++ b/src/Bitmask/bitmasks.ezfio_config @@ -4,5 +4,5 @@ bitmasks N_mask_gen integer generators integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen) N_mask_ref integer - reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_ref) + reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,2,bitmasks_N_mask_ref) diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 353d6a3b..dbf3304c 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -175,10 +175,6 @@ BEGIN_PROVIDER [ integer(bit_kind), reference_bitmask, (N_int,2,2,N_reference_bi else reference_bitmask(:,:,s_hole ,1) = HF_bitmask reference_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - reference_bitmask(:,:,d_hole1,1) = HF_bitmask - reference_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) - reference_bitmask(:,:,d_hole2,1) = HF_bitmask - reference_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) endif END_PROVIDER diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index a16698ec..801d00a5 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -1,4 +1,5 @@ use bitmasks +use omp_lib type H_apply_buffer_type integer :: N_det @@ -11,7 +12,8 @@ end type H_apply_buffer_type type(H_apply_buffer_type), pointer :: H_apply_buffer(:) -BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] + BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] +&BEGIN_PROVIDER [ integer(omp_lock_kind), H_apply_buffer_lock, (64,0:nproc-1) ] use omp_lib implicit none BEGIN_DOC @@ -24,7 +26,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] allocate(H_apply_buffer(0:nproc-1)) iproc = 0 !$OMP PARALLEL PRIVATE(iproc) DEFAULT(NONE) & - !$OMP SHARED(H_apply_buffer,N_int,sze,N_states) + !$OMP SHARED(H_apply_buffer,N_int,sze,N_states,H_apply_buffer_lock) !$ iproc = omp_get_thread_num() H_apply_buffer(iproc)%N_det = 0 H_apply_buffer(iproc)%sze = sze @@ -36,6 +38,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] H_apply_buffer(iproc)%det = 0_bit_kind H_apply_buffer(iproc)%coef = 0.d0 H_apply_buffer(iproc)%e2 = 0.d0 + call omp_init_lock(H_apply_buffer_lock(1,iproc)) !$OMP END PARALLEL endif @@ -56,6 +59,7 @@ subroutine resize_H_apply_buffer(new_size,iproc) ASSERT (iproc >= 0) ASSERT (iproc < nproc) + call omp_set_lock(H_apply_buffer_lock(1,iproc)) allocate ( buffer_det(N_int,2,new_size), & buffer_coef(new_size,N_states), & buffer_e2(new_size,N_states) ) @@ -89,6 +93,7 @@ subroutine resize_H_apply_buffer(new_size,iproc) H_apply_buffer(iproc)%sze = new_size H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det) + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) end @@ -174,6 +179,7 @@ subroutine copy_H_apply_buffer_to_wf H_apply_buffer(j)%N_det = 0 !$OMP END PARALLEL call normalize(psi_coef,N_det) + SOFT_TOUCH N_det psi_det psi_coef end @@ -194,6 +200,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) if (new_size > H_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc) endif + call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) @@ -216,6 +223,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) enddo + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) end diff --git a/src/Dets/utils.irp.f b/src/Dets/utils.irp.f index 27022965..22faee83 100644 --- a/src/Dets/utils.irp.f +++ b/src/Dets/utils.irp.f @@ -1,7 +1,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] implicit none BEGIN_DOC - ! H matrix on the basis of the slater deter;inants defined by psi_det + ! H matrix on the basis of the slater determinants defined by psi_det END_DOC integer :: i,j double precision :: hij diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index 2c823baf..5e144cc6 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -3,7 +3,7 @@ program H_CORE_guess character*(64) :: label mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef - label = "H_CORE_GUESS" + label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),label) call save_mos diff --git a/src/Makefile.common b/src/Makefile.common index 0aa47245..5b6a53ce 100644 --- a/src/Makefile.common +++ b/src/Makefile.common @@ -1,3 +1,9 @@ +# Required EZFIO version +EZFIO_VERSION=1.1 + +# Required IRPF90 version +IRPF90_VERSION=1.4 + # Check if QPACKAGE_ROOT is defined ifndef QPACKAGE_ROOT @@ -28,10 +34,12 @@ include $(QPACKAGE_ROOT)/src/Makefile.config # Check the version of IRPF90 -IRP_VERSION_OK=$(shell $(IRPF90) -v | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= 1.3") +IRP_VERSION_OK=$(shell $(IRPF90) -v | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= $(IRPF90_VERSION)") ifeq ($(IRP_VERSION_OK),False) $(info -------------------- Error --------------------) -$(info IRPF90 version >= 1.3 is required) +$(info IRPF90 version >= $(IRPF90_VERSION) is required) +$(info To upgrade IRPF90, run : ) +$(info $(QPACKAGE_ROOT)/scripts/upgrade_irpf90.sh ) $(info -----------------------------------------------) $(error ) endif @@ -85,12 +93,26 @@ $(info -----------------------------------------------) endif -# Define the Makefile common variables and rules +# Define the Makefile common variables EZFIO_DIR=$(QPACKAGE_ROOT)/EZFIO EZFIO=$(EZFIO_DIR)/lib/libezfio_irp.a INCLUDE_DIRS=$(NEEDED_MODULES) include +# Check EZFIO version + +EZFIO_VERSION_OK=$(shell cat $(EZFIO_DIR)/version | cut -d '=' -f 2 | python -c "import sys ; print float(sys.stdin.readline().rsplit('.',1)[0]) >= $(EZFIO_VERSION)") +ifeq ($(EZFIO_VERSION_OK),False) +$(info -------------------- Error --------------------) +$(info EZFIO version >= $(EZFIO_VERSION) is required ) +$(info To upgrade EZFIO, run : ) +$(info $(QPACKAGE_ROOT)/scripts/upgrade_ezfio.sh ) +$(info -----------------------------------------------) +$(error ) +endif + +# Define the EZFIO rules + $(EZFIO): $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ROOT)/src/*/*.ezfio_config) @echo Building EZFIO library @cp $(wildcard $(QPACKAGE_ROOT)/src/*.ezfio_config) $(wildcard $(QPACKAGE_ROOT)/src/*/*.ezfio_config) $(EZFIO_DIR)/config diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 40e3f67c..8751b818 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -26,6 +26,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c if (new_size > h_apply_buffer(iproc)%sze) then call resize_h_apply_buffer(max(h_apply_buffer(iproc)%sze*2,new_size),iproc) endif + call omp_set_lock(H_apply_buffer_lock(1,iproc)) do i=1,H_apply_buffer(iproc)%N_det ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) @@ -63,6 +64,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) enddo + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) !$OMP CRITICAL selection_criterion = max(selection_criterion,smax) selection_criterion_min = min(selection_criterion_min,smin) @@ -146,6 +148,7 @@ subroutine make_s2_eigenfunction integer, parameter :: bufsze = 1000 logical, external :: is_in_wavefunction + print *, irp_here ! !TODO DEBUG ! do i=1,N_det ! do j=i+1,N_det @@ -172,9 +175,11 @@ subroutine make_s2_eigenfunction do i=1,N_occ_pattern call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int) + s += 1 if (s > smax) then deallocate(d) allocate ( d(N_int,2,s) ) + smax = s endif call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) do j=1,s