diff --git a/bin/qp_set_frozen_core b/bin/qp_set_frozen_core index f9761144..d2821bd9 100755 --- a/bin/qp_set_frozen_core +++ b/bin/qp_set_frozen_core @@ -83,6 +83,7 @@ def main(arguments): elif charge <= 118: n_frozen += 43 elif arguments["--small"]: + for charge in ezfio.nuclei_nucl_charge: if charge <= 4: pass elif charge <= 18: n_frozen += 1 elif charge <= 36: n_frozen += 5 diff --git a/docs/source/appendix/contributors.rst b/docs/source/appendix/contributors.rst index e3574e5a..74837282 100644 --- a/docs/source/appendix/contributors.rst +++ b/docs/source/appendix/contributors.rst @@ -46,6 +46,7 @@ The following people have contributed to this project (by alphabetical order): * Nicolas Renon * Lorenzo Tenti * Julien Toulouse +* Diata Traoré * Mikaël Véril diff --git a/docs/source/index.rst b/docs/source/index.rst index 4231b1f8..273582d4 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -39,9 +39,10 @@ programmers_guide/programming programmers_guide/ezfio programmers_guide/plugins + programmers_guide/plugins_tuto_intro + programmers_guide/plugins_tuto_I programmers_guide/new_ks programmers_guide/index - programmers_guide/plugins .. toctree:: @@ -52,5 +53,6 @@ appendix/benchmarks appendix/license appendix/contributors + appendix/references diff --git a/docs/source/programmers_guide/plugins_tuto_I.rst b/docs/source/programmers_guide/plugins_tuto_I.rst new file mode 100644 index 00000000..27864487 --- /dev/null +++ b/docs/source/programmers_guide/plugins_tuto_I.rst @@ -0,0 +1 @@ +.. include:: ../../../plugins/local/tuto_plugins/tuto_I/tuto_I.rst diff --git a/docs/source/programmers_guide/plugins_tuto_intro.rst b/docs/source/programmers_guide/plugins_tuto_intro.rst new file mode 100644 index 00000000..63482462 --- /dev/null +++ b/docs/source/programmers_guide/plugins_tuto_intro.rst @@ -0,0 +1 @@ +.. include:: ../../../plugins/README.rst diff --git a/external/irpf90 b/external/irpf90 index ba1a2837..4ab1b175 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit ba1a2837aa61cb8f9892860cec544d7c6659badd +Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 diff --git a/plugins/README.rst b/plugins/README.rst new file mode 100644 index 00000000..3dc50873 --- /dev/null +++ b/plugins/README.rst @@ -0,0 +1,131 @@ +============================== +Tutorial for creating a plugin +============================== + +Introduction: what is a plugin, and what tutorial will be about ? +================================================================= + +The |QP| is split into two kinds of routines/global variables (i.e. *providers*): + 1) the **core modules** locatedin qp2/src/, which contains all the bulk of a quantum chemistry software (integrals, matrix elements between Slater determinants, linear algebra routines, DFT stuffs etc..) + 2) the **plugins** which are external routines/*providers* connected to the qp2/src/ routines/*providers*. + +More precisely, a **plugin** of the |QP| is a directory where you can create routines, +providers and executables that use all the global variables/functions/routines already created +in the modules of qp2/src or in other plugins. + +Instead of giving a theoretical lecture on what is a plugin, +we will go through a series of examples that allow you to do the following thing: + +1) print out **one- and two-electron integrals** on the AO/MO basis, creates two providers which manipulate these objects, print out these providers, + +2) browse the **Slater determinants stored** in the |EZFIO| wave function and compute their matrix elements, + +3) build the **Hamiltonian matrix** and **diagonalize** it either with **Lapack or Davidson**, + +4) print out the **one- and two-electron rdms**, + +5) obtain the **AOs** and **MOs** on the **DFT grid**, together with the **density**, + +How the tutorial will be done +----------------------------- + +This tuto is as follows: + + 1) you **READ THIS FILE UNTIL THE END** in order to get the big picture and vocabulary, + + 2) you go to the directory :file:`qp2/plugins/tuto_plugins/` and you will find detailed tutorials for each of the 5 examples. + +Creating a plugin: the basic +---------------------------- + +The first thing to do is to be in the QPSH mode: you execute the qp2/bin/qpsh script that essentially loads all +the environement variables and allows for the completion of command lines in bash (that is an AMAZING feature :) + +Then, you need to known **where** you want to create your plugin, and what is the **name** of the plugin. + +.. important:: + + The plugins are **NECESSARILY** located in qp2/plugins/, and from there you can create any structures of directories. + + +Ex: If you want to create a plugin named "my_fancy_plugin" in the directory plugins/plugins_test/, +this goes with the command + +.. code:: bash + + qp plugins create -n my_fancy_plugin -r plugins_test/ + +Then, to create the plugin of your dreams, the two questions you need to answer are the following: + +1) What do I **need** to compute what I want, which means what are the **objects** that I need ? + + There are two kind of objects: + + + the *routines/functions*: + + Ex: Linear algebra routines, integration routines etc ... + + + the global variables which are called the *providers*: + + Ex: one-electron integrals, Slater determinants, density matrices etc ... + +2) **Where do I find** these objects ? + + The objects (routines/functions/providers) are necessarily created in other *modules/plugins*. + +.. seealso:: + + The routine :c:func:`lapack_diagd` (which diagonalises a real hermitian matrix) is located in the file + :file:`qp2/src/utils/linear_algebra.irp.f` + therefore it "belongs" to the module :ref:`module_utils` + + The routine :c:func:`ao_to_mo` (which converts a given matrix A from the AO basis to the MO basis) is located in the file + :file:`qp2/src/mo_one_e_ints/ao_to_mo.irp.f` + therefore it "belongs" to the module :ref:`module_mo_one_e_ints` + + The provider :c:data:`ao_one_e_integrals` (which is the integrals of one-body part of H on the AO basis) is located in the file + :file:`qp2/src/ao_one_e_ints/ao_one_e_ints.irp.f` + therefore it belongs to the module :ref:`module_ao_one_e_ints` + + The provider :c:data:`one_e_dm_mo_beta_average` (which is the state average beta density matrix on the MO basis) is located in the file + :file:`qp2/src/determinants/density_matrix.irp.f` + therefore it belongs to the module :ref:`module_determinants` + +To import all the variables that you need, you just need to write the name of the plugins in the :file:`NEED` file . + +To import all the variables/routines of the module :ref:`module_utils`, :ref:`module_determinants` and :ref:`module_mo_one_e_ints`, the :file:`NEED` file you will need is simply the following: + +.. code:: bash + + cat NEED + + utils + determinants + mo_one_e_ints + + +.. important:: + + There are **many** routines/providers in the core modules of QP. + + Nevertheless, as everything is coded with the |IRPF90|, you can use the following amazing tools: :command:`irpman` + + :command:`irpman` can be used in command line in bash to obtain all the info on a routine or variable ! + + +Example: execute the following command line : + +.. code:: bash + + irpman ao_one_e_integrals + +Then all the information you need on :c:data:`ao_one_e_integrals` will appear on the screen. +This includes + + - **where** the provider is created, (*i.e.* the actual file where the provider is designed) + - the **type** of the provider (*i.e.* a logical, integer etc ...) + - the **dimension** if it is an array, + - what other *providers* are **needed** to build this provider, + - what other *providers* **need** this provider. + + diff --git a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f index e0637fa5..12163e06 100644 --- a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f @@ -960,7 +960,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d ! endif e_pert(istate) = 0.25 * val / delta_E ! e_pert(istate) = 0.5d0 * (tmp - delta_E) - if(dsqrt(dabs(tmp)).gt.1.d-4.and.dabs(alpha_h_psi).gt.1.d-4)then + if(dsqrt(tmp).gt.1.d-4.and.dabs(psi_h_alpha).gt.1.d-4)then coef(istate) = e_pert(istate) / psi_h_alpha else coef(istate) = alpha_h_psi / delta_E diff --git a/plugins/local/tuto_plugins/H2.xyz b/plugins/local/tuto_plugins/H2.xyz new file mode 100644 index 00000000..7af12291 --- /dev/null +++ b/plugins/local/tuto_plugins/H2.xyz @@ -0,0 +1,6 @@ +2 +H2, equilibrium geometry +H 0.0 0.0 0. +H 0.0 0.0 0.74 + + diff --git a/plugins/local/tuto_plugins/n2.xyz b/plugins/local/tuto_plugins/n2.xyz new file mode 100644 index 00000000..016732d8 --- /dev/null +++ b/plugins/local/tuto_plugins/n2.xyz @@ -0,0 +1,4 @@ +2 +N2 Geo: Experiment Mult: 1 symmetry: 14 +N 0.0 0.0 0.5488 +N 0.0 0.0 -0.5488 diff --git a/plugins/local/tuto_plugins/tuto_I/print_one_e_h.irp.f b/plugins/local/tuto_plugins/tuto_I/print_one_e_h.irp.f new file mode 100644 index 00000000..5d8dc1e7 --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/print_one_e_h.irp.f @@ -0,0 +1,20 @@ +program my_program_to_print_stuffs + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + integer :: i,j + print*,'AO integrals ' + do i = 1, ao_num + do j = 1, ao_num + print*,j,i,ao_one_e_integrals(j,i) + enddo + enddo + + print*,'MO integrals ' + do i = 1, mo_num + do j = 1, mo_num + print*,j,i,mo_one_e_integrals(j,i) + enddo + enddo +end diff --git a/plugins/local/tuto_plugins/tuto_I/print_traces_on_e.irp.f b/plugins/local/tuto_plugins/tuto_I/print_traces_on_e.irp.f new file mode 100644 index 00000000..2bf3b86b --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/print_traces_on_e.irp.f @@ -0,0 +1,24 @@ +program my_program + implicit none + BEGIN_DOC +! This program is there essentially to show how one can use providers in programs + END_DOC + integer :: i,j + double precision :: accu + print*,'Trace on the AO basis ' + print*,trace_ao_one_e_ints + print*,'Trace on the AO basis after projection on the MO basis' + print*,trace_ao_one_e_ints_from_mo + print*,'Trace of MO integrals ' + print*,trace_mo_one_e_ints + print*,'ao_num = ',ao_num + print*,'mo_num = ',mo_num + if(ao_num .ne. mo_num)then + print*,'The AO basis and MO basis are different ...' + print*,'Trace on the AO basis should not be the same as Trace of MO integrals' + print*,'Only the second one must be equal to the trace on the MO integrals' + else + print*,'The AO basis and MO basis are the same !' + print*,'All traces should coincide ' + endif +end diff --git a/plugins/local/tuto_plugins/tuto_I/print_two_e_h.irp.f b/plugins/local/tuto_plugins/tuto_I/print_two_e_h.irp.f new file mode 100644 index 00000000..eaeb6c98 --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/print_two_e_h.irp.f @@ -0,0 +1,32 @@ +program my_program_to_print_stuffs + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + integer :: i,j,k,l + double precision :: integral + double precision :: get_ao_two_e_integral, get_two_e_integral ! declaration of the functions + print*,'AO integrals, physicist notations : ' + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + print*,i,j,k,l,integral + enddo + enddo + enddo + enddo + + print*,'MO integrals, physicist notations : ' + do i = 1, mo_num + do j = 1, mo_num + do k = 1, mo_num + do l = 1, mo_num + integral = get_two_e_integral(i, j, k, l, mo_integrals_map) + print*,i,j,k,l,integral + enddo + enddo + enddo + enddo +end diff --git a/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f b/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f new file mode 100644 index 00000000..e71d49fc --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f @@ -0,0 +1,111 @@ + +! This file is an example of the kind of manipulations that you can do with providers +! + +!!!!!!!!!!!!!!!!!!!!!!!!!! Main providers useful for the program !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!!! type name +BEGIN_PROVIDER [ double precision, trace_mo_one_e_ints] + implicit none + BEGIN_DOC +! trace_mo_one_e_ints = Trace of the one-electron integrals on the MO basis +! +! = sum_i mo_one_e_integrals(i,i) + END_DOC + integer :: i + trace_mo_one_e_ints = 0.d0 + do i = 1, mo_num + trace_mo_one_e_ints += mo_one_e_integrals(i,i) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, trace_ao_one_e_ints] + implicit none + BEGIN_DOC +! trace_ao_one_e_ints = Trace of the one-electron integrals on the AO basis taking into account the non orthogonality +! +! Be aware that the trace of an operator in a non orthonormal basis is Tr(A S^{-1}) = \sum_{m,n}(A_mn S^{-1}_mn) +! +! WARNING: it is equal to the trace on the MO basis if and only if the AO basis and MO basis +! have the same number of functions + END_DOC + integer :: i,j + double precision, allocatable :: inv_overlap_times_integrals(:,:) ! = h S^{-1} + allocate(inv_overlap_times_integrals(ao_num,ao_num)) + ! routine that computes the product of two matrices, you can check it with + ! irpman get_AB_prod + call get_AB_prod(ao_one_e_integrals,ao_num,ao_num,s_inv,ao_num,inv_overlap_times_integrals) + ! Tr(inv_overlap_times_integrals) = Tr(h S^{-1}) + trace_ao_one_e_ints = 0.d0 + do i = 1, ao_num + trace_ao_one_e_ints += inv_overlap_times_integrals(i,i) + enddo + ! + ! testing the formula Tr(A S^{-1}) = \sum_{m,n}(A_mn S^{-1}_mn) + double precision :: test + test = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + test += ao_one_e_integrals(j,i) * s_inv(i,j) + enddo + enddo + if(dabs(accu - trace_ao_one_e_ints).gt.1.d-12)then + print*,'Warning ! ' + print*,'Something is wrong because Tr(AB) \ne sum_{mn}A_mn B_nm' + endif +END_PROVIDER + +BEGIN_PROVIDER [ double precision, trace_ao_one_e_ints_from_mo] + implicit none + BEGIN_DOC +! trace_ao_one_e_ints_from_mo = Trace of the one-electron integrals on the AO basis after projection on the MO basis +! +! = Tr([SC h {SC}^+] S^{-1}) +! +! = Be aware that the trace of an operator in a non orthonormal basis is = Tr(A S^{-1}) where S is the metric +! Must be equal to the trace_mo_one_e_ints + END_DOC + integer :: i + double precision, allocatable :: inv_overlap_times_integrals(:,:) + allocate(inv_overlap_times_integrals(ao_num,ao_num)) + ! Using the provider ao_one_e_integrals_from_mo = [SC h {SC}^+] + call get_AB_prod(ao_one_e_integrals_from_mo,ao_num,ao_num,s_inv,ao_num,inv_overlap_times_integrals) + ! inv_overlap_times_integrals = [SC h {SC}^+] S^{-1} + trace_ao_one_e_ints_from_mo = 0.d0 + ! Computing the trace + do i = 1, ao_num + trace_ao_one_e_ints_from_mo += inv_overlap_times_integrals(i,i) + enddo +END_PROVIDER + +!!!!!!!!!!!!!!!!!!!!!!!!!!! Additional providers to check some stuffs !!!!!!!!!!!!!!!!!!!!!!!!! + +BEGIN_PROVIDER [ double precision, ao_one_e_int_no_ov_from_mo, (ao_num, ao_num) ] + BEGIN_DOC + ! ao_one_e_int_no_ov_from_mo = C mo_one_e_integrals C^T + ! + ! WARNING : NON EQUAL TO ao_one_e_integrals due to the non orthogonality + END_DOC + call mo_to_ao_no_overlap(mo_one_e_integrals,mo_num,ao_one_e_int_no_ov_from_mo,ao_num) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_one_e_int_no_ov_from_mo_ov_ov, (ao_num, ao_num)] + BEGIN_DOC + ! ao_one_e_int_no_ov_from_mo_ov_ov = S ao_one_e_int_no_ov_from_mo S = SC mo_one_e_integrals (SC)^T + ! + ! EQUAL TO ao_one_e_integrals ONLY IF ao_num = mo_num + END_DOC + double precision, allocatable :: tmp(:,:) + allocate(tmp(ao_num, ao_num)) + call get_AB_prod(ao_overlap,ao_num,ao_num,ao_one_e_int_no_ov_from_mo,ao_num,tmp) + call get_AB_prod(tmp,ao_num,ao_num,ao_overlap,ao_num,ao_one_e_int_no_ov_from_mo_ov_ov) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, c_t_s_c, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! C^T S C = should be the identity + END_DOC + call get_AB_prod(mo_coef_transp,mo_num,ao_num,S_mo_coef,mo_num,c_t_s_c) +END_PROVIDER + diff --git a/plugins/local/tuto_plugins/tuto_I/tuto_I.rst b/plugins/local/tuto_plugins/tuto_I/tuto_I.rst new file mode 100644 index 00000000..43b4af0b --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/tuto_I.rst @@ -0,0 +1,218 @@ +============================================= +Tuto I: One- and two-e integrals (20 minutes) +============================================= + +Requirements +------------ +1) You know how to create an |EZFIO| file and run calculations with |QP| (check the tuto: ``_), + +2) You have an |EZFIO| file with MOs created (with the :ref:`scf` executable for instance). As we are going to print out some integrals, don't take a too large system/basis (Ex: H2, cc-pVDZ is ok :) + +3) You made an qp set_file YOUR_EZFIO_FILE_FOR_H2 in order to work on that ezfio folder. + +4) You have READ the :file:`qp2/plugins/README.rst` file to HAVE THE **VOCABULARY**. + +Our goals: +---------- +We want to create a plugin to do the following things: + 1) print out one- and two-electron integrals on the AO/MO basis, + + 2) creates two providers which manipulate these objects, + + 3) print out these providers. + +I) Getting started: creating the plugin +--------------------------------------- +We will go step-by-step through these plugins. + +We will create a plugin named "plugin_I", and its location will be in "tuto_plugins". +Therefore to create the plugin, we do: + +.. code:: bash + + qp plugins create -n plugin_I -r tuto_plugins + +Then do an "ls" in qp2/plugins/tuto_plugins/ and you will find a directory called "plugin_I". + +In that directory you will find: + +1) a :file:`NEED` file that will eventually contain all the other modules/plugins needed by our "plugin_I", + +2) a :file:`README.rst` file that you can and **SHOULD** modify in order to **DOCUMENT** what is doing the plugin, + +3) a :file:`plugin_I.irp.f` file that is a program to be compiled and just printing "Hello world" + +II) Specifying the dependencies +------------------------------- +The next step is to know what are the other modules/plugins that we need to do what we want. +We need here + +a) the one-electron integrals on the AO basis, which are computed in :file:`qp2/src/ao_one_e_ints/` + +b) the one-electron integrals on the MO basis, which are computed in :file:`qp2/src/mo_one_e_ints/` + +c) the two-electron integrals on the AO basis, which are computed in :file:`qp2/src/ao_two_e_ints/` + +d) the two-electron integrals on the MO basis, which are computed in :file:`qp2/src/mo_two_e_ints/` + +Therefore, we will need the following four modules: + + a) ao_one_e_ints + b) mo_one_e_ints + c) ao_two_e_ints + d) mo_two_e_ints + +You can then create the following "NEED" file by executing the following command + +.. code:: bash + + cat < NEED + ao_one_e_ints + mo_one_e_ints + ao_two_e_ints + mo_two_e_ints + EOF + +II) Installing the plugin +------------------------- +Now that we have specified the various depenencies we need now to INSTALL the plugin, which means to create the equivalent of a Makefile for the compilation. + +To do it we simply do + +.. code:: bash + + qp plugins install plugin_I + + +III) Compiling the void plugin +------------------------------ +It is customary to compile first your "void" plugin, void in the sense that it does not contain anything else than the program printing "Hello world". + +To do so, just go in the plugin and execute the following command + +.. code:: bash + + ninja + +It does a lot of stuffs, but it must conclude with something like + +.. code:: bash + + make: Leaving directory 'SOME_PATH_TOWARD_YOUR_QP2_DIRECTORY/qp2/ocaml' + + +Since that it has compiled, an executable "plugin_I" has been created. + +Also, if you make "ls" in the "plugin_I" you will notice that many symbolink links have been created, and among which the four modules that you included in the NEED file. + +All the other modules (Ex::ref:`module_ao_basis`, :ref:`module_utils`) are here because they are need by some of the four modules that you need. +The variables that we need are + +:data:`ao_one_e_integrals` + +:data:`mo_one_e_integrals` + +You can check them with + +.. code:: bash + + irpman ao_one_e_integrals + + +.. code:: bash + + irpman mo_one_e_integrals + +in order to get some information on where they are created, and many more information. +We will now create an executable such that it prints out the integrals. + + +IV) Printing out the one-electron integrals +-------------------------------------------- +We will now create a program that will print out the one-electron integrals on the AO and MO basis. + +You can then copy the file :file:`qp2/plugins/tuto_plugins/tuto_I/print_one_e_h.irp.f` in your plugin. + +In this file you will see that we simply browse the two arrays :data:`ao_one_e_integrals` and :data:`mo_one_e_integrals`, which are the providers and we browse them until either :data:`ao_num` or :data:`mo_num` which are also providers representing the number of AOs or MOs. + + +.. seealso:: + + You can check these variables with :command:`irpman` ! + +If you recompile using |ninja| as before, and another executable has been created "print_one_e_h". +Then, you can run the program on the ezfio file by doing + +.. code:: bash + + qp run print_one_e_h + +and will print out the data you need :) + +By the way, as the file :file:`plugin_I.irp.f` contains nothing but a "Hello world" print, you can simply remove it if you want. + +V) Printing out the two-electron integrals +------------------------------------------ +We will now create a file that prints out the two-electron integrals in the AO and MO basis. +These can be accessed with the following subroutines : + +1- :c:func:`get_ao_two_e_integral` for the AO basis + +2- :c:func:`get_two_e_integral` for the MO basis + + +.. seealso:: + + check them with irpman ! + +To print the two-electron integrals, you can copy the file :file:`qp2/plugins/tuto_plugins/tuto_I/print_two_e_h.irp.f` in your plugin and recompile with |ninja|. +Then just run the program + +.. code:: bash + + qp run print_two_e_h + +and it will print all the things you want :) + +VI) Creating new providers and a program to print them +------------------------------------------------------ +We will now create new providers that manipulates the objects that we just printed. +As an example, we will compute the trace of the one electron integrals in the AO and MO basis. +In the file :file:`qp2/plugins/tuto_plugins/tuto_I/traces_one_e.irp.f` you will find the several new providers among which + + 1- :c:data:`trace_mo_one_e_ints` : simply the sum of the diagonal matrix element of the one-electron integrals + + 2- :c:data:`trace_ao_one_e_ints` : the corresponding trace on the AO basis + .. math:: + + \text{Tr}({\bf h}{\bf S}^{-1}) = \sum_{m,n} S^{-1}_{mn} h_{mn} + + + 3- :c:data:`trace_ao_one_e_ints_from_mo` : the trace on the AO basis with the integrals obtained first from the MO basis + .. math:: + + \text{Tr}({\bf \tilde{h}}{\bf S}^{-1}) = \text{Tr}\big({\bf SC h}({\bf SC }^T){\bf S}^{-1}\big) + +Just copy the :file:`qp2/plugins/tuto_plugins/tuto_I/traces_one_e.irp.f` in your plugin and recompile. + +.. seealso:: + + Once it has compiled, check your new providers with :command:`irpman` ! + +As explained in the files :file:`qp2/plugins/tuto_plugins/tuto_I/traces_one_e.irp.f` and :file:`qp2/plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f`, :c:data:`trace_mo_one_e_ints` is equal to :c:data:`trace_ao_one_e_ints` only if the number of AO basis functions is equal to the number of MO basis functions, which means if you work with cartesian functions. + + +.. seealso:: + + You can check with :command:`qp create_ezfio -h` for the option to create an |EZFIO| with cartesian basis functions + +In the file :file:`qp2/plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f` you will find an example of executable that prints out the various providers. +Copy these two files in your plugin and recompile to execute it. + +Execute the program print_traces_on_e and check for the results with + +.. code:: bash + + qp run print_traces_on_e + +The code in :file:`qp2/plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f` should be easy to read, I let the reader interpret it. diff --git a/src/ao_one_e_ints/ao_one_e_ints.irp.f b/src/ao_one_e_ints/ao_one_e_ints.irp.f index 65981dc9..9b914dee 100644 --- a/src/ao_one_e_ints/ao_one_e_ints.irp.f +++ b/src/ao_one_e_ints/ao_one_e_ints.irp.f @@ -45,3 +45,13 @@ BEGIN_PROVIDER [ double precision, ao_one_e_integrals_imag,(ao_num,ao_num)] END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_one_e_integrals_from_mo, (ao_num, ao_num)] + implicit none + BEGIN_DOC +! Integrals of the one e hamiltonian obtained from the integrals on the MO basis +! +! WARNING : this is equal to ao_one_e_integrals only if the AO and MO basis have the same number of functions + END_DOC + call mo_to_ao(mo_one_e_integrals,mo_num,ao_one_e_integrals_from_mo,ao_num) +END_PROVIDER diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 1ead9d78..3513f215 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -522,6 +522,84 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ enddo endif + if (state_following) then + if (.not. only_expected_s2) then + print*,'' + print*,'!!! State following only available with only_expected_s2 = .True. !!!' + STOP + endif + endif + + if (state_following) then + + integer :: state(N_st), idx + double precision :: omax + logical :: used + logical, allocatable :: ok(:) + double precision, allocatable :: overlp(:,:) + + allocate(overlp(shift2,N_st),ok(shift2)) + + overlp = 0d0 + do j = 1, shift2-1, N_st_diag + + ! Computes some states from the guess vectors + ! Psi(:,j:j+N_st_diag) = U y(:,j:j+N_st_diag) and put them + ! in U(1,shift2+1:shift2+1+N_st_diag) as temporary array + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y(1,j), size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + + ! Overlap + do l = 1, N_st + do k = 1, N_st_diag + do i = 1, sze + overlp(k+j-1,l) += U(i,l) * U(i,shift2+k) + enddo + enddo + enddo + + enddo + + state = 0 + do l = 1, N_st + + omax = 0d0 + idx = 0 + do k = 1, shift2 + + ! Already used ? + used = .False. + do i = 1, N_st + if (state(i) == k) then + used = .True. + endif + enddo + + ! Maximum overlap + if (dabs(overlp(k,l)) > omax .and. .not. used .and. state_ok(k)) then + omax = dabs(overlp(k,l)) + idx = k + endif + enddo + + state(l) = idx + enddo + + ! tmp array before setting state_ok + ok = .False. + do l = 1, N_st + ok(state(l)) = .True. + enddo + + do k = 1, shift2 + if (.not. ok(k)) then + state_ok(k) = .False. + endif + enddo + + deallocate(overlp,ok) + endif + do k=1,shift2 if (.not. state_ok(k)) then do l=k+1,shift2 @@ -537,46 +615,46 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ endif enddo - if (state_following) then - - overlap = -1.d0 - do k=1,shift2 - do i=1,shift2 - overlap(k,i) = dabs(y(k,i)) - enddo - enddo - do k=1,N_st - cmax = -1.d0 - do i=1,N_st - if (overlap(i,k) > cmax) then - cmax = overlap(i,k) - order(k) = i - endif - enddo - do i=1,N_st_diag - overlap(order(k),i) = -1.d0 - enddo - enddo - overlap = y - do k=1,N_st - l = order(k) - if (k /= l) then - y(1:shift2,k) = overlap(1:shift2,l) - endif - enddo - do k=1,N_st - overlap(k,1) = lambda(k) - overlap(k,2) = s2(k) - enddo - do k=1,N_st - l = order(k) - if (k /= l) then - lambda(k) = overlap(l,1) - s2(k) = overlap(l,2) - endif - enddo - - endif +! if (state_following) then +! +! overlap = -1.d0 +! do k=1,shift2 +! do i=1,shift2 +! overlap(k,i) = dabs(y(k,i)) +! enddo +! enddo +! do k=1,N_st +! cmax = -1.d0 +! do i=1,N_st +! if (overlap(i,k) > cmax) then +! cmax = overlap(i,k) +! order(k) = i +! endif +! enddo +! do i=1,N_st_diag +! overlap(order(k),i) = -1.d0 +! enddo +! enddo +! overlap = y +! do k=1,N_st +! l = order(k) +! if (k /= l) then +! y(1:shift2,k) = overlap(1:shift2,l) +! endif +! enddo +! do k=1,N_st +! overlap(k,1) = lambda(k) +! overlap(k,2) = s2(k) +! enddo +! do k=1,N_st +! l = order(k) +! if (k /= l) then +! lambda(k) = overlap(l,1) +! s2(k) = overlap(l,2) +! endif +! enddo +! +! endif ! Express eigenvectors of h in the determinant basis diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 46ad8f78..59c8313a 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -123,6 +123,7 @@ END_PROVIDER endif enddo + if (N_states_diag > N_states_diag_save) then N_states_diag = N_states_diag_save TOUCH N_states_diag @@ -133,24 +134,101 @@ END_PROVIDER print *, 'Diagonalization of H using Lapack' allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) allocate (eigenvalues(N_det)) + if (s2_eig) then + double precision, parameter :: alpha = 0.1d0 allocate (H_prime(N_det,N_det) ) + H_prime(1:N_det,1:N_det) = H_matrix_all_dets(1:N_det,1:N_det) + & alpha * S2_matrix_all_dets(1:N_det,1:N_det) + do j=1,N_det H_prime(j,j) = H_prime(j,j) - alpha*expected_s2 enddo + call lapack_diag(eigenvalues,eigenvectors,H_prime,size(H_prime,1),N_det) call nullify_small_elements(N_det,N_det,eigenvectors,size(eigenvectors,1),1.d-12) + CI_electronic_energy(:) = 0.d0 i_state = 0 + allocate (s2_eigvalues(N_det)) allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& N_det,size(eigenvectors,1)) - if (only_expected_s2) then + + if (state_following) then + if (.not. only_expected_s2) then + print*,'' + print*,'!!! State following only available with only_expected_s2 = .True. !!!' + STOP + endif + if (N_det < N_states) then + print*,'' + print*,'!!! State following requires at least N_states determinants to be activated !!!' + STOP + endif + endif + + if (state_following .and. only_expected_s2) then + + integer :: state(N_states), idx,l + double precision :: omax + double precision, allocatable :: overlp(:) + logical :: used + logical, allocatable :: ok(:) + + allocate(overlp(N_det), ok(N_det)) + + i_state = 0 + state = 0 + do l = 1, N_states + + ! Overlap wrt each state + overlp = 0d0 + do k = 1, N_det + do i = 1, N_det + overlp(k) = overlp(k) + psi_coef(i,l) * eigenvectors(i,k) + enddo + enddo + + ! Idx of the state with the maximum overlap not already "used" + omax = 0d0 + idx = 0 + do k = 1, N_det + + ! Already used ? + used = .False. + do i = 1, N_states + if (state(i) == k) then + used = .True. + endif + enddo + + ! Maximum overlap + if (dabs(overlp(k)) > omax .and. .not. used) then + if (dabs(s2_eigvalues(k)-expected_s2) > 0.5d0) cycle + omax = dabs(overlp(k)) + idx = k + endif + enddo + + state(l) = idx + i_state +=1 + enddo + + deallocate(overlp, ok) + + do i = 1, i_state + index_good_state_array(i) = state(i) + good_state_array(i) = .True. + enddo + + else if (only_expected_s2) then + do j=1,N_det ! Select at least n_states states with S^2 values closed to "expected_s2" if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then @@ -158,17 +236,23 @@ END_PROVIDER index_good_state_array(i_state) = j good_state_array(j) = .True. endif + if(i_state.eq.N_states) then exit endif enddo + else + do j=1,N_det index_good_state_array(j) = j good_state_array(j) = .True. enddo + endif + if(i_state .ne.0)then + ! Fill the first "i_state" states that have a correct S^2 value do j = 1, i_state do i=1,N_det @@ -177,6 +261,7 @@ END_PROVIDER CI_electronic_energy(j) = eigenvalues(index_good_state_array(j)) CI_s2(j) = s2_eigvalues(index_good_state_array(j)) enddo + i_other_state = 0 do j = 1, N_det if(good_state_array(j))cycle @@ -201,6 +286,7 @@ END_PROVIDER print*,' as the CI_eigenvectors' print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' print*,'' + do j=1,min(N_states_diag,N_det) do i=1,N_det CI_eigenvectors(i,j) = eigenvectors(i,j) @@ -209,14 +295,18 @@ END_PROVIDER CI_s2(j) = s2_eigvalues(j) enddo endif + deallocate(index_good_state_array,good_state_array) deallocate(s2_eigvalues) + else + call lapack_diag(eigenvalues,eigenvectors, & H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) CI_electronic_energy(:) = 0.d0 call u_0_S2_u_0(CI_s2,eigenvectors,N_det,psi_det,N_int, & min(N_det,N_states_diag),size(eigenvectors,1)) + ! Select the "N_states_diag" states of lowest energy do j=1,min(N_det,N_states_diag) do i=1,N_det @@ -224,7 +314,9 @@ END_PROVIDER enddo CI_electronic_energy(j) = eigenvalues(j) enddo + endif + do k=1,N_states_diag CI_electronic_energy(k) = 0.d0 do j=1,N_det @@ -235,6 +327,7 @@ END_PROVIDER enddo enddo enddo + deallocate(eigenvectors,eigenvalues) endif diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index 6054b99c..269a441b 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -166,6 +166,10 @@ if(frozen_orb_scf)then integer :: iorb,jorb + ! active|core|active + !active | | 0 | + !core | 0 | | 0 + !active | | 0 | do i = 1, n_core_orb iorb = list_core(i) do j = 1, n_act_orb diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 26e390b7..20386b30 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -2041,3 +2041,22 @@ subroutine get_A_squared(A,n,A2) double precision, intent(out):: A2(n,n) call dgemm('N','N',n,n,n,1.d0,A,size(A,1),A,size(A,1),0.d0,A2,size(A2,1)) end + +subroutine get_AB_prod(A,n,m,B,l,AB) + implicit none + BEGIN_DOC +! AB = A B where A is n x m, B is m x l. Use the dgemm routine + END_DOC + double precision, intent(in) :: A(n,m),B(m,l) + integer, intent(in) :: n,m,l + double precision, intent(out):: AB(n,l) + if(size(A,2).ne.m.or.size(B,1).ne.m)then + print*,'error in get_AB_prod ! ' + print*,'matrices do not have the good dimension ' + print*,'size(A,2) = ',size(A,2) + print*,'size(B,1) = ',size(B,1) + print*,'m = ',m + stop + endif + call dgemm('N','N',n,l,m,1.d0,A,size(A,1),B,size(B,1),0.d0,AB,size(AB,1)) +end