From 9d3743e530f2b7d342778a32bc2ca89e36f97044 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 22 Mar 2024 14:56:39 +0100 Subject: [PATCH] added some providers and the first tutorial for plugins --- plugins/README.rst | 4 +- .../tuto_I/print_traces_on_e.irp.f | 24 ++++ .../tuto_plugins/tuto_I/print_two_e_h.irp.f | 32 +++++ .../tuto_plugins/tuto_I/traces_one_e.irp.f | 111 ++++++++++++++++++ plugins/tuto_plugins/tuto_I/tuto_I.rst | 65 +++++++--- src/ao_one_e_ints/ao_one_e_ints.irp.f | 10 ++ src/scf_utils/fock_matrix.irp.f | 4 + src/utils/linear_algebra.irp.f | 19 +++ 8 files changed, 250 insertions(+), 19 deletions(-) create mode 100644 plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f create mode 100644 plugins/tuto_plugins/tuto_I/print_two_e_h.irp.f create mode 100644 plugins/tuto_plugins/tuto_I/traces_one_e.irp.f diff --git a/plugins/README.rst b/plugins/README.rst index 7f3f3c75..7fc011a3 100644 --- a/plugins/README.rst +++ b/plugins/README.rst @@ -22,6 +22,8 @@ we will go through a series of examples that allow you to do the following thing IV) print out the one- and two-electron rdms, V) obtain the AOs and MOs on the DFT grid, together with the density, +How the tutorial will be done +----------------------------- This tuto is as follows: i) you READ THIS FILE UNTIL THE END in order to get the big picture and vocabulary, ii) you go to the directory qp2/plugins/tuto_plugins/ and you will find detailed tuto there for each of the 5 examples. @@ -32,7 +34,7 @@ The first thing to do is to be in the QPSH mode: you execute the qp2/bin/qpsh sc 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. -!!!! WARINING: The plugins are NECESSARILY located in qp2/plugins/ !!!! +!!!! WARNING: The plugins are NECESSARILY located in qp2/plugins/ !!!! Ex: If you want to create a plugin named "my_fancy_plugin" in the directory plugins/plugins_test/, this goes with the command qp plugins create -n my_fancy_plugin -r plugins_test/ diff --git a/plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f b/plugins/tuto_plugins/tuto_I/print_traces_on_e.irp.f new file mode 100644 index 00000000..2bf3b86b --- /dev/null +++ b/plugins/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/tuto_plugins/tuto_I/print_two_e_h.irp.f b/plugins/tuto_plugins/tuto_I/print_two_e_h.irp.f new file mode 100644 index 00000000..eaeb6c98 --- /dev/null +++ b/plugins/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/tuto_plugins/tuto_I/traces_one_e.irp.f b/plugins/tuto_plugins/tuto_I/traces_one_e.irp.f new file mode 100644 index 00000000..e71d49fc --- /dev/null +++ b/plugins/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/tuto_plugins/tuto_I/tuto_I.rst b/plugins/tuto_plugins/tuto_I/tuto_I.rst index 05db8635..fea07e3d 100644 --- a/plugins/tuto_plugins/tuto_I/tuto_I.rst +++ b/plugins/tuto_plugins/tuto_I/tuto_I.rst @@ -1,14 +1,15 @@ -====================================== -Tutorial for plugin I: One-e integrals -====================================== +===================================================================== +Tutorial for plugin I: One-e integrals (duration: 20 minutes at most) +===================================================================== -!!! Requirements: - a) you know how to create an EZFIO file and run calculations with QP +Requirements +------------ + a) You know how to create an EZFIO file and run calculations with QP (check the tuto: ``), - b) you have an EZFIO file in the sto-3g from the file H2.xyz in plugins/tuto_plugins, - and you have run an HF calculation giving an energy of -1.116759 a.u., - c) you made an qp set_file YOUR_EZFIO_FILE_FOR_H2 in order to be, - d) you have READ the ../README.rst file to HAVE THE VOCABULARY. + b) You have an EZFIO file with MOs created (with the '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 :) + c) You made an qp set_file YOUR_EZFIO_FILE_FOR_H2 in order to work on that ezfio folder, + d) You have READ the ../README.rst file to HAVE THE VOCABULARY. Our goals: ---------- @@ -22,14 +23,14 @@ I) Starting: creating the plugin We will go step-by-step through these plugins. The name of the plugin will be "plugin_I", and its location is in "tuto_plugins". -Therefore to create the plugin, we do +Therefore to create the plugin, we do: -$ qp plugins create -n plugin_I -r tuto_plugins -Then to an "ls" in qp2/plugins/tuto_plugins/ -and you will find a directory called "plugin_I". +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: - i) a "NEED" file that will eventually contain all the other modules/plugins needed by our "plugin_I" - ii) a "README.rst" file that you can AND SHOULD modify in order to document what is doing the plugin. + i) a "NEED" file that will eventually contain all the other modules/plugins needed by our "plugin_I" + ii) a "README.rst" file that you can AND SHOULD modify in order to document what is doing the plugin. iii) a "plugin_I.irp.f" file that is a program to be compiled and just printing "Hello world" II) Specifying the dependencies @@ -78,8 +79,8 @@ The variables that we need are ao_one_e_integrals mo_one_e_integrals You can check them with -irpman ao_one_e_integral -irpman mo_one_e_integral +irpman ao_one_e_integrals +irpman mo_one_e_integrals in order to get some information on where they are created, and many more information. We will modify the executable such that it prints out the integrals. @@ -87,7 +88,7 @@ We will modify the executable such that it prints out the integrals. IV) Printing out the one-electron integrals -------------------------------------------- We will create a program that will print out the one-electron integrals on the AO and MO basis. -You can then copy the file "print_one_e_h.irp.f" in your plugin. +You can then copy the file "print_one_e_h.irp.f" located in "plugins/tuto_plugins/tuto_I" in your plugin. In the file you will see that we simply browse the two arrays "ao_one_e_integrals" and "mo_one_e_integrals", which are global variables (providers) and we browse them until either "ao_num" or "mo_num" which are also providers representing the number of AOs or MOs. You can check these variables with irpman ! If you recompile using "ninja" as before, and another executable has been created "print_one_e_h". @@ -95,3 +96,31 @@ Then, you can run the program on the ezfio file by doing qp run print_one_e_h and will print out the data you need :) +By the way, as the 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 : ++) get_ao_two_e_integral for the AO basis ++) get_two_e_integral for the MO basis +check them with irpman ! +To print the two-electron integrals, you can copy the file "print_two_e_h.irp.f" in your plugin and recompile. +Then just run the program +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 "traces_one_e.irp.f" you will find the several new providers among which + a) trace_mo_one_e_ints : simply the sum of the diagonal matrix element of the one-electron integrals + b) trace_ao_one_e_ints : the corresponding trace on the AO basis : Sum(m,n) S^{-1}_{mn} h_{mn} + c) trace_ao_one_e_ints_from_mo : the trace on the AO basis with the integrals obtained first from the MO basis +As explained in these files, "trace_mo_one_e_ints" is equal to "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. +(You can check with "qp create_ezfio -h" for the option to create an EZFIO with cartesian basis functions) + +In the file "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 ! 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/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index 1942e542..c8fa8333 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