PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 12:23:43 +01:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
Anthony Scemama 2024-03-27 14:18:55 +01:00
commit a4db5a87e0
20 changed files with 800 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
bin/qp_set_frozen_core
View File

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

1
docs/source/appendix/contributors.rst
View File

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

4
docs/source/index.rst
View File

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

1
docs/source/programmers_guide/plugins_tuto_I.rst Normal file
View File

@ -0,0 +1 @@
.. include:: ../../../plugins/local/tuto_plugins/tuto_I/tuto_I.rst

1
docs/source/programmers_guide/plugins_tuto_intro.rst Normal file
View File

@ -0,0 +1 @@
.. include:: ../../../plugins/README.rst

2
external/irpf90 vendored

@ -1 +1 @@
Subproject commit ba1a2837aa61cb8f9892860cec544d7c6659badd
Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6

131
plugins/README.rst Normal file
View File

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

2
plugins/local/cipsi_tc_bi_ortho/selection.irp.f
View File

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

6
plugins/local/tuto_plugins/H2.xyz Normal file
View File

@ -0,0 +1,6 @@
2
H2, equilibrium geometry
H 0.0 0.0 0.
H 0.0 0.0 0.74

4
plugins/local/tuto_plugins/n2.xyz Normal file
View File

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

20
plugins/local/tuto_plugins/tuto_I/print_one_e_h.irp.f Normal file
View File

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

24
plugins/local/tuto_plugins/tuto_I/print_traces_on_e.irp.f Normal file
View File

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

32
plugins/local/tuto_plugins/tuto_I/print_two_e_h.irp.f Normal file
View File

@ -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 : <i j|k l>'
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 : <i j|k l>'
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

111
plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f Normal file
View File

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

218
plugins/local/tuto_plugins/tuto_I/tuto_I.rst Normal file
View File

@ -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: `<https://quantumpackage.github.io/qp2/post/hartree-fock/>`_),
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 <<EOF > 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.

10
src/ao_one_e_ints/ao_one_e_ints.irp.f
View File

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

158
src/davidson/diagonalization_hs2_dressed.irp.f
View File

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

95
src/davidson/diagonalize_ci.irp.f
View File

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

4
src/scf_utils/fock_matrix.irp.f
View File

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

19
src/utils/linear_algebra.irp.f
View File

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