LCOV - code coverage report
Current view: top level - src - atom_electronic_structure.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 321 376 85.4 %
Date: 2024-03-29 07:50:05 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : MODULE atom_electronic_structure
       9             :    USE atom_optimization,               ONLY: atom_history_init,&
      10             :                                               atom_history_release,&
      11             :                                               atom_history_type,&
      12             :                                               atom_history_update,&
      13             :                                               atom_opt_fmat
      14             :    USE atom_output,                     ONLY: atom_print_energies,&
      15             :                                               atom_print_iteration,&
      16             :                                               atom_print_state,&
      17             :                                               atom_print_zmp_iteration
      18             :    USE atom_types,                      ONLY: &
      19             :         atom_type, create_opgrid, create_opmat, ecp_pseudo, gth_pseudo, lmat, no_pseudo, &
      20             :         opgrid_type, opmat_type, release_opgrid, release_opmat, set_atom, setup_hf_section, &
      21             :         sgp_pseudo, upf_pseudo
      22             :    USE atom_utils,                      ONLY: &
      23             :         atom_denmat, atom_density, atom_read_external_density, atom_read_external_vxc, &
      24             :         atom_read_zmp_restart, atom_solve, atom_trace, ceri_contract, coulomb_potential_analytic, &
      25             :         coulomb_potential_numeric, eeri_contract, err_matrix, exchange_numeric, &
      26             :         exchange_semi_analytic, numpot_matrix, slater_density, wigner_slater_functional
      27             :    USE atom_xc,                         ONLY: calculate_atom_ext_vxc,&
      28             :                                               calculate_atom_vxc_lda,&
      29             :                                               calculate_atom_vxc_lsd,&
      30             :                                               calculate_atom_zmp
      31             :    USE input_constants,                 ONLY: &
      32             :         do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
      33             :         do_numeric, do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_semi_analytic, &
      34             :         do_uhf_atom, do_uks_atom, do_zoramp_atom
      35             :    USE input_section_types,             ONLY: section_vals_get,&
      36             :                                               section_vals_get_subs_vals,&
      37             :                                               section_vals_type
      38             :    USE kinds,                           ONLY: dp
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             :    PRIVATE
      43             :    PUBLIC  :: calculate_atom
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_electronic_structure'
      46             : 
      47             : CONTAINS
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief General routine to perform electronic structure atomic calculations.
      51             : !> \param atom       information about the atomic kind
      52             : !> \param iw         output file unit
      53             : !> \param noguess    skip initial guess
      54             : !> \param converged  whether SCF iterations have been converged
      55             : !> \par History
      56             : !>    * 06.2017 disable XC input options [Juerg Hutter]
      57             : !>    * 11.2009 created from the subroutine calculate_atom() [Juerg Hutter]
      58             : ! **************************************************************************************************
      59       11116 :    SUBROUTINE calculate_atom(atom, iw, noguess, converged)
      60             :       TYPE(atom_type), POINTER                           :: atom
      61             :       INTEGER, INTENT(IN)                                :: iw
      62             :       LOGICAL, INTENT(IN), OPTIONAL                      :: noguess
      63             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: converged
      64             : 
      65             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calculate_atom'
      66             : 
      67             :       INTEGER                                            :: handle, method
      68             :       LOGICAL                                            :: explicit
      69             :       TYPE(section_vals_type), POINTER                   :: sub_section
      70             : 
      71       11116 :       CALL timeset(routineN, handle)
      72             : 
      73             :       ! test for not supported methods and functionals
      74       11116 :       IF (ASSOCIATED(atom%xc_section)) THEN
      75             :          !
      76        2930 :          sub_section => section_vals_get_subs_vals(atom%xc_section, "ADIABATIC_RESCALING")
      77        2930 :          CALL section_vals_get(sub_section, explicit=explicit)
      78        2930 :          IF (explicit) CALL cp_abort(__LOCATION__, "ADIABATIC_RESCALING not supported in ATOM code")
      79             :          !
      80        2930 :          sub_section => section_vals_get_subs_vals(atom%xc_section, "VDW_POTENTIAL")
      81        2930 :          CALL section_vals_get(sub_section, explicit=explicit)
      82        2930 :          IF (explicit) CALL cp_abort(__LOCATION__, "VDW_POTENTIAL not supported in ATOM code")
      83             :          !
      84        2930 :          sub_section => section_vals_get_subs_vals(atom%xc_section, "XC_POTENTIAL")
      85        2930 :          CALL section_vals_get(sub_section, explicit=explicit)
      86        2930 :          IF (explicit) CALL cp_abort(__LOCATION__, "XC_POTENTIAL not supported in ATOM code")
      87             :          !
      88        2930 :          sub_section => section_vals_get_subs_vals(atom%xc_section, "WF_CORRELATION")
      89        2930 :          CALL section_vals_get(sub_section, explicit=explicit)
      90        2930 :          IF (explicit) CALL cp_abort(__LOCATION__, "WF_CORRELATION methods not supported in ATOM code")
      91             :          !
      92             :       END IF
      93             : 
      94       11116 :       method = atom%method_type
      95       10999 :       SELECT CASE (method)
      96             :       CASE (do_rks_atom, do_rhf_atom)
      97       10999 :          CALL calculate_atom_restricted(atom, iw, noguess, converged)
      98             :       CASE (do_uks_atom, do_uhf_atom)
      99         117 :          CALL calculate_atom_unrestricted(atom, iw, noguess, converged)
     100             :       CASE (do_rohf_atom)
     101           0 :          CPABORT("")
     102             :       CASE DEFAULT
     103       11116 :          CPABORT("")
     104             :       END SELECT
     105             : 
     106       11116 :       CALL timestop(handle)
     107             : 
     108       11116 :    END SUBROUTINE calculate_atom
     109             : 
     110             : ! **************************************************************************************************
     111             : !> \brief Perform restricted (closed shell) electronic structure atomic calculations.
     112             : !> \param atom       information about the atomic kind
     113             : !> \param iw         output file unit
     114             : !> \param noguess    skip initial guess
     115             : !> \param converged  whether SCF iterations have been converged
     116             : !> \par History
     117             : !>    * 02.2016 support UPF files and ECP potentials [Juerg Hutter]
     118             : !>    * 09.2015 Polarized Atomic Orbitals (PAO) method [Ole Schuett]
     119             : !>    * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
     120             : !>    * 07.2013 scaled ZORA with model potential [Juerg Hutter]
     121             : !>    * 11.2009 split into three subroutines calculate_atom(), calculate_atom_restricted(),
     122             : !>              and calculate_atom_unrestricted() [Juerg Hutter]
     123             : !>    * 12.2008 refactored and renamed to calculate_atom() [Juerg Hutter]
     124             : !>    * 08.2008 created as atom_electronic_structure() [Juerg Hutter]
     125             : ! **************************************************************************************************
     126       10999 :    SUBROUTINE calculate_atom_restricted(atom, iw, noguess, converged)
     127             :       TYPE(atom_type), POINTER                           :: atom
     128             :       INTEGER, INTENT(IN)                                :: iw
     129             :       LOGICAL, INTENT(IN), OPTIONAL                      :: noguess
     130             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: converged
     131             : 
     132             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atom_restricted'
     133             : 
     134             :       INTEGER                                            :: handle, i, iter, l, max_iter, method, &
     135             :                                                             reltyp
     136             :       LOGICAL                                            :: do_hfx, doguess, is_converged, need_vxc, &
     137             :                                                             need_x, need_xc, need_zmp
     138             :       REAL(KIND=dp)                                      :: deps, eps_scf, hf_frac
     139       10999 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ext_density, ext_vxc, tmp_dens
     140             :       TYPE(atom_history_type)                            :: history
     141             :       TYPE(opgrid_type), POINTER                         :: cpot, density
     142             :       TYPE(opmat_type), POINTER                          :: fmat, hcore, jmat, kmat, xcmat
     143             :       TYPE(section_vals_type), POINTER                   :: xc_section
     144             : 
     145       10999 :       CALL timeset(routineN, handle)
     146             : 
     147       10999 :       IF (PRESENT(noguess)) THEN
     148         329 :          doguess = .NOT. noguess
     149             :       ELSE
     150       10670 :          doguess = .TRUE.
     151             :       END IF
     152             : 
     153       10999 :       CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
     154             : 
     155       10999 :       method = atom%method_type
     156       10999 :       max_iter = atom%optimization%max_iter
     157       10999 :       eps_scf = atom%optimization%eps_scf
     158             : 
     159           0 :       SELECT CASE (method)
     160             :       CASE DEFAULT
     161           0 :          CPABORT("")
     162             :       CASE (do_rks_atom)
     163        9770 :          need_x = do_hfx
     164        9770 :          need_xc = .TRUE.
     165             :       CASE (do_uks_atom)
     166           0 :          CPABORT("")
     167             :       CASE (do_rhf_atom)
     168        1229 :          need_x = .TRUE.
     169        1229 :          need_xc = .FALSE.
     170        1229 :          hf_frac = 1._dp
     171             :       CASE (do_uhf_atom)
     172           0 :          CPABORT("")
     173             :       CASE (do_rohf_atom)
     174           0 :          need_x = .TRUE.
     175           0 :          need_xc = .FALSE.
     176           0 :          hf_frac = 1._dp
     177       10999 :          CPABORT("")
     178             :       END SELECT
     179             : 
     180             :       ! ZMP starting to read external density for the zmp calculation
     181       10999 :       need_zmp = atom%do_zmp
     182       10999 :       IF (need_zmp) THEN
     183           0 :          need_x = .FALSE.
     184           0 :          need_xc = .FALSE.
     185           0 :          ALLOCATE (ext_density(atom%basis%grid%nr))
     186           0 :          ext_density = 0._dp
     187           0 :          CALL atom_read_external_density(ext_density, atom, iw)
     188             :       END IF
     189             : 
     190             :       ! ZMP starting to read external vxc for the zmp calculation
     191       10999 :       need_vxc = atom%read_vxc
     192             : 
     193       10999 :       IF (need_vxc) THEN
     194           0 :          need_x = .FALSE.
     195           0 :          need_xc = .FALSE.
     196           0 :          need_zmp = .FALSE.
     197           0 :          ALLOCATE (ext_vxc(atom%basis%grid%nr))
     198           0 :          ext_vxc = 0._dp
     199           0 :          CALL atom_read_external_vxc(ext_vxc, atom, iw)
     200             :       END IF
     201             : 
     202             :       ! check for relativistic method
     203       10999 :       reltyp = atom%relativistic
     204             : 
     205       10999 :       IF (iw > 0) CALL atom_print_state(atom%state, iw)
     206             : 
     207       10999 :       NULLIFY (hcore)
     208       10999 :       CALL create_opmat(hcore, atom%basis%nbas)
     209             :       ! Pseudopotentials
     210       10999 :       SELECT CASE (atom%potential%ppot_type)
     211             :       CASE DEFAULT
     212           0 :          CPABORT("")
     213             :       CASE (NO_PSEUDO)
     214        8657 :          SELECT CASE (reltyp)
     215             :          CASE DEFAULT
     216           0 :             CPABORT("")
     217             :          CASE (do_nonrel_atom)
     218     2561772 :             hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
     219             :          CASE (do_zoramp_atom, do_sczoramp_atom)
     220      242752 :             hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
     221             :          CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
     222      641486 :             hcore%op = atom%integrals%hdkh
     223             :          END SELECT
     224             :       CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
     225     2636395 :          hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
     226             :       END SELECT
     227             :       ! add confinement potential (not included in relativistic transformations)
     228       10999 :       IF (atom%potential%confinement) THEN
     229     2914143 :          hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
     230             :       END IF
     231             : 
     232       10999 :       NULLIFY (fmat, jmat, kmat, xcmat)
     233       10999 :       CALL create_opmat(fmat, atom%basis%nbas)
     234       10999 :       CALL create_opmat(jmat, atom%basis%nbas)
     235       10999 :       CALL create_opmat(kmat, atom%basis%nbas)
     236       10999 :       CALL create_opmat(xcmat, atom%basis%nbas)
     237             : 
     238       10999 :       NULLIFY (density, cpot)
     239       10999 :       CALL create_opgrid(density, atom%basis%grid)
     240       10999 :       CALL create_opgrid(cpot, atom%basis%grid)
     241             : 
     242             :       ! ZMP reading the file to restart
     243       10999 :       IF (atom%doread) CALL atom_read_zmp_restart(atom, doguess, iw)
     244             : 
     245       10999 :       IF (doguess) THEN
     246             :          ! initial guess
     247       32046 :          ALLOCATE (tmp_dens(SIZE(density%op)))
     248       10682 :          CALL slater_density(density%op, tmp_dens, atom%z, atom%state, atom%basis%grid)
     249     4286252 :          density%op = density%op + tmp_dens
     250       10682 :          DEALLOCATE (tmp_dens)
     251       10682 :          CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
     252       10682 :          CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     253       10682 :          CALL wigner_slater_functional(density%op, cpot%op)
     254       10682 :          CALL numpot_matrix(xcmat%op, cpot%op, atom%basis, 0)
     255     5049938 :          fmat%op = hcore%op + jmat%op + xcmat%op
     256             :          CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
     257       10682 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     258             :       END IF
     259             :       CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
     260       10999 :                        atom%state%maxl_occ, atom%state%maxn_occ)
     261             : 
     262             :       ! wavefunction history
     263       10999 :       NULLIFY (history%dmat, history%hmat)
     264       10999 :       CALL atom_history_init(history, atom%optimization, fmat%op)
     265       10999 :       is_converged = .FALSE.
     266       10999 :       iter = 0
     267       65173 :       DO !SCF Loop
     268             : 
     269             :          ! Kinetic energy
     270       76172 :          atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
     271             : 
     272             :          ! Band energy
     273       76172 :          atom%energy%eband = 0._dp
     274      533204 :          DO l = 0, lmat
     275     1092716 :             DO i = 1, MIN(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
     276     1016544 :                atom%energy%eband = atom%energy%eband + atom%orbitals%ener(i, l)*atom%state%occupation(l, i)
     277             :             END DO
     278             :          END DO
     279             : 
     280             :          ! Pseudopotential energy
     281       76172 :          SELECT CASE (atom%potential%ppot_type)
     282             :          CASE DEFAULT
     283           0 :             CPABORT("")
     284             :          CASE (no_pseudo)
     285       26314 :             atom%energy%eploc = 0._dp
     286       26314 :             atom%energy%epnl = 0._dp
     287             :          CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
     288       49858 :             atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
     289       76172 :             atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
     290             :          END SELECT
     291       76172 :          atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
     292             : 
     293             :          ! Core energy
     294       76172 :          atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
     295             : 
     296             :          ! Confinement energy
     297       76172 :          IF (atom%potential%confinement) THEN
     298       49602 :             atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
     299             :          ELSE
     300       26570 :             atom%energy%econfinement = 0._dp
     301             :          END IF
     302             : 
     303             :          ! Hartree Term
     304    30439268 :          jmat%op = 0._dp
     305       76172 :          SELECT CASE (atom%coulomb_integral_type)
     306             :          CASE DEFAULT
     307           0 :             CPABORT("")
     308             :          CASE (do_analytic)
     309        1505 :             CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
     310             :          CASE (do_semi_analytic)
     311             :             CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
     312          38 :                                             atom%state%maxl_occ)
     313          38 :             CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     314             :          CASE (do_numeric)
     315       74629 :             CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
     316       74629 :             CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
     317      150801 :             CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     318             :          END SELECT
     319       76172 :          atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
     320             : 
     321             :          ! Exchange Term
     322       76172 :          IF (need_x) THEN
     323     5399476 :             kmat%op = 0._dp
     324       21016 :             SELECT CASE (atom%exchange_integral_type)
     325             :             CASE DEFAULT
     326           0 :                CPABORT("")
     327             :             CASE (do_analytic)
     328         535 :                CALL eeri_contract(kmat%op, atom%integrals%eeri, atom%orbitals%pmat, atom%integrals%n)
     329             :             CASE (do_semi_analytic)
     330         191 :                CALL exchange_semi_analytic(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
     331             :             CASE (do_numeric)
     332       21016 :                CALL exchange_numeric(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
     333             :             END SELECT
     334       21016 :             atom%energy%eexchange = hf_frac*0.5_dp*atom_trace(kmat%op, atom%orbitals%pmat)
     335     5399476 :             kmat%op = hf_frac*kmat%op
     336             :          ELSE
     337    25039792 :             kmat%op = 0._dp
     338       55156 :             atom%energy%eexchange = 0._dp
     339             :          END IF
     340             : 
     341             :          ! XC
     342       76172 :          IF (need_xc) THEN
     343    26500752 :             xcmat%op = 0._dp
     344       55728 :             CALL calculate_atom_vxc_lda(xcmat, atom, xc_section)
     345             :             ! ZMP added options for the zmp calculations, building external density and vxc potential
     346       20444 :          ELSEIF (need_zmp) THEN
     347           0 :             xcmat%op = 0._dp
     348           0 :             CALL calculate_atom_zmp(ext_density=ext_density, atom=atom, lprint=.FALSE., xcmat=xcmat)
     349       20444 :          ELSEIF (need_vxc) THEN
     350           0 :             xcmat%op = 0._dp
     351           0 :             CALL calculate_atom_ext_vxc(vxc=ext_vxc, atom=atom, lprint=.FALSE., xcmat=xcmat)
     352             :          ELSE
     353     3938516 :             xcmat%op = 0._dp
     354       20444 :             atom%energy%exc = 0._dp
     355             :          END IF
     356             : 
     357             :          ! Zero this contribution
     358       76172 :          atom%energy%elsd = 0._dp
     359             : 
     360             :          ! Total energy
     361       76172 :          atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
     362             : 
     363             :          ! Potential energy
     364       76172 :          atom%energy%epot = atom%energy%etot - atom%energy%ekin
     365             : 
     366             :          ! Total HF/KS matrix
     367    60802364 :          fmat%op = hcore%op + jmat%op + kmat%op + xcmat%op
     368             : 
     369             :          ! calculate error matrix
     370             :          CALL err_matrix(jmat%op, deps, fmat%op, atom%orbitals%pmat, atom%integrals%utrans, &
     371       76172 :                          atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
     372             : 
     373       76172 :          iter = iter + 1
     374             : 
     375       76172 :          IF (iw > 0) THEN
     376        9211 :             IF (need_zmp) THEN
     377           0 :                CALL atom_print_zmp_iteration(iter, deps, atom, iw)
     378             :             ELSE
     379        9211 :                CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
     380             :             END IF
     381             :          END IF
     382             : 
     383       76172 :          IF (deps < eps_scf) THEN
     384             :             is_converged = .TRUE.
     385             :             EXIT
     386             :          END IF
     387       65173 :          IF (iter >= max_iter) THEN
     388           0 :             IF (iw > 0) THEN
     389           0 :                WRITE (iw, "(A)") " No convergence within maximum number of iterations "
     390             :             END IF
     391             :             EXIT
     392             :          END IF
     393             : 
     394             :          ! update history container and extrapolate KS matrix
     395       65173 :          CALL atom_history_update(history, atom%orbitals%pmat, fmat%op, jmat%op, atom%energy%etot, deps)
     396       65173 :          CALL atom_opt_fmat(fmat%op, history, deps)
     397             : 
     398             :          ! Solve HF/KS equations
     399             :          CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
     400       65173 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     401             :          CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
     402       65173 :                           atom%state%maxl_occ, atom%state%maxn_occ)
     403             : 
     404             :       END DO !SCF Loop
     405             : 
     406       10999 :       IF (iw > 0) THEN
     407        2131 :          CALL atom_print_energies(atom, iw)
     408             :       END IF
     409             : 
     410       10999 :       CALL atom_history_release(history)
     411             : 
     412             :       ! conserve fmat within atom_type
     413       10999 :       CALL set_atom(atom, fmat=fmat)
     414       10999 :       CALL release_opmat(jmat)
     415       10999 :       CALL release_opmat(kmat)
     416       10999 :       CALL release_opmat(xcmat)
     417       10999 :       CALL release_opmat(hcore)
     418             : 
     419       10999 :       CALL release_opgrid(density)
     420       10999 :       CALL release_opgrid(cpot)
     421             : 
     422             :       ! ZMP deallocating ext_density ext_vxc
     423       10999 :       IF (need_zmp) DEALLOCATE (ext_density)
     424       10999 :       IF (need_vxc) DEALLOCATE (ext_vxc)
     425             : 
     426       10999 :       IF (PRESENT(converged)) THEN
     427         310 :          converged = is_converged
     428             :       END IF
     429             : 
     430       10999 :       CALL timestop(handle)
     431             : 
     432       21998 :    END SUBROUTINE calculate_atom_restricted
     433             : 
     434             : ! **************************************************************************************************
     435             : !> \brief Perform unrestricted (spin-polarised) electronic structure atomic calculations.
     436             : !> \param atom       information about the atomic kind
     437             : !> \param iw         output file unit
     438             : !> \param noguess    skip initial guess
     439             : !> \param converged  whether SCF iterations have been converged
     440             : !> \par History
     441             : !>    * identical to the subroutine calculate_atom_restricted()
     442             : ! **************************************************************************************************
     443         234 :    SUBROUTINE calculate_atom_unrestricted(atom, iw, noguess, converged)
     444             :       TYPE(atom_type), POINTER                           :: atom
     445             :       INTEGER, INTENT(IN)                                :: iw
     446             :       LOGICAL, INTENT(IN), OPTIONAL                      :: noguess
     447             :       LOGICAL, INTENT(OUT), OPTIONAL                     :: converged
     448             : 
     449             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atom_unrestricted'
     450             : 
     451             :       INTEGER                                            :: handle, i, iter, k, l, max_iter, method, &
     452             :                                                             reltyp
     453             :       LOGICAL                                            :: do_hfx, doguess, is_converged, lsdpot, &
     454             :                                                             need_x, need_xc
     455             :       REAL(KIND=dp)                                      :: deps, depsa, depsb, eps_scf, hf_frac, &
     456             :                                                             ne, nm
     457             :       TYPE(atom_history_type)                            :: historya, historyb
     458             :       TYPE(opgrid_type), POINTER                         :: cpot, density, rhoa, rhob
     459             :       TYPE(opmat_type), POINTER                          :: fmata, fmatb, hcore, hlsd, jmat, kmata, &
     460             :                                                             kmatb, xcmata, xcmatb
     461             :       TYPE(section_vals_type), POINTER                   :: xc_section
     462             : 
     463         117 :       CALL timeset(routineN, handle)
     464             : 
     465         117 :       IF (PRESENT(noguess)) THEN
     466          22 :          doguess = .NOT. noguess
     467             :       ELSE
     468             :          doguess = .TRUE.
     469             :       END IF
     470             : 
     471         117 :       CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
     472             : 
     473         117 :       method = atom%method_type
     474         117 :       max_iter = atom%optimization%max_iter
     475         117 :       eps_scf = atom%optimization%eps_scf
     476             : 
     477           0 :       SELECT CASE (method)
     478             :       CASE DEFAULT
     479           0 :          CPABORT("")
     480             :       CASE (do_rks_atom)
     481           0 :          CPABORT("")
     482             :       CASE (do_uks_atom)
     483         113 :          need_x = do_hfx
     484         113 :          need_xc = .TRUE.
     485             :       CASE (do_rhf_atom)
     486           0 :          CPABORT("")
     487             :       CASE (do_uhf_atom)
     488           4 :          need_x = .TRUE.
     489           4 :          need_xc = .FALSE.
     490           4 :          hf_frac = 1._dp
     491             :       CASE (do_rohf_atom)
     492           0 :          need_x = .TRUE.
     493           0 :          need_xc = .FALSE.
     494           0 :          hf_frac = 1._dp
     495         117 :          CPABORT("")
     496             :       END SELECT
     497             : 
     498             :       ! set alpha and beta occupations
     499        8307 :       IF (SUM(ABS(atom%state%occa) + ABS(atom%state%occb)) == 0.0_dp) THEN
     500         240 :          DO l = 0, 3
     501         192 :             nm = REAL((2*l + 1), KIND=dp)
     502         342 :             DO k = 1, 10
     503         294 :                ne = atom%state%occupation(l, k)
     504         294 :                IF (ne == 0._dp) THEN !empty shell
     505             :                   EXIT !assume there are no holes
     506         102 :                ELSEIF (ne == 2._dp*nm) THEN !closed shell
     507          48 :                   atom%state%occa(l, k) = nm
     508          48 :                   atom%state%occb(l, k) = nm
     509          54 :                ELSEIF (atom%state%multiplicity == -2) THEN !High spin case
     510          40 :                   atom%state%occa(l, k) = MIN(ne, nm)
     511          40 :                   atom%state%occb(l, k) = MAX(0._dp, ne - nm)
     512             :                ELSE
     513          14 :                   atom%state%occa(l, k) = 0.5_dp*(ne + atom%state%multiplicity - 1._dp)
     514          14 :                   atom%state%occb(l, k) = ne - atom%state%occa(l, k)
     515             :                END IF
     516             :             END DO
     517             :          END DO
     518             :       END IF
     519             :       ! check for relativistic method
     520         117 :       reltyp = atom%relativistic
     521             : 
     522         117 :       IF (iw > 0) CALL atom_print_state(atom%state, iw)
     523             : 
     524         117 :       NULLIFY (hcore, hlsd)
     525         117 :       CALL create_opmat(hcore, atom%basis%nbas)
     526         117 :       CALL create_opmat(hlsd, atom%basis%nbas)
     527      275739 :       hlsd%op = 0._dp
     528             :       ! Pseudopotentials
     529         117 :       lsdpot = .FALSE.
     530         117 :       SELECT CASE (atom%potential%ppot_type)
     531             :       CASE DEFAULT
     532           0 :          CPABORT("")
     533             :       CASE (no_pseudo)
     534          93 :          SELECT CASE (reltyp)
     535             :          CASE DEFAULT
     536           0 :             CPABORT("")
     537             :          CASE (do_nonrel_atom)
     538       86856 :             hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
     539             :          CASE (do_zoramp_atom, do_sczoramp_atom)
     540           0 :             hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
     541             :          CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
     542          24 :             hcore%op = atom%integrals%hdkh
     543             :          END SELECT
     544             :       CASE (gth_pseudo)
     545      464505 :          hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
     546          93 :          IF (atom%potential%gth_pot%lsdpot) THEN
     547           0 :             lsdpot = .TRUE.
     548           0 :             hlsd%op = atom%integrals%clsd
     549             :          END IF
     550             :       CASE (upf_pseudo, sgp_pseudo, ecp_pseudo)
     551         117 :          hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
     552             :       END SELECT
     553             :       ! add confinement potential (not included in relativistic transformations)
     554         117 :       IF (atom%potential%confinement) THEN
     555      496413 :          hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
     556             :       END IF
     557             : 
     558         117 :       NULLIFY (fmata, fmatb, jmat, kmata, kmatb, xcmata, xcmatb)
     559         117 :       CALL create_opmat(fmata, atom%basis%nbas)
     560         117 :       CALL create_opmat(fmatb, atom%basis%nbas)
     561         117 :       CALL create_opmat(jmat, atom%basis%nbas)
     562         117 :       CALL create_opmat(kmata, atom%basis%nbas)
     563         117 :       CALL create_opmat(kmatb, atom%basis%nbas)
     564         117 :       CALL create_opmat(xcmata, atom%basis%nbas)
     565         117 :       CALL create_opmat(xcmatb, atom%basis%nbas)
     566             : 
     567         117 :       NULLIFY (density, rhoa, rhob, cpot)
     568         117 :       CALL create_opgrid(density, atom%basis%grid)
     569         117 :       CALL create_opgrid(rhoa, atom%basis%grid)
     570         117 :       CALL create_opgrid(rhob, atom%basis%grid)
     571         117 :       CALL create_opgrid(cpot, atom%basis%grid)
     572             : 
     573         117 :       IF (doguess) THEN
     574             :          ! initial guess
     575          96 :          CALL slater_density(rhoa%op, rhob%op, atom%z, atom%state, atom%basis%grid)
     576       78496 :          density%op = rhoa%op + rhob%op
     577          96 :          CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
     578          96 :          CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     579             :          ! alpha spin
     580       78496 :          density%op = 2._dp*rhoa%op
     581          96 :          CALL wigner_slater_functional(density%op, cpot%op)
     582          96 :          CALL numpot_matrix(xcmata%op, cpot%op, atom%basis, 0)
     583      137808 :          fmata%op = hcore%op + hlsd%op + jmat%op + xcmata%op
     584             :          CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
     585          96 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     586             :          ! beta spin
     587       78496 :          density%op = 2._dp*rhob%op
     588          96 :          CALL wigner_slater_functional(density%op, cpot%op)
     589          96 :          CALL numpot_matrix(xcmatb%op, cpot%op, atom%basis, 0)
     590      137808 :          fmatb%op = hcore%op - hlsd%op + jmat%op + xcmatb%op
     591             :          CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
     592          96 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     593             :       END IF
     594             :       CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
     595         117 :                        atom%state%maxl_occ, atom%state%maxn_occ)
     596             :       CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
     597         117 :                        atom%state%maxl_occ, atom%state%maxn_occ)
     598      275739 :       atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
     599             : 
     600             :       ! wavefunction history
     601         117 :       NULLIFY (historya%dmat, historya%hmat)
     602         117 :       CALL atom_history_init(historya, atom%optimization, fmata%op)
     603         117 :       NULLIFY (historyb%dmat, historyb%hmat)
     604         117 :       CALL atom_history_init(historyb, atom%optimization, fmatb%op)
     605             : 
     606         117 :       is_converged = .FALSE.
     607         117 :       iter = 0
     608             :       DO !SCF Loop
     609             : 
     610             :          ! Kinetic energy
     611        1351 :          atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
     612             : 
     613             :          ! Band energy
     614        1351 :          atom%energy%eband = 0._dp
     615        6755 :          DO l = 0, 3
     616       17231 :             DO i = 1, MIN(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
     617       10476 :                atom%energy%eband = atom%energy%eband + atom%orbitals%enera(i, l)*atom%state%occa(l, i)
     618       15880 :                atom%energy%eband = atom%energy%eband + atom%orbitals%enerb(i, l)*atom%state%occb(l, i)
     619             :             END DO
     620             :          END DO
     621             : 
     622             :          ! Pseudopotential energy
     623        1351 :          SELECT CASE (atom%potential%ppot_type)
     624             :          CASE DEFAULT
     625           0 :             CPABORT("")
     626             :          CASE (no_pseudo)
     627         148 :             atom%energy%eploc = 0._dp
     628         148 :             atom%energy%epnl = 0._dp
     629             :          CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
     630        1203 :             atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
     631        1351 :             atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
     632             :          END SELECT
     633        1351 :          atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
     634             : 
     635             :          ! Core energy
     636        1351 :          atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
     637             : 
     638             :          ! Confinement energy
     639        1351 :          IF (atom%potential%confinement) THEN
     640         747 :             atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
     641             :          ELSE
     642         604 :             atom%energy%econfinement = 0._dp
     643             :          END IF
     644             : 
     645             :          ! Hartree Term
     646     3696193 :          jmat%op = 0._dp
     647        1351 :          SELECT CASE (atom%coulomb_integral_type)
     648             :          CASE DEFAULT
     649           0 :             CPABORT("")
     650             :          CASE (do_analytic)
     651          10 :             CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
     652             :          CASE (do_semi_analytic)
     653             :             CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
     654           0 :                                             atom%state%maxl_occ)
     655           0 :             CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     656             :          CASE (do_numeric)
     657        1341 :             CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
     658        1341 :             CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
     659        2692 :             CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
     660             :          END SELECT
     661        1351 :          atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
     662             : 
     663             :          ! Exchange Term
     664        1351 :          IF (need_x) THEN
     665      295546 :             kmata%op = 0._dp
     666      295546 :             kmatb%op = 0._dp
     667         118 :             SELECT CASE (atom%exchange_integral_type)
     668             :             CASE DEFAULT
     669           0 :                CPABORT("")
     670             :             CASE (do_analytic)
     671           4 :                CALL eeri_contract(kmata%op, atom%integrals%eeri, atom%orbitals%pmata, atom%integrals%n)
     672           4 :                CALL eeri_contract(kmatb%op, atom%integrals%eeri, atom%orbitals%pmatb, atom%integrals%n)
     673             :             CASE (do_semi_analytic)
     674           0 :                CALL exchange_semi_analytic(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
     675           0 :                CALL exchange_semi_analytic(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
     676             :             CASE (do_numeric)
     677         114 :                CALL exchange_numeric(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
     678         232 :                CALL exchange_numeric(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
     679             :             END SELECT
     680             :             atom%energy%eexchange = hf_frac*(atom_trace(kmata%op, atom%orbitals%pmata) + &
     681         118 :                                              atom_trace(kmatb%op, atom%orbitals%pmatb))
     682      295546 :             kmata%op = 2._dp*hf_frac*kmata%op
     683      295546 :             kmatb%op = 2._dp*hf_frac*kmatb%op
     684             :          ELSE
     685     3400647 :             kmata%op = 0._dp
     686     3400647 :             kmatb%op = 0._dp
     687        1233 :             atom%energy%eexchange = 0._dp
     688             :          END IF
     689             : 
     690             :          ! XC
     691        1351 :          IF (need_xc) THEN
     692     3558249 :             xcmata%op = 0._dp
     693     3558249 :             xcmatb%op = 0._dp
     694        1335 :             CALL calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
     695             :          ELSE
     696      137944 :             xcmata%op = 0._dp
     697      137944 :             xcmatb%op = 0._dp
     698          16 :             atom%energy%exc = 0._dp
     699             :          END IF
     700             : 
     701        1351 :          IF (lsdpot) THEN
     702             :             atom%energy%elsd = atom_trace(hlsd%op, atom%orbitals%pmata) - &
     703           0 :                                atom_trace(hlsd%op, atom%orbitals%pmatb)
     704           0 :             atom%energy%epseudo = atom%energy%epseudo + atom%energy%elsd
     705           0 :             atom%energy%ecore = atom%energy%ecore + atom%energy%elsd
     706             :          ELSE
     707        1351 :             atom%energy%elsd = 0._dp
     708             :          END IF
     709             : 
     710             :          ! Total energy
     711        1351 :          atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
     712             : 
     713             :          ! Potential energy
     714        1351 :          atom%energy%epot = atom%energy%etot - atom%energy%ekin
     715             : 
     716             :          ! Total HF/KS matrix
     717     7391035 :          fmata%op = hcore%op + hlsd%op + jmat%op + kmata%op + xcmata%op
     718     7391035 :          fmatb%op = hcore%op - hlsd%op + jmat%op + kmatb%op + xcmatb%op
     719             : 
     720             :          ! calculate error matrix
     721             :          CALL err_matrix(xcmata%op, depsa, fmata%op, atom%orbitals%pmata, atom%integrals%utrans, &
     722        1351 :                          atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
     723             :          CALL err_matrix(xcmatb%op, depsb, fmatb%op, atom%orbitals%pmatb, atom%integrals%utrans, &
     724        1351 :                          atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
     725        1351 :          deps = 2._dp*MAX(depsa, depsb)
     726             : 
     727        1351 :          iter = iter + 1
     728             : 
     729        1351 :          IF (iw > 0) THEN
     730         504 :             CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
     731             :          END IF
     732             : 
     733        1351 :          IF (deps < eps_scf) THEN
     734             :             is_converged = .TRUE.
     735             :             EXIT
     736             :          END IF
     737        1234 :          IF (iter >= max_iter) THEN
     738           0 :             IF (iw > 0) THEN
     739           0 :                WRITE (iw, "(A)") " No convergence within maximum number of iterations "
     740             :             END IF
     741             :             EXIT
     742             :          END IF
     743             : 
     744             :          ! update history container and extrapolate KS matrix
     745        1234 :          CALL atom_history_update(historya, atom%orbitals%pmata, fmata%op, xcmata%op, atom%energy%etot, deps)
     746        1234 :          CALL atom_history_update(historyb, atom%orbitals%pmatb, fmatb%op, xcmatb%op, atom%energy%etot, deps)
     747        1234 :          CALL atom_opt_fmat(fmata%op, historya, depsa)
     748        1234 :          CALL atom_opt_fmat(fmatb%op, historyb, depsb)
     749             : 
     750             :          ! Solve HF/KS equations
     751             :          CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
     752        1234 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     753             :          CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
     754        1234 :                           atom%state%maxl_occ, atom%state%maxn_occ)
     755             :          CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
     756        1234 :                          atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
     757             :          CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
     758        1234 :                           atom%state%maxl_occ, atom%state%maxn_occ)
     759     3420454 :          atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
     760             : 
     761             :       END DO !SCF Loop
     762             : 
     763         117 :       IF (iw > 0) THEN
     764          45 :          CALL atom_print_energies(atom, iw)
     765             :       END IF
     766             : 
     767         117 :       CALL atom_history_release(historya)
     768         117 :       CALL atom_history_release(historyb)
     769             : 
     770         117 :       CALL release_opgrid(density)
     771         117 :       CALL release_opgrid(rhoa)
     772         117 :       CALL release_opgrid(rhob)
     773         117 :       CALL release_opgrid(cpot)
     774             : 
     775             :       ! conserve fmata as fmat within atom_type.
     776         117 :       CALL set_atom(atom, fmat=fmata)
     777         117 :       CALL release_opmat(fmatb)
     778         117 :       CALL release_opmat(jmat)
     779         117 :       CALL release_opmat(kmata)
     780         117 :       CALL release_opmat(kmatb)
     781         117 :       CALL release_opmat(xcmata)
     782         117 :       CALL release_opmat(xcmatb)
     783         117 :       CALL release_opmat(hlsd)
     784         117 :       CALL release_opmat(hcore)
     785             : 
     786         117 :       IF (PRESENT(converged)) THEN
     787          22 :          converged = is_converged
     788             :       END IF
     789             : 
     790         117 :       CALL timestop(handle)
     791             : 
     792         117 :    END SUBROUTINE calculate_atom_unrestricted
     793             : 
     794             : END MODULE atom_electronic_structure

Generated by: LCOV version 1.15