Commit b05e753c authored by Eirik F. Kjønstad's avatar Eirik F. Kjønstad
Browse files

Merge branch 'project_out_SC_merge' into 'development'

redundancies projection SC

See merge request !1557
parents 20674985 da87ebe9
Loading
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -258,6 +258,7 @@ some roots were converged and others below the linear dependence threshold. eT-p
- Optimized ``get_eri_1der_libcint`` to remove redundancy in derivative integrals. eT-program/eT!1511
- Optimized low memory CC2 Jacobian by replacing c1 transformations of the integrals, by c1 transformations of the Cholesky vectors. eT-program/eT!1518, eT-program/eT!1519, eT-program/eT!1520
- Optimized cube file generator. eT-program/eT!1543
- SC-QED-HF convergence is more robust by projecting out eta redundancies. eT-program/eT!!1557

### Documentation
- Updates to README. eT-program/eT!1247, eT-program/eT!1487
+2 −0
Original line number Diff line number Diff line
@@ -184,6 +184,8 @@ add_eT_runtest(restart_qed_ccsd_es_energy "eT;short;qed;ccsd;g
add_eT_runtest(sc_qed_hf_energy                             "eT;short;sc-qed;hf")
add_eT_runtest(restart_sc_qed_hf                            "eT;short;restart;sc-qed;hf")
#
add_eT_runtest(sc_qed_hf_intensivity_symmetry               "eT;short;sc-qed;hf")
#
add_eT_runtest(mp2_energy                                   "eT;short;hf;mp2")
add_eT_runtest(mp2_energy_qmmmnopol                         "eT;short;hf;mp2;qmmm;mmnopol")
add_eT_runtest(mp2_energy_qmfq                              "eT;short;hf;mp2;qmmm;fq")
+151 −15
Original line number Diff line number Diff line
@@ -1145,25 +1145,45 @@ contains
      !!
      !! Constructs the eta-eta Hessian matrix (E^(2)):
      !! E^(2)_rs = 1/omega^2 <R|[E_ss(b-b^+),[E_rr(b-b^+),H_SC]]|R>.
      !! Multiply the invert of it ((E^(2))^-1) to minus the gradient.
      !! The result is the second order step for the eta parameters.
      !!
      !! eta_step = - ((E^(2))^-1)*eta_grad
      !! The eigenvalue problem is projected in the space orthogonal
      !! to the one spanned by the eta redundancies:
      !! E^(2) * eta_step = - eta_grad -->
      !! --> P*E^(2) * P*eta_step = - P*eta_grad.
      !!
      !! The projector P is obtaiened after the diagonalization of
      !! overlap matrix S_pq = <R|E_pp E_qq|R>.
      !!
      !! The projected Hessian P*E^(2) is level-shifed so to ensure
      !! invertibility in case of singular values due to redundancies.
      !!
      !! The projected step P*eta_step is obtained solving:
      !! P*eta_step = - (P*E^(2) + level_shift*I)^-1 * P*eta_grad.
      !!
      use array_utilities, only: invert
      use array_initialization, only: identity_array
      use parameters

      implicit none

      class(sc_qed_hf) :: wf

      real(dp), dimension(:,:), allocatable :: hessian, D, h, hessian_inverse
      real(dp), dimension(:,:), allocatable :: D, h, work, overlap, P, identity

      real(dp), dimension(:,:), allocatable :: hessian, projected_hessian, hessian_inverse

      real(dp), dimension(wf%n_mo), intent(out) :: step

      real(dp) :: gaussian_factor, common_addend
      real(dp), dimension(:), allocatable :: eigenvalues, projected_gradient

      real(dp) :: gaussian_factor, common_addend, ddot

      integer :: r, s, q, t
      real(dp), parameter ::  threshold = 1.0d-7
      real(dp), parameter ::  level_shift = 1.0d-6

      logical :: regularize_Hessian = .False.

      integer :: r, s, q, t, info

      call mem%alloc(D, wf%n_mo, wf%n_mo)
      call wf%construct_density_dipole_basis(D)
@@ -1174,6 +1194,8 @@ contains

      call mem%alloc(hessian, wf%n_mo, wf%n_mo)

      ! Calculation of the Hessian.

      !$omp parallel do private(r, s, q, t, gaussian_factor, common_addend)
      do r = 1, wf%n_mo

@@ -1218,12 +1240,126 @@ contains
      !$omp end parallel do

      call mem%dealloc(h)

      call mem%alloc(overlap, wf%n_mo, wf%n_mo, set_zero=.True.)

      ! Overlap matrix S_pq = <R|E_pp E_qq|R> = d_ppqq + D_pq delta_pq.

      !$omp parallel do private(q)
      do q = 1, wf%n_mo
         overlap(q,q) = D(q,q)
      enddo
      !$omp end parallel do

      !$omp parallel do private(q, r)
      do q = 1, wf%n_mo
         do r = 1, wf%n_mo
            overlap(q,r) = overlap(q,r) + (D(q,q)*D(r,r)-half*D(q,r)*D(r,q))
         enddo
      enddo
      !$omp end parallel do

      call mem%dealloc(D)

      call mem%alloc(work, wf%n_mo, wf%n_mo)
      call mem%alloc(eigenvalues, wf%n_mo)

      call dsyev('V', 'U',     &
                  wf%n_mo,     &
                  overlap,     &
                  wf%n_mo,     &
                  eigenvalues, &
                  work,        &
                  wf%n_mo**2,  &
                  info)

      call mem%dealloc(work)

      call mem%alloc(P, wf%n_mo, wf%n_mo)

      ! Projector P = I - sum_i vT_i*v_i where vi are eigenvector
      ! associated to zero eigenvalues of the overlap matrix S.

      call identity_array(P, wf%n_mo)

      do q = 1, wf%n_mo

         if (abs(eigenvalues(q)) < threshold) then
            regularize_Hessian = .True.

            call   dger(wf%n_mo,       &
                        wf%n_mo,       &
                        -one,          &
                        overlap(:,q),  &
                        1,             &
                        overlap(:,q),  &
                        1,             &
                        P,             &
                        wf%n_mo)
         end if

      end do

      call mem%dealloc(eigenvalues)
      call mem%dealloc(overlap)

      call mem%alloc(projected_gradient, wf%n_mo)
      call mem%alloc(projected_hessian, wf%n_mo, wf%n_mo)

      ! Projecting the gradient eta_grad --> P*eta_grad.

      call dgemm('N','N',                &
                  wf%n_mo,               &
                  1,                     &
                  wf%n_mo,               &
                  one,                   &
                  P,                     &
                  wf%n_mo,               &
                  wf%eta_gradient,       &
                  wf%n_mo,               &
                  zero,                  &
                  projected_gradient,     &
                  wf%n_mo)

      ! Projecting the Hessian E^(2) --> P*E^(2).

      call dgemm('N','N',                     &
                  wf%n_mo,                    &
                  wf%n_mo,                    &
                  wf%n_mo,                    &
                  one,                        &
                  P,                          &
                  wf%n_mo,                    &
                  hessian,                    &
                  wf%n_mo,                    &
                  zero,                       &
                  projected_hessian,               &
                  wf%n_mo)

      call dcopy(wf%n_mo, projected_gradient, 1, wf%eta_gradient, 1)
      call dcopy(wf%n_mo**2, projected_hessian, 1, hessian, 1)

      call mem%dealloc(projected_gradient)
      call mem%dealloc(projected_hessian)
      call mem%dealloc(P)

      ! If Hessian is singular, level-shift is done.

      if (regularize_Hessian) then
         call mem%alloc(identity, wf%n_mo, wf%n_mo)
         call identity_array(identity, wf%n_mo)

         hessian = hessian + identity * level_shift

         call mem%dealloc(identity)
      endif

      call mem%alloc(hessian_inverse, wf%n_mo, wf%n_mo)

      call invert(hessian_inverse, hessian, wf%n_mo)

      ! Calculating the projected step P*eta_step.

      call dgemm('N','N',                &
                  wf%n_mo,               &
                  1,                     &
+435 −0

File added.

Preview size limit exceeded, changes collapsed.

+249 −0
Original line number Diff line number Diff line



                     eT 2.0 - an electronic structure program

  ------------------------------------------------------------------------
   Author list in alphabetical order:
  ------------------------------------------------------------------------
   J. H. Andersen, S. Angelico, A. Balbi, A. Barlini, A. Bianchi, C. 
   Cappelli, M. Castagnola, S. Coriani, S. D. Folkestad, Y. El Moutaoukal, 
   T. Giovannini, L. Goletto, T. S. Haugland, A. Hutcheson, I-M. Høyvik, 
   E. F. Kjønstad, H. Koch, M. T. Lexander, G. Marrazzini, T. Moitra, 
   R. H. Myhre, Y. Os, A. C. Paul, R. Paul, J. Pedersen, M. Rinaldi, 
   R. R. Riso, S. Roet, E. Ronca, F. Rossi, B. S. Sannes, M. Scavino, 
   A. K. Schnack-Petersen, A. S. Skeidsvoll, J. H. M. Trabski, Å. H. 
   Tveten
  ------------------------------------------------------------------------
   J. Chem. Phys. 152, 184103 (2020); https://doi.org/10.1063/5.0004713


   This is eT 2.0.0 A (development)
  -------------------------------------------------------------
  Configuration date:  2025-04-14 11:11:32 UTC +02:00
  Git branch:          project_out_SC_merge
  Git hash:            9742e869970d1b82624a1bb5de4dc3bb3f38aa6d
  Fortran compiler:    GNU 11.4.0
  C compiler:          GNU 11.4.0
  C++ compiler:        GNU 11.4.0
  LAPACK type:         SYSTEM_NATIVE
  BLAS type:           SYSTEM_NATIVE
  OEI type:            libcint
  ERI type:            libcint
  64-bit integers:     OFF
  OpenMP:              ON
  PCM:                 OFF
  Forced batching:     OFF
  Runtime checks:      OFF
  Initializing to NAN: OFF
  -------------------------------------------------------------

  Calculation start:2025-04-14 11:27:52 UTC +02:00

  Running on 12 OMP threads


  :: Input file
  =============

     Note: geometry section is excluded from this print.

     - system
        charge: 0

     - do
        ground state

     - memory
        available: 8

     - solver cholesky
     threshold: 1.0d-12

     - solver scf
        algorithm:          scf-diis
        energy threshold:   1.0d-10
        residual threshold: 1.0d-10
        max iterations: 1000
        write orbitals

     - method
      sc-qed-hf

     - boson
        modes:        1
        frequency:    {0.5}
        polarization: {0.270, -0.892, 0.363}
        coupling:     {0.1}

  Memory available for calculation: 8.000000 GB


  :: SC-QED-RHF wavefunction
  ==========================

     ================================================================================
                                    Geometry (angstrom)
     ================================================================================
     Atom           X                  Y                  Z         #internal  #input
     ================================================================================
     Basis: sto-3g
     O        1.151306358200     0.596196864600    -0.528619881900          1       1
     C        2.335856168400    -0.050626472200    -0.155984189800          2       2
     C        0.000000000000     0.000000000000     0.000000000000          3       3
     H        2.483504302300    -0.053830165400     0.945131400800          4       4
     H        3.176554004200     0.488613979100    -0.617545858200          5       5
     H        2.364079626500    -1.106062129300    -0.501707806000          6       6
     H        0.000000000000     0.000000000000     1.110975264200          7       7
     H       -0.119418084700    -1.052237046200    -0.335861183600          8       8
     H       -0.870532905900     0.576327038500    -0.347292217100          9       9
     O     1001.151306358200     0.596196864600    -0.528619881900         10      10
     C     1002.335856168400    -0.050626472200    -0.155984189800         11      11
     C     1000.000000000000     0.000000000000     0.000000000000         12      12
     H     1002.483504302300    -0.053830165400     0.945131400800         13      13
     H     1003.176554004200     0.488613979100    -0.617545858200         14      14
     H     1002.364079626500    -1.106062129300    -0.501707806000         15      15
     H     1000.000000000000     0.000000000000     1.110975264200         16      16
     H      999.880581915300    -1.052237046200    -0.335861183600         17      17
     H      999.129467094100     0.576327038500    -0.347292217100         18      18
     ================================================================================

     ================================================================================
                                    Geometry (a.u.)
     ================================================================================
     Atom           X                  Y                  Z         #internal  #input
     ================================================================================
     Basis: sto-3g
     O        2.175653702468     1.126648790418    -0.998946800791          1       1
     C        4.414128424652    -0.095670167111    -0.294767398484          2       2
     C        0.000000000000     0.000000000000     0.000000000000          3       3
     H        4.693142960526    -0.101724269846     1.786039499239          4       4
     H        6.002817087828     0.923346601133    -1.166992541357          5       5
     H        4.467463030749    -2.090154501130    -0.948090347896          6       6
     H        0.000000000000     0.000000000000     2.099438980504          7       7
     H       -0.225667474403    -1.988439835439    -0.634685652876          8       8
     H       -1.645068774573     1.089100260947    -0.656287175512          9       9
     O     1891.901778267530     1.126648790418    -0.998946800791         10      10
     C     1894.140252989714    -0.095670167111    -0.294767398484         11      11
     C     1889.726124565062     0.000000000000     0.000000000000         12      12
     H     1894.419267525588    -0.101724269846     1.786039499239         13      13
     H     1895.728941652890     0.923346601133    -1.166992541357         14      14
     H     1894.193587595811    -2.090154501130    -0.948090347896         15      15
     H     1889.726124565062     0.000000000000     2.099438980504         16      16
     H     1889.500457090659    -1.988439835439    -0.634685652876         17      17
     H     1888.081055790489     1.089100260947    -0.656287175512         18      18
     ================================================================================

  - Cholesky decomposition of AO overlap to get linearly independent AOs:

     Linear dependence threshold:             0.10E-05
     Number of atomic orbitals:               42
     Number of orthonormal atomic orbitals:   42

  - Molecular orbital details:

     Number of occupied orbitals:        26
     Number of virtual orbitals:         16
     Number of molecular orbitals:       42


  Generating initial SAD density
  ==============================


  Determining reference state
  ===========================

  - Setting initial AO density to sad

     Energy of initial guess:              -302.647732621831
     Number of electrons in guess:           52.000000000000

  - Screening and integral thresholds:

     Coulomb screening threshold:      0.1000E-15
     Exchange screening threshold:     0.1000E-13
     Cumulative Fock threshold:        0.0000E+00
     ERI precision:                    0.1000E-15
     One-electron integral precision:  0.1000E-20

  Self-consistent field solver
  ----------------------------

  A Self-consistent field solver that solves an eigenvalue equations: 
  F(C) C = C e.

  - SCF solver settings:

     Maximum iterations:                  1000
     Acceleration type:                   diis

  - Convergence thresholds:
     Residual threshold:            0.1000E-09
     Energy threshold:              0.1000E-09

  - DIIS tool settings:

     DIIS dimension:   8
     Storage (solver scf_errors): memory
     Storage (solver scf_parameters): memory

  Iteration       Energy (a.u.)      Max(grad.)    Delta E (a.u.)
  ---------------------------------------------------------------
     1          -304.058899719938     0.2273E+00     0.3041E+03
     2          -304.214475785147     0.1726E+00     0.1556E+00
     3          -304.221639876684     0.3194E-01     0.7164E-02
     4          -304.222964060590     0.3542E-02     0.1324E-02
     5          -304.223123421708     0.4696E-02     0.1594E-03
     6          -304.223137388731     0.2865E-03     0.1397E-04
     7          -304.223137571511     0.1805E-04     0.1828E-06
     8          -304.223137577751     0.1779E-04     0.6241E-08
     9          -304.223137578153     0.1038E-05     0.4020E-09
    10          -304.223137578175     0.6553E-06     0.2115E-10
    11          -304.223137578175     0.4674E-07     0.7390E-12
    12          -304.223137578175     0.5930E-08     0.2274E-12
    13          -304.223137578176     0.1895E-08     0.5684E-12
    14          -304.223137578174     0.5708E-09     0.1137E-11
    15          -304.223137578176     0.2121E-09     0.1137E-11
    16          -304.223137578175     0.1010E-10     0.7390E-12
  ---------------------------------------------------------------
  Convergence criterion met in 16 iterations!

  - Summary of SC-QED-RHF wavefunction energetics (a.u.):

     HOMO-LUMO gap:                    0.909366882695
     Nuclear repulsion energy:       168.084142730147
     Electronic energy:             -472.307280308322
     Total energy:                  -304.223137578175

  - Photon parameters and properties

     Optimize photons displacement?: True
     Complete basis approximation?:  False

     Mode 1:
     --------
      Frequency:              0.500000000000
      Coherent state z:       0.000000000000

      Polarization:           0.27000000  -0.89200000   0.36300000

      Coupling
         sqrt(1/eps V):       0.100000000000
         Bilinear:            0.050000000000
         Quadratic:           0.005000000000

  Peak memory usage during the execution of eT: 45.796064 MB

  Total wall time in eT (s): 6.768
  Total cpu time in eT (s):  64.006

  Calculation end:2025-04-14 11:27:58 UTC +02:00

  - Implementation references:

     eT: https://doi.org/10.1063/5.0004713
     SC-QED-HF: https://doi.org/https://doi.org/10.1038/s41467-022-29003-2
     Cholesky decomposition of ERIs: https://doi.org/10.1063/1.5083802

  eT terminated successfully!
Loading