LCOV - code coverage report
Current view: top level - src - atom_output.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 260 431 60.3 %
Date: 2025-02-18 08:24:35 Functions: 9 13 69.2 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines that print various information about an atomic kind.
      10             : ! **************************************************************************************************
      11             : MODULE atom_output
      12             :    USE atom_types,                      ONLY: &
      13             :         atom_basis_type, atom_gthpot_type, atom_potential_type, atom_state, atom_type, cgto_basis, &
      14             :         ecp_pseudo, gth_pseudo, gto_basis, lmat, no_pseudo, num_basis, sgp_pseudo, sto_basis, &
      15             :         upf_pseudo
      16             :    USE atom_utils,                      ONLY: get_maxl_occ,&
      17             :                                               get_maxn_occ,&
      18             :                                               get_rho0
      19             :    USE cp_files,                        ONLY: close_file,&
      20             :                                               open_file
      21             :    USE input_constants,                 ONLY: &
      22             :         barrier_conf, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
      23             :         do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_uhf_atom, do_uks_atom, &
      24             :         do_zoramp_atom, poly_conf, xc_none
      25             :    USE input_cp2k_check,                ONLY: xc_functionals_expand
      26             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      27             :                                               section_vals_get_subs_vals2,&
      28             :                                               section_vals_type,&
      29             :                                               section_vals_val_get
      30             :    USE kinds,                           ONLY: default_string_length,&
      31             :                                               dp
      32             :    USE mathconstants,                   ONLY: dfac,&
      33             :                                               pi,&
      34             :                                               rootpi
      35             :    USE periodic_table,                  ONLY: ptable
      36             :    USE physcon,                         ONLY: evolt
      37             :    USE xc_derivatives,                  ONLY: xc_functional_get_info
      38             :    USE xc_libxc,                        ONLY: libxc_check_existence_in_libxc,&
      39             :                                               libxc_get_reference_length
      40             :    USE xmgrace,                         ONLY: xm_graph_data,&
      41             :                                               xm_graph_info,&
      42             :                                               xm_write_defaults,&
      43             :                                               xm_write_frame,&
      44             :                                               xm_write_frameport
      45             : #include "./base/base_uses.f90"
      46             : 
      47             :    IMPLICIT NONE
      48             : 
      49             :    PRIVATE
      50             : 
      51             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_output'
      52             : 
      53             :    PUBLIC :: atom_print_state, atom_print_energies, atom_print_iteration, &
      54             :              atom_print_basis, atom_print_method, atom_print_info, atom_print_potential, &
      55             :              atom_print_basis_file, atom_write_pseudo_param, atom_print_orbitals, &
      56             :              atom_print_zmp_iteration
      57             : 
      58             : CONTAINS
      59             : 
      60             : ! **************************************************************************************************
      61             : !> \brief Print an information string related to the atomic kind.
      62             : !> \param zval  atomic number
      63             : !> \param info  information string
      64             : !> \param iw    output file unit
      65             : !> \par History
      66             : !>    * 09.2008 created [Juerg Hutter]
      67             : ! **************************************************************************************************
      68         180 :    SUBROUTINE atom_print_info(zval, info, iw)
      69             :       INTEGER, INTENT(IN)                                :: zval
      70             :       CHARACTER(len=*), INTENT(IN)                       :: info
      71             :       INTEGER, INTENT(IN)                                :: iw
      72             : 
      73             :       WRITE (iw, '(/," ",A,T40,A," [",A,"]",T62,"Atomic number:",T78,I3,/)') &
      74         180 :          ADJUSTL(TRIM(info)), TRIM(ptable(zval)%name), TRIM(ptable(zval)%symbol), zval
      75             : 
      76         180 :    END SUBROUTINE atom_print_info
      77             : 
      78             : ! **************************************************************************************************
      79             : !> \brief Print information about electronic state.
      80             : !> \param state  electronic state
      81             : !> \param iw     output file unit
      82             : !> \par History
      83             : !>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
      84             : !>    * 11.2009 print multiplicity [Juerg Hutter]
      85             : !>    * 08.2008 created [Juerg Hutter]
      86             : ! **************************************************************************************************
      87        2258 :    SUBROUTINE atom_print_state(state, iw)
      88             :       TYPE(atom_state)                                   :: state
      89             :       INTEGER, INTENT(IN)                                :: iw
      90             : 
      91             :       CHARACTER(LEN=1), DIMENSION(0:7), PARAMETER :: &
      92             :          label = (/"S", "P", "D", "F", "G", "H", "I", "K"/)
      93             : 
      94             :       INTEGER                                            :: j, l, mc, mlc, mlo, mm(0:lmat), mo
      95             : 
      96             :       CPASSERT(lmat <= 7)
      97        2258 :       WRITE (iw, '(/,T2,A)') "Electronic structure"
      98      160318 :       WRITE (iw, '(T5,A,T71,F10.2)') "Total number of core electrons", SUM(state%core)
      99      160318 :       WRITE (iw, '(T5,A,T71,F10.2)') "Total number of valence electrons", SUM(state%occ)
     100      160318 :       WRITE (iw, '(T5,A,T71,F10.2)') "Total number of electrons", SUM(state%occ + state%core)
     101        4471 :       SELECT CASE (state%multiplicity)
     102             :       CASE (-1)
     103        2213 :          WRITE (iw, '(T5,A,T68,A)') "Multiplicity", "not specified"
     104             :       CASE (-2)
     105          10 :          WRITE (iw, '(T5,A,T72,A)') "Multiplicity", "high spin"
     106             :       CASE (-3)
     107           0 :          WRITE (iw, '(T5,A,T73,A)') "Multiplicity", "low spin"
     108             :       CASE (1)
     109          22 :          WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "singlet"
     110             :       CASE (2)
     111          10 :          WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "doublet"
     112             :       CASE (3)
     113           3 :          WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "triplet"
     114             :       CASE (4)
     115           0 :          WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quartet"
     116             :       CASE (5)
     117           0 :          WRITE (iw, '(T5,A,T74,A)') "Multiplicity", "quintet"
     118             :       CASE (6)
     119           0 :          WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "sextet"
     120             :       CASE (7)
     121        2258 :          WRITE (iw, '(T5,A,T75,A)') "Multiplicity", "septet"
     122             :       CASE DEFAULT
     123             :       END SELECT
     124             : 
     125        2258 :       mlo = get_maxl_occ(state%occ)
     126        2258 :       mlc = get_maxl_occ(state%core)
     127        2258 :       mm = get_maxn_occ(state%core)
     128             : 
     129        2258 :       IF (state%multiplicity == -1) THEN
     130        5722 :          DO l = 0, MAX(mlo, mlc)
     131        3509 :             mo = state%maxn_occ(l)
     132       40812 :             IF (SUM(state%core(l, :)) == 0) THEN
     133        5771 :                WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occ(l, j), j=1, mo)
     134             :             ELSE
     135        1086 :                mc = mm(l)
     136        2448 :                CPASSERT(SUM(state%occ(l, 1:mc)) == 0)
     137        2448 :                WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (state%core(l, j), j=1, mc)
     138        2141 :                WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occ(l, j), j=mc + 1, mc + mo)
     139             :             END IF
     140             :          END DO
     141             :       ELSE
     142          45 :          WRITE (iw, '(T5,A)') "Alpha Electrons"
     143         122 :          DO l = 0, MAX(mlo, mlc)
     144          77 :             mo = state%maxn_occ(l)
     145         892 :             IF (SUM(state%core(l, :)) == 0) THEN
     146         104 :                WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occa(l, j), j=1, mo)
     147             :             ELSE
     148          35 :                mc = mm(l)
     149          86 :                WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
     150          60 :                WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occa(l, j), j=1, mo)
     151             :             END IF
     152             :          END DO
     153          45 :          WRITE (iw, '(T5,A)') "Beta Electrons"
     154         122 :          DO l = 0, MAX(mlo, mlc)
     155          77 :             mo = state%maxn_occ(l)
     156         892 :             IF (SUM(state%core(l, :)) == 0) THEN
     157         104 :                WRITE (iw, '(A5,T10,10F6.2)') label(l), (state%occb(l, j), j=1, mo)
     158             :             ELSE
     159          35 :                mc = mm(l)
     160          86 :                WRITE (iw, ADVANCE="no", FMT='(A5,T9,A1,10F6.2)') label(l), "[", (0.5_dp*state%core(l, j), j=1, mc)
     161          60 :                WRITE (iw, FMT='(A1,F5.2,10F6.2)') "]", (state%occb(l, j), j=1, mo)
     162             :             END IF
     163             :          END DO
     164             :       END IF
     165        2258 :       WRITE (iw, *)
     166             : 
     167        2258 :    END SUBROUTINE atom_print_state
     168             : 
     169             : ! **************************************************************************************************
     170             : !> \brief Print energy components.
     171             : !> \param atom  information about the atomic kind
     172             : !> \param iw    output file unit
     173             : !> \par History
     174             : !>    * 05.2010 print virial coefficient [Juerg Hutter]
     175             : !>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
     176             : !>    * 09.2008 print orbital energies [Juerg Hutter]
     177             : !>    * 08.2008 created [Juerg Hutter]
     178             : ! **************************************************************************************************
     179        4516 :    SUBROUTINE atom_print_energies(atom, iw)
     180             :       TYPE(atom_type)                                    :: atom
     181             :       INTEGER, INTENT(IN)                                :: iw
     182             : 
     183             :       INTEGER                                            :: i, l, n
     184             :       REAL(KIND=dp)                                      :: drho
     185             : 
     186        2258 :       WRITE (iw, '(/,A,T36,A,T61,F20.12)') " Energy components [Hartree]", &
     187        4516 :          "    Total Energy ::", atom%energy%etot
     188        2258 :       WRITE (iw, '(T36,A,T61,F20.12)') "     Band Energy ::", atom%energy%eband
     189        2258 :       WRITE (iw, '(T36,A,T61,F20.12)') "  Kinetic Energy ::", atom%energy%ekin
     190        2258 :       WRITE (iw, '(T36,A,T61,F20.12)') "Potential Energy ::", atom%energy%epot
     191        2258 :       IF (atom%energy%ekin /= 0.0_dp) THEN
     192        2245 :          WRITE (iw, '(T36,A,T61,F20.12)') "   Virial (-V/T) ::", -atom%energy%epot/atom%energy%ekin
     193             :       END IF
     194        2258 :       WRITE (iw, '(T36,A,T61,F20.12)') "     Core Energy ::", atom%energy%ecore
     195        2258 :       IF (atom%energy%exc /= 0._dp) &
     196        2227 :          WRITE (iw, '(T36,A,T61,F20.12)') "       XC Energy ::", atom%energy%exc
     197        2258 :       WRITE (iw, '(T36,A,T61,F20.12)') "  Coulomb Energy ::", atom%energy%ecoulomb
     198        2258 :       IF (atom%energy%eexchange /= 0._dp) &
     199          46 :          WRITE (iw, '(T34,A,T61,F20.12)') "HF Exchange Energy ::", atom%energy%eexchange
     200        2258 :       IF (atom%potential%ppot_type /= NO_PSEUDO) THEN
     201        1868 :          WRITE (iw, '(T20,A,T61,F20.12)') "    Total Pseudopotential Energy ::", atom%energy%epseudo
     202        1868 :          WRITE (iw, '(T20,A,T61,F20.12)') "    Local Pseudopotential Energy ::", atom%energy%eploc
     203        1868 :          IF (atom%energy%elsd /= 0._dp) &
     204           0 :             WRITE (iw, '(T20,A,T61,F20.12)') "     Local Spin-potential Energy ::", atom%energy%elsd
     205        1868 :          WRITE (iw, '(T20,A,T61,F20.12)') " Nonlocal Pseudopotential Energy ::", atom%energy%epnl
     206             :       END IF
     207        2258 :       IF (atom%potential%confinement) THEN
     208        1862 :          WRITE (iw, '(T36,A,T61,F20.12)') "     Confinement ::", atom%energy%econfinement
     209             :       END IF
     210             : 
     211        2258 :       IF (atom%state%multiplicity == -1) THEN
     212        2213 :          WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T49,A,T71,A,/)') " Orbital energies", &
     213        4426 :             "State", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
     214        5832 :          DO l = 0, atom%state%maxl_calc
     215        3619 :             n = atom%state%maxn_calc(l)
     216        8137 :             DO i = 1, n
     217             :                WRITE (iw, '(T23,I2,T30,I1,T36,F10.3,T46,F15.6,T66,F15.6)') &
     218        8137 :                   i, l, atom%state%occupation(l, i), atom%orbitals%ener(i, l), atom%orbitals%ener(i, l)*evolt
     219             :             END DO
     220        5832 :             IF (n > 0) WRITE (iw, *)
     221             :          END DO
     222             :       ELSE
     223          45 :          WRITE (iw, '(/,A,T20,A,T30,A,T36,A,T42,A,T55,A,T71,A,/)') " Orbital energies", &
     224          90 :             "State", "Spin", "L", "Occupation", "Energy[a.u.]", "Energy[eV]"
     225         162 :          DO l = 0, atom%state%maxl_calc
     226         117 :             n = atom%state%maxn_calc(l)
     227         208 :             DO i = 1, n
     228             :                WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
     229         208 :                   i, "alpha", l, atom%state%occa(l, i), atom%orbitals%enera(i, l), atom%orbitals%enera(i, l)*evolt
     230             :             END DO
     231         208 :             DO i = 1, n
     232             :                WRITE (iw, '(T23,I2,T29,A,T36,I1,T42,F10.3,T52,F15.6,T68,F13.6)') &
     233         208 :                   i, " beta", l, atom%state%occb(l, i), atom%orbitals%enerb(i, l), atom%orbitals%enerb(i, l)*evolt
     234             :             END DO
     235         162 :             IF (n > 0) WRITE (iw, *)
     236             :          END DO
     237             :       END IF
     238             : 
     239        2258 :       CALL get_rho0(atom, drho)
     240        2258 :       WRITE (iw, '(/,A,T66,F15.6)') " Total Electron Density at R=0: ", drho
     241             : 
     242        2258 :    END SUBROUTINE atom_print_energies
     243             : 
     244             : ! **************************************************************************************************
     245             : !> \brief Printing of the atomic iterations when ZMP is active.
     246             : !> \param iter  current iteration number
     247             : !> \param deps  convergence
     248             : !> \param atom  intormation about the atomic kind
     249             : !> \param iw    output file unit
     250             : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
     251             : ! **************************************************************************************************
     252           0 :    SUBROUTINE atom_print_zmp_iteration(iter, deps, atom, iw)
     253             :       INTEGER, INTENT(IN)                                :: iter
     254             :       REAL(dp), INTENT(IN)                               :: deps
     255             :       TYPE(atom_type), INTENT(IN)                        :: atom
     256             :       INTEGER, INTENT(IN)                                :: iw
     257             : 
     258           0 :       IF (iter == 1) THEN
     259             :          WRITE (iw, '(/," ",79("*"),/,T33,"Integral",T48,"Integral",/,T3,A,T16,A,T33,A,T46,A,T69,A/," ",79("*"))') &
     260           0 :             "Iteration", "Convergence", "rho diff.", "rho*v_xc[au]", "Energy[au]"
     261             :       END IF
     262           0 :       WRITE (iw, '(T3,I9,T15,G13.6,T30,G13.6,T46,G13.6,T61,F20.12)') iter, deps, atom%rho_diff_integral, &
     263           0 :          atom%energy%exc, atom%energy%etot
     264             : 
     265           0 :    END SUBROUTINE atom_print_zmp_iteration
     266             : 
     267             : ! **************************************************************************************************
     268             : !> \brief Print convergence information.
     269             : !> \param iter  current iteration number
     270             : !> \param deps  convergency
     271             : !> \param etot  total energy
     272             : !> \param iw    output file unit
     273             : !> \par History
     274             : !>    * 08.2008 created [Juerg Hutter]
     275             : ! **************************************************************************************************
     276       10025 :    SUBROUTINE atom_print_iteration(iter, deps, etot, iw)
     277             :       INTEGER, INTENT(IN)                                :: iter
     278             :       REAL(dp), INTENT(IN)                               :: deps, etot
     279             :       INTEGER, INTENT(IN)                                :: iw
     280             : 
     281       10025 :       IF (iter == 1) THEN
     282             :          WRITE (iw, '(/," ",79("*"),/,T19,A,T38,A,T70,A,/," ",79("*"))') &
     283        2258 :             "Iteration", "Convergence", "Energy [au]"
     284             :       END IF
     285       10025 :       WRITE (iw, '(T20,i8,T34,G14.6,T61,F20.12)') iter, deps, etot
     286             : 
     287       10025 :    END SUBROUTINE atom_print_iteration
     288             : 
     289             : ! **************************************************************************************************
     290             : !> \brief Print atomic basis set.
     291             : !> \param atom_basis  atomic basis set
     292             : !> \param iw          output file unit
     293             : !> \param title       header to print on top of the basis set
     294             : !> \par History
     295             : !>    * 09.2008 created [Juerg Hutter]
     296             : ! **************************************************************************************************
     297          22 :    SUBROUTINE atom_print_basis(atom_basis, iw, title)
     298             :       TYPE(atom_basis_type)                              :: atom_basis
     299             :       INTEGER, INTENT(IN)                                :: iw
     300             :       CHARACTER(len=*)                                   :: title
     301             : 
     302             :       INTEGER                                            :: i, j, l
     303             : 
     304          22 :       WRITE (iw, '(/,A)') TRIM(title)
     305          39 :       SELECT CASE (atom_basis%basis_type)
     306             :       CASE (GTO_BASIS)
     307          17 :          IF (atom_basis%geometrical) THEN
     308          15 :             WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
     309          15 :             WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
     310          30 :                " Proportionality factor: ", atom_basis%cval
     311             :          ELSE
     312           2 :             WRITE (iw, '(/," ",21("*"),A,21("*"))') " Uncontracted Gaussian Type Orbitals "
     313             :          END IF
     314         119 :          DO l = 0, lmat
     315         119 :             IF (atom_basis%nbas(l) > 0) THEN
     316          14 :                SELECT CASE (l)
     317             :                CASE DEFAULT
     318             :                   WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
     319         280 :                      "X Exponents: ", (i, atom_basis%am(i, l), i=1, atom_basis%nbas(l))
     320             :                CASE (0)
     321             :                   WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
     322         440 :                      "s Exponents: ", (i, atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
     323             :                CASE (1)
     324             :                   WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
     325         408 :                      "p Exponents: ", (i, atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
     326             :                CASE (2)
     327             :                   WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
     328         376 :                      "d Exponents: ", (i, atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
     329             :                CASE (3)
     330             :                   WRITE (iw, '(/,T2,A,(T30,I5,T51,F30.8))') &
     331         395 :                      "f Exponents: ", (i, atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
     332             :                END SELECT
     333             :             END IF
     334             :          END DO
     335          17 :          WRITE (iw, '(" ",79("*"))')
     336             :       CASE (CGTO_BASIS)
     337           1 :          WRITE (iw, '(/," ",22("*"),A,22("*"))') " Contracted Gaussian Type Orbitals "
     338           7 :          DO l = 0, lmat
     339           7 :             IF (atom_basis%nbas(l) > 0) THEN
     340           2 :                IF (l == 0) WRITE (iw, '(A)') " s Functions"
     341           2 :                IF (l == 1) WRITE (iw, '(A)') " p Functions"
     342           2 :                IF (l == 2) WRITE (iw, '(A)') " d Functions"
     343           2 :                IF (l == 3) WRITE (iw, '(A)') " f Functions"
     344           2 :                IF (l >= 3) WRITE (iw, '(A)') " x Functions"
     345          11 :                DO i = 1, atom_basis%nprim(l)
     346             :                   WRITE (iw, '(F15.6,5(T21,6F10.6,/))') &
     347          35 :                      atom_basis%am(i, l), (atom_basis%cm(i, j, l), j=1, atom_basis%nbas(l))
     348             :                END DO
     349             :             END IF
     350             :          END DO
     351           1 :          WRITE (iw, '(" ",79("*"))')
     352             :       CASE (STO_BASIS)
     353           4 :          WRITE (iw, '(/," ",28("*"),A,29("*"))') " Slater Type Orbitals "
     354          28 :          DO l = 0, lmat
     355          39 :             DO i = 1, atom_basis%nbas(l)
     356          24 :                SELECT CASE (l)
     357             :                CASE DEFAULT
     358           0 :                   WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, l), "X Exponent :", atom_basis%as(i, l)
     359             :                CASE (0)
     360           8 :                   WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 0), "S Exponent :", atom_basis%as(i, 0)
     361             :                CASE (1)
     362           3 :                   WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 1), "P Exponent :", atom_basis%as(i, 1)
     363             :                CASE (2)
     364           0 :                   WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 2), "D Exponent :", atom_basis%as(i, 2)
     365             :                CASE (3)
     366          11 :                   WRITE (iw, '(T10,I1,A,T40,F25.12)') atom_basis%ns(i, 3), "F Exponent :", atom_basis%as(i, 3)
     367             :                END SELECT
     368             :             END DO
     369             :          END DO
     370           4 :          WRITE (iw, '(" ",79("*"))')
     371             :       CASE (NUM_BASIS)
     372           0 :          CPABORT("")
     373             :       CASE DEFAULT
     374          22 :          CPABORT("")
     375             :       END SELECT
     376             : 
     377          22 :    END SUBROUTINE atom_print_basis
     378             : 
     379             : ! **************************************************************************************************
     380             : !> \brief Print the optimized atomic basis set into a file.
     381             : !> \param atom_basis  atomic basis set
     382             : !> \param wfn ...
     383             : !> \par History
     384             : !>    * 11.2016 revised output format [Matthias Krack]
     385             : !>    * 11.2011 Slater basis functions [Juerg Hutter]
     386             : !>    * 03.2011 created [Juerg Hutter]
     387             : !> \note  The basis set is stored as the file 'OPT_BASIS' inside the current working directory.
     388             : !>        It may be a good idea, however, to specify the name of this file via some input section.
     389             : ! **************************************************************************************************
     390           5 :    SUBROUTINE atom_print_basis_file(atom_basis, wfn)
     391             :       TYPE(atom_basis_type)                              :: atom_basis
     392             :       REAL(KIND=dp), DIMENSION(:, :, 0:), OPTIONAL       :: wfn
     393             : 
     394             :       INTEGER                                            :: i, im, iw, l
     395             :       REAL(KIND=dp)                                      :: expzet, prefac, zeta
     396             : 
     397           5 :       CALL open_file(file_name="OPT_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
     398           7 :       SELECT CASE (atom_basis%basis_type)
     399             :       CASE (GTO_BASIS)
     400           2 :          IF (atom_basis%geometrical) THEN
     401           0 :             WRITE (iw, '(/," ",21("*"),A,22("*"))') " Geometrical Gaussian Type Orbitals "
     402           0 :             WRITE (iw, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", atom_basis%aval, &
     403           0 :                " Proportionality factor: ", atom_basis%cval
     404             :          ELSE
     405           2 :             WRITE (iw, '(T3,A)') "BASIS_TYPE GAUSSIAN"
     406             :          END IF
     407          14 :          DO l = 0, lmat
     408          14 :             IF (atom_basis%nbas(l) > 0) THEN
     409           0 :                SELECT CASE (l)
     410             :                CASE DEFAULT
     411             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     412           0 :                      "X_EXPONENTS ", (atom_basis%am(i, l), i=1, atom_basis%nbas(l))
     413             :                CASE (0)
     414             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     415          14 :                      "S_EXPONENTS ", (atom_basis%am(i, 0), i=1, atom_basis%nbas(0))
     416             :                CASE (1)
     417             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     418          14 :                      "P_EXPONENTS ", (atom_basis%am(i, 1), i=1, atom_basis%nbas(1))
     419             :                CASE (2)
     420             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     421          14 :                      "D_EXPONENTS ", (atom_basis%am(i, 2), i=1, atom_basis%nbas(2))
     422             :                CASE (3)
     423             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     424           6 :                      "F_EXPONENTS ", (atom_basis%am(i, 3), i=1, atom_basis%nbas(3))
     425             :                END SELECT
     426             :             END IF
     427             :          END DO
     428             :       CASE (CGTO_BASIS)
     429           0 :          CPABORT("")
     430             :       CASE (STO_BASIS)
     431           3 :          WRITE (iw, '(T3,A)') "BASIS_TYPE SLATER"
     432          21 :          DO l = 0, lmat
     433          21 :             IF (atom_basis%nbas(l) > 0) THEN
     434           0 :                SELECT CASE (l)
     435             :                CASE DEFAULT
     436             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     437           0 :                      "X_EXPONENTS ", (atom_basis%as(i, l), i=1, atom_basis%nbas(l))
     438             :                   WRITE (iw, '(T3,A,60I3)') &
     439           0 :                      "X_QUANTUM_NUMBERS ", (atom_basis%ns(i, l), i=1, atom_basis%nbas(l))
     440             :                CASE (0)
     441             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     442          10 :                      "S_EXPONENTS ", (atom_basis%as(i, 0), i=1, atom_basis%nbas(0))
     443             :                   WRITE (iw, '(T3,A,60I3)') &
     444          10 :                      "S_QUANTUM_NUMBERS ", (atom_basis%ns(i, 0), i=1, atom_basis%nbas(0))
     445             :                CASE (1)
     446             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     447           3 :                      "P_EXPONENTS ", (atom_basis%as(i, 1), i=1, atom_basis%nbas(1))
     448             :                   WRITE (iw, '(T3,A,60I3)') &
     449           3 :                      "P_QUANTUM_NUMBERS ", (atom_basis%ns(i, 1), i=1, atom_basis%nbas(1))
     450             :                CASE (2)
     451             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     452           0 :                      "D_EXPONENTS ", (atom_basis%as(i, 2), i=1, atom_basis%nbas(2))
     453             :                   WRITE (iw, '(T3,A,60I3)') &
     454           0 :                      "D_QUANTUM_NUMBERS ", (atom_basis%ns(i, 2), i=1, atom_basis%nbas(2))
     455             :                CASE (3)
     456             :                   WRITE (iw, '(T3,A,(T15,F20.8,:," \"))') &
     457           0 :                      "F_EXPONENTS ", (atom_basis%as(i, 3), i=1, atom_basis%nbas(3))
     458             :                   WRITE (iw, '(T3,A,60I3)') &
     459           4 :                      "F_QUANTUM_NUMBERS ", (atom_basis%ns(i, 3), i=1, atom_basis%nbas(3))
     460             :                END SELECT
     461             :             END IF
     462             :          END DO
     463             :       CASE (NUM_BASIS)
     464           0 :          CPABORT("")
     465             :       CASE DEFAULT
     466           5 :          CPABORT("")
     467             :       END SELECT
     468             : 
     469           5 :       IF (PRESENT(wfn)) THEN
     470           7 :          SELECT CASE (atom_basis%basis_type)
     471             :          CASE DEFAULT
     472             :          CASE (GTO_BASIS)
     473           5 :             IF (.NOT. atom_basis%geometrical) THEN
     474           2 :                WRITE (iw, '(/,T3,A)') "ORBITAL COEFFICENTS (Quickstep normalization)"
     475           2 :                im = MIN(6, SIZE(wfn, 2))
     476          14 :                DO l = 0, lmat
     477          14 :                   IF (atom_basis%nbas(l) > 0) THEN
     478           6 :                      WRITE (iw, '(T3,A,I3)') "L Quantum Number:", l
     479             :                      ! Quickstep normalization
     480           6 :                      expzet = 0.25_dp*REAL(2*l + 3, dp)
     481           6 :                      prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
     482          42 :                      DO i = 1, atom_basis%nbas(l)
     483          36 :                         zeta = (2._dp*atom_basis%am(i, l))**expzet
     484          78 :                         WRITE (iw, '(T5,F14.8,2x,6F12.8)') atom_basis%am(i, l), wfn(i, 1:im, l)*prefac/zeta
     485             :                      END DO
     486             :                   END IF
     487             :                END DO
     488             :             END IF
     489             :          END SELECT
     490             :       END IF
     491             : 
     492           5 :       CALL close_file(unit_number=iw)
     493             : 
     494           5 :    END SUBROUTINE atom_print_basis_file
     495             : 
     496             : ! **************************************************************************************************
     497             : !> \brief Print information about the electronic structure method in use.
     498             : !> \param atom  information about the atomic kind
     499             : !> \param iw    output file unit
     500             : !> \par History
     501             : !>    * 09.2015 direct use of the LibXC Fortran interface [Andreas Gloess]
     502             : !>    * 10.2012 LibXC interface [Fabien Tran]
     503             : !>    * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
     504             : !>    * 04.2009 print geometrical Gaussian type orbitals [Juerg Hutter]
     505             : !>    * 09.2008 new subroutine's prototype; print relativistic methods [Juerg Hutter]
     506             : !>    * 09.2008 created [Juerg Hutter]
     507             : ! **************************************************************************************************
     508         366 :    SUBROUTINE atom_print_method(atom, iw)
     509             :       TYPE(atom_type)                                    :: atom
     510             :       INTEGER, INTENT(IN)                                :: iw
     511             : 
     512             :       CHARACTER(len=160)                                 :: shortform
     513             :       CHARACTER(LEN=20)                                  :: tmpStr
     514         366 :       CHARACTER(len=:), ALLOCATABLE                      :: reference
     515             :       INTEGER                                            :: ifun, il, meth, myfun, reltyp
     516             :       LOGICAL                                            :: lsd
     517             :       TYPE(section_vals_type), POINTER                   :: xc_fun, xc_fun_section, xc_section
     518             : 
     519         366 :       NULLIFY (xc_fun, xc_fun_section, xc_section)
     520             : 
     521         366 :       meth = atom%method_type
     522             : 
     523         366 :       xc_section => atom%xc_section
     524         366 :       xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     525           0 :       SELECT CASE (meth)
     526             :       CASE DEFAULT
     527           0 :          CPABORT("")
     528             :       CASE (do_rks_atom)
     529         328 :          CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
     530             :       CASE (do_uks_atom)
     531          34 :          CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
     532             :       CASE (do_rhf_atom)
     533          34 :          myfun = xc_none
     534             :       CASE (do_uhf_atom)
     535           4 :          myfun = xc_none
     536             :       CASE (do_rohf_atom)
     537         366 :          myfun = xc_none
     538             :       END SELECT
     539             : 
     540           0 :       SELECT CASE (meth)
     541             :       CASE DEFAULT
     542           0 :          CPABORT("")
     543             :       CASE (do_rks_atom)
     544         294 :          IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Kohn-Sham Calculation')")
     545             :       CASE (do_uks_atom)
     546          34 :          IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Unrestricted Kohn-Sham Calculation')")
     547             :       CASE (do_rhf_atom)
     548          34 :          IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Hartree-Fock Calculation')")
     549             :       CASE (do_uhf_atom)
     550           4 :          IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Unrestricted Hartree-Fock Calculation')")
     551             :       CASE (do_rohf_atom)
     552         328 :          IF (iw > 0) WRITE (iw, fmt="(/,' METHOD    | Restricted Open-Shell Kohn-Sham Calculation')")
     553             :       END SELECT
     554             : 
     555             :       ! zmp
     556         366 :       IF (atom%do_zmp) THEN
     557           0 :          IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Method on atomic radial density')")
     558           0 :          IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Lambda : ',F5.1)") atom%lambda
     559           0 :          IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Reading external density : ',A20)") atom%ext_file
     560           0 :          IF (atom%dm) THEN
     561           0 :             IF (iw > 0) WRITE (iw, fmt="(' ZMP       | The file is in the form of a density matrix')")
     562             :          ELSE
     563           0 :             IF (iw > 0) WRITE (iw, fmt="(' ZMP       | The file is in the form of a linear density')")
     564             :          END IF
     565           0 :          IF (atom%doread) THEN
     566           0 :             IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Restarting calculation from ',A20,' file if present')") atom%zmp_restart_file
     567             :          END IF
     568         366 :       ELSE IF (atom%read_vxc) THEN
     569           0 :          IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Calculating density from external V_xc')")
     570           0 :          IF (iw > 0) WRITE (iw, fmt="(' ZMP       | Reading external v_xc file : ',A20)") atom%ext_vxc_file
     571             :       END IF
     572             : 
     573         366 :       IF (atom%pp_calc) THEN
     574          78 :          IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Nonrelativistic Calculation')")
     575             :       ELSE
     576         288 :          reltyp = atom%relativistic
     577             : 
     578           0 :          SELECT CASE (reltyp)
     579             :          CASE DEFAULT
     580           0 :             CPABORT("")
     581             :          CASE (do_nonrel_atom)
     582         220 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Nonrelativistic Calculation')")
     583             :          CASE (do_zoramp_atom)
     584          14 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using ZORA(MP)')")
     585             :          CASE (do_sczoramp_atom)
     586           2 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using scaled ZORA(MP)')")
     587             :          CASE (do_dkh0_atom)
     588           2 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 0th order')")
     589           2 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using kietic energy scaling')")
     590             :          CASE (do_dkh1_atom)
     591           2 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 1st order')")
     592           2 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Foldy-Wouthuysen transformation')")
     593             :          CASE (do_dkh2_atom)
     594          16 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 2nd order')")
     595             :          CASE (do_dkh3_atom)
     596         288 :             IF (iw > 0) WRITE (iw, fmt="(' METHOD    | Relativistic Calculation using Douglas-Kroll 3rd order')")
     597             :          END SELECT
     598             :       END IF
     599             : 
     600         366 :       lsd = (meth == do_uks_atom)
     601             : 
     602         366 :       IF (myfun /= xc_none) THEN
     603         326 :          CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", c_val=tmpStr)
     604         326 :          IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| ROUTINE=',a)") TRIM(tmpStr)
     605         326 :          CALL xc_functionals_expand(xc_fun_section, xc_section)
     606         326 :          IF (iw > 0) THEN
     607         163 :             ifun = 0
     608         232 :             DO
     609         395 :                ifun = ifun + 1
     610         395 :                xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
     611         395 :                IF (.NOT. ASSOCIATED(xc_fun)) EXIT
     612         232 :                IF (libxc_check_existence_in_libxc(xc_fun)) THEN
     613           3 :                   ALLOCATE (CHARACTER(LEN=libxc_get_reference_length(xc_fun, lsd)) :: reference)
     614             :                ELSE
     615         229 :                   ALLOCATE (CHARACTER(LEN=20*default_string_length) :: reference)
     616             :                END IF
     617         232 :                CALL xc_functional_get_info(xc_fun, lsd=lsd, reference=reference, shortform=shortform)
     618             :                WRITE (iw, fmt="(' FUNCTIONAL| ',a,':')") &
     619         232 :                   TRIM(xc_fun%section%name)
     620         630 :                DO il = 1, LEN_TRIM(reference), 67
     621         630 :                   WRITE (iw, fmt="(' FUNCTIONAL| ',a67)") reference(il:)
     622             :                END DO
     623         395 :                DEALLOCATE (reference)
     624             :             END DO
     625             :          END IF
     626             :       ELSE
     627          40 :          IF (iw > 0) WRITE (iw, fmt="(' FUNCTIONAL| NO EXCHANGE-CORRELATION FUNCTIONAL USED.')")
     628             :       END IF
     629             : 
     630         366 :    END SUBROUTINE atom_print_method
     631             : 
     632             : ! **************************************************************************************************
     633             : !> \brief Print information about the pseudo-potential.
     634             : !> \param potential pseudo-potential
     635             : !> \param iw        output file unit
     636             : !> \par History
     637             : !>    * 05.2017 SGP pseudo-potentials [Juerg Hutter]
     638             : !>    * 02.2016 pseudo-potential in Quantum Espresso UPF format [Juerg Hutter]
     639             : !>    * 01.2016 new confinement potential form [Juerg Hutter]
     640             : !>    * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
     641             : !>    * 05.2009 GTH pseudo-potential [Juerg Hutter]
     642             : !>    * 09.2008 created [Juerg Hutter]
     643             : ! **************************************************************************************************
     644           3 :    SUBROUTINE atom_print_potential(potential, iw)
     645             :       TYPE(atom_potential_type)                          :: potential
     646             :       INTEGER, INTENT(IN)                                :: iw
     647             : 
     648             :       CHARACTER(len=60)                                  :: pline
     649             :       INTEGER                                            :: i, j, k, l
     650             : 
     651           3 :       SELECT CASE (potential%ppot_type)
     652             :       CASE (no_pseudo)
     653           0 :          WRITE (iw, '(/," ",28("*"),A,27("*"))') " All Electron Potential "
     654             :       CASE (gth_pseudo)
     655           0 :          WRITE (iw, '(/," ",29("*"),A,29("*"))') " GTH Pseudopotential "
     656           0 :          WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%gth_pot%zion
     657           0 :          WRITE (iw, '(T10,A,T66,F15.6)') " Rc ", potential%gth_pot%rc
     658           0 :          WRITE (pline, '(5F12.6)') (potential%gth_pot%cl(i), i=1, potential%gth_pot%ncl)
     659           0 :          WRITE (iw, '(T10,A,T21,A60)') " C1 C2 ... ", ADJUSTR(pline)
     660           0 :          IF (potential%gth_pot%lpotextended) THEN
     661           0 :             DO k = 1, potential%gth_pot%nexp_lpot
     662           0 :                WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LPot: rc=", potential%gth_pot%alpha_lpot(k), &
     663           0 :                   "CX=", (potential%gth_pot%cval_lpot(i, k), i=1, potential%gth_pot%nct_lpot(k))
     664             :             END DO
     665             :          END IF
     666           0 :          IF (potential%gth_pot%nlcc) THEN
     667           0 :             DO k = 1, potential%gth_pot%nexp_nlcc
     668           0 :                WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_nlcc(k), &
     669           0 :                   "CX=", (potential%gth_pot%cval_nlcc(i, k)*4.0_dp*pi, i=1, potential%gth_pot%nct_nlcc(k))
     670             :             END DO
     671             :          END IF
     672           0 :          IF (potential%gth_pot%lsdpot) THEN
     673           0 :             DO k = 1, potential%gth_pot%nexp_lsd
     674           0 :                WRITE (iw, '(T10,A,F10.6,T38,A,4F10.6)') " LSDPot: rc=", potential%gth_pot%alpha_lsd(k), &
     675           0 :                   "CX=", (potential%gth_pot%cval_lsd(i, k), i=1, potential%gth_pot%nct_lsd(k))
     676             :             END DO
     677             :          END IF
     678           0 :          DO l = 0, lmat
     679           0 :             IF (potential%gth_pot%nl(l) > 0) THEN
     680           0 :                WRITE (iw, '(T10,A,T76,I5)') " Angular momentum ", l
     681           0 :                WRITE (iw, '(T10,A,T66,F15.6)') " Rcnl ", potential%gth_pot%rcnl(l)
     682           0 :                WRITE (iw, '(T10,A,T76,I5)') " Nl ", potential%gth_pot%nl(l)
     683           0 :                WRITE (pline, '(5F12.6)') (potential%gth_pot%hnl(1, j, l), j=1, potential%gth_pot%nl(l))
     684           0 :                WRITE (iw, '(T10,A,T21,A60)') " Hnl ", ADJUSTR(pline)
     685           0 :                DO i = 2, potential%gth_pot%nl(l)
     686           0 :                   WRITE (pline, '(T21,5F12.6)') (potential%gth_pot%hnl(i, j, l), j=i, potential%gth_pot%nl(l))
     687           0 :                   WRITE (iw, '(T21,A60)') ADJUSTR(pline)
     688             :                END DO
     689             :             END IF
     690             :          END DO
     691             :       CASE (upf_pseudo)
     692           0 :          WRITE (iw, '(/," ",29("*"),A,29("*"))') " UPF Pseudopotential "
     693           0 :          DO k = 1, potential%upf_pot%maxinfo
     694           0 :             WRITE (iw, '(A80)') potential%upf_pot%info(k)
     695             :          END DO
     696             :       CASE (sgp_pseudo)
     697           0 :          WRITE (iw, '(/," ",29("*"),A,29("*"))') " SGP Pseudopotential "
     698           0 :          WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%sgp_pot%zion
     699             :       CASE (ecp_pseudo)
     700           3 :          WRITE (iw, '(/," ",26("*"),A,27("*"))') " Effective Core Potential "
     701           3 :          WRITE (iw, '(T10,A,T76,F5.1)') " Core Charge ", potential%ecp_pot%zion
     702           6 :          DO k = 1, potential%ecp_pot%nloc
     703           6 :             IF (k == 1) THEN
     704           3 :                WRITE (iw, '(T10,A,T40,I3,T49,2F16.8)') " Local Potential ", potential%ecp_pot%nrloc(k), &
     705           6 :                   potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
     706             :             ELSE
     707           0 :                WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrloc(k), &
     708           0 :                   potential%ecp_pot%bloc(k), potential%ecp_pot%aloc(k)
     709             :             END IF
     710             :          END DO
     711          14 :          DO l = 0, potential%ecp_pot%lmax
     712          11 :             WRITE (iw, '(T10,A,I3)') " ECP l-value ", l
     713          29 :             DO k = 1, potential%ecp_pot%npot(l)
     714          15 :                WRITE (iw, '(T40,I3,T49,2F16.8)') potential%ecp_pot%nrpot(k, l), &
     715          41 :                   potential%ecp_pot%bpot(k, l), potential%ecp_pot%apot(k, l)
     716             :             END DO
     717             :          END DO
     718             :       CASE DEFAULT
     719           3 :          CPABORT("")
     720             :       END SELECT
     721           3 :       IF (potential%confinement) THEN
     722           0 :          IF (potential%conf_type == poly_conf) THEN
     723             :             WRITE (iw, '(/,T10,A,T51,F12.6," * (R /",F6.2,")**",F6.2)') &
     724           0 :                " Confinement Potential ", potential%acon, potential%rcon, potential%scon
     725           0 :          ELSE IF (potential%conf_type == barrier_conf) THEN
     726           0 :             WRITE (iw, '(/,T10,A)') " Confinement Potential s*F[(r-ron)/w] "
     727           0 :             WRITE (iw, '(T57,A,F12.6,A)') "s     =", potential%acon, "   Ha"
     728           0 :             WRITE (iw, '(T57,A,F12.6,A)') "w     =", potential%rcon, " Bohr"
     729           0 :             WRITE (iw, '(T57,A,F12.6,A)') "ron   =", potential%scon, " Bohr"
     730             :          ELSE
     731           0 :             CPABORT("")
     732             :          END IF
     733             :       ELSE
     734           3 :          WRITE (iw, '(/,T10,A)') " No Confinement Potential is applied "
     735             :       END IF
     736           3 :       WRITE (iw, '(" ",79("*"))')
     737             : 
     738           3 :    END SUBROUTINE atom_print_potential
     739             : 
     740             : ! **************************************************************************************************
     741             : !> \brief Print GTH pseudo-potential parameters.
     742             : !> \param gthpot  pseudo-potential
     743             : !> \param iunit   output file unit
     744             : !> \par History
     745             : !>    * 09.2012 created [Juerg Hutter]
     746             : !> \note  The pseudo-potential is written into the 'iunit' file unit or as the file 'GTH-PARAMETER'
     747             : !>        inside the current working directory if the I/O unit is not given explicitly.
     748             : ! **************************************************************************************************
     749          37 :    SUBROUTINE atom_write_pseudo_param(gthpot, iunit)
     750             :       TYPE(atom_gthpot_type), INTENT(INOUT)              :: gthpot
     751             :       INTEGER, INTENT(IN), OPTIONAL                      :: iunit
     752             : 
     753             :       INTEGER                                            :: i, iw, j, k, n
     754             : 
     755          37 :       IF (PRESENT(iunit)) THEN
     756           2 :          iw = iunit
     757             :       ELSE
     758          35 :          CALL open_file(file_name="GTH-PARAMETER", file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
     759             :       END IF
     760          37 :       WRITE (iw, '(A2,A)') gthpot%symbol, ADJUSTL(TRIM(gthpot%pname))
     761          37 :       WRITE (iw, '(4I5)') gthpot%econf(0:3)
     762         111 :       WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rc, gthpot%ncl, (gthpot%cl(i), i=1, gthpot%ncl)
     763          37 :       IF (gthpot%lpotextended) THEN
     764           0 :          WRITE (iw, '(A,I5)') "  LPOT", gthpot%nexp_lpot
     765           0 :          DO i = 1, gthpot%nexp_lpot
     766           0 :             WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lpot(i), gthpot%nct_lpot(i), &
     767           0 :                (gthpot%cval_lpot(j, i), j=1, gthpot%nct_lpot(i))
     768             :          END DO
     769             :       END IF
     770          37 :       IF (gthpot%lsdpot) THEN
     771           0 :          WRITE (iw, '(A,I5)') "  LSD ", gthpot%nexp_lsd
     772           0 :          DO i = 1, gthpot%nexp_lsd
     773           0 :             WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_lsd(i), gthpot%nct_lsd(i), &
     774           0 :                (gthpot%cval_lsd(j, i), j=1, gthpot%nct_lsd(i))
     775             :          END DO
     776             :       END IF
     777          37 :       IF (gthpot%nlcc) THEN
     778          24 :          WRITE (iw, '(A,I5)') " NLCC ", gthpot%nexp_nlcc
     779          48 :          DO i = 1, gthpot%nexp_nlcc
     780          24 :             WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%alpha_nlcc(i), gthpot%nct_nlcc(i), &
     781          96 :                (gthpot%cval_nlcc(j, i)*4.0_dp*pi, j=1, gthpot%nct_nlcc(i))
     782             :          END DO
     783             :       END IF
     784          37 :       n = 0
     785         201 :       DO i = lmat, 0, -1
     786         201 :          IF (gthpot%nl(i) > 0) THEN
     787          34 :             n = i + 1
     788          34 :             EXIT
     789             :          END IF
     790             :       END DO
     791          37 :       WRITE (iw, '(I8)') n
     792          95 :       DO i = 0, n - 1
     793         116 :          WRITE (iw, '(F20.14,I8,5F20.14)') gthpot%rcnl(i), gthpot%nl(i), (gthpot%hnl(1, k, i), k=1, gthpot%nl(i))
     794          37 :          SELECT CASE (gthpot%nl(i))
     795             :          CASE (2)
     796           0 :             WRITE (iw, '(T49,F20.14)') gthpot%hnl(2, 2, i)
     797             :          CASE (3)
     798           0 :             WRITE (iw, '(T49,2F20.14)') gthpot%hnl(2, 2, i), gthpot%hnl(2, 3, i)
     799           0 :             WRITE (iw, '(T69,F20.14)') gthpot%hnl(3, 3, i)
     800             :          CASE DEFAULT
     801          58 :             DO j = 2, gthpot%nl(i)
     802          58 :                WRITE (iw, '(T29,5F20.14)') (gthpot%hnl(j, k, i), k=j, gthpot%nl(i))
     803             :             END DO
     804             :          END SELECT
     805             :       END DO
     806          37 :       IF (.NOT. PRESENT(iunit)) CALL close_file(unit_number=iw)
     807             : 
     808          37 :    END SUBROUTINE atom_write_pseudo_param
     809             : 
     810             : ! **************************************************************************************************
     811             : !> \brief Print atomic orbitals.
     812             : !> \param atom  information about the atomic kind
     813             : !> \param iw    output file unit
     814             : !> \param xmgrace ...
     815             : !> \par History
     816             : !>    * 04.2013 created [Juerg Hutter]
     817             : ! **************************************************************************************************
     818           0 :    SUBROUTINE atom_print_orbitals(atom, iw, xmgrace)
     819             :       TYPE(atom_type), POINTER                           :: atom
     820             :       INTEGER, INTENT(IN)                                :: iw
     821             :       LOGICAL, INTENT(IN), OPTIONAL                      :: xmgrace
     822             : 
     823             :       CHARACTER(LEN=40)                                  :: fnbody
     824             :       INTEGER                                            :: z
     825             :       LOGICAL                                            :: graph
     826             : 
     827           0 :       SELECT CASE (atom%method_type)
     828             :       CASE DEFAULT
     829           0 :          CPABORT("")
     830             :       CASE (do_rks_atom)
     831           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
     832             :       CASE (do_uks_atom)
     833           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
     834           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
     835             :       CASE (do_rhf_atom)
     836           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfn, "", iw)
     837             :       CASE (do_uhf_atom)
     838           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfna, "Alpha", iw)
     839           0 :          CALL atom_print_orbitals_helper(atom, atom%orbitals%wfnb, "Beta", iw)
     840             :       CASE (do_rohf_atom)
     841           0 :          CPABORT("")
     842             :       END SELECT
     843             : 
     844           0 :       graph = .FALSE.
     845           0 :       IF (PRESENT(xmgrace)) graph = xmgrace
     846           0 :       IF (graph .AND. iw > 0) THEN
     847           0 :          z = atom%z
     848           0 :          fnbody = TRIM(ptable(z)%symbol)//"_PPorbital"
     849           0 :          SELECT CASE (atom%method_type)
     850             :          CASE DEFAULT
     851           0 :             CPABORT("")
     852             :          CASE (do_rks_atom)
     853           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
     854             :          CASE (do_uks_atom)
     855           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
     856           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
     857             :          CASE (do_rhf_atom)
     858           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfn, fnbody)
     859             :          CASE (do_uhf_atom)
     860           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfna, TRIM(fnbody)//"alpha")
     861           0 :             CALL atom_orbitals_grace(atom, atom%orbitals%wfnb, TRIM(fnbody)//"beta")
     862             :          CASE (do_rohf_atom)
     863           0 :             CPABORT("")
     864             :          END SELECT
     865             :       END IF
     866             : 
     867           0 :    END SUBROUTINE atom_print_orbitals
     868             : 
     869             : ! **************************************************************************************************
     870             : !> \brief Print atomic orbitals of the given spin.
     871             : !> \param atom         information about the atomic kind
     872             : !> \param wfn          atomic orbitals
     873             : !> \param description  description string
     874             : !> \param iw           output file unit
     875             : !> \par History
     876             : !>    * 04.2013 created [Juerg Hutter]
     877             : ! **************************************************************************************************
     878           0 :    SUBROUTINE atom_print_orbitals_helper(atom, wfn, description, iw)
     879             :       TYPE(atom_type), POINTER                           :: atom
     880             :       REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT)  :: wfn
     881             :       CHARACTER(len=*), INTENT(IN)                       :: description
     882             :       INTEGER, INTENT(IN)                                :: iw
     883             : 
     884             :       INTEGER                                            :: b, l, maxl, nb, nv, v
     885             : 
     886           0 :       WRITE (iw, '(/,A,A,A)') " Atomic orbital expansion coefficients [", description, "]"
     887             : 
     888           0 :       maxl = atom%state%maxl_calc
     889           0 :       DO l = 0, maxl
     890             : 
     891           0 :          nb = atom%basis%nbas(l)
     892           0 :          nv = atom%state%maxn_calc(l)
     893           0 :          IF (nb > 0 .AND. nv > 0) THEN
     894           0 :             nv = MIN(nv, SIZE(wfn, 2))
     895           0 :             DO v = 1, nv
     896           0 :                WRITE (iw, '(/,"    ORBITAL      L = ",I1,"      State = ",I3)') l, v
     897           0 :                DO b = 1, nb
     898           0 :                   WRITE (iw, '("      ",ES23.15)') wfn(b, v, l)
     899             :                END DO
     900             :             END DO
     901             :          END IF
     902             :       END DO
     903           0 :    END SUBROUTINE atom_print_orbitals_helper
     904             : 
     905             : ! **************************************************************************************************
     906             : !> \brief Print atomic orbitals of the given spin.
     907             : !> \param atom         information about the atomic kind
     908             : !> \param wfn          atomic orbitals
     909             : !> \param fnbody       body of file name
     910             : !> \par History
     911             : !>    * 02.2025 created [Juerg Hutter]
     912             : ! **************************************************************************************************
     913           0 :    SUBROUTINE atom_orbitals_grace(atom, wfn, fnbody)
     914             :       TYPE(atom_type), POINTER                           :: atom
     915             :       REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT)  :: wfn
     916             :       CHARACTER(len=*), INTENT(IN)                       :: fnbody
     917             : 
     918             :       CHARACTER(LEN=1), DIMENSION(0:8)                   :: lname
     919             :       CHARACTER(LEN=1), DIMENSION(1:9)                   :: wnum
     920             :       CHARACTER(LEN=40)                                  :: fname, legend
     921             :       INTEGER                                            :: b, i, iw, l, m, maxl, nb, nv, v
     922           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gdata, wfnr
     923             :       REAL(KIND=dp), DIMENSION(4)                        :: world_coord
     924             : 
     925           0 :       lname = ['s', 'p', 'd', 'f', 'g', 'h', 'j', 'k', 'l']
     926           0 :       wnum = ['1', '2', '3', '4', '5', '6', '7', '8', '9']
     927           0 :       m = atom%basis%grid%nr
     928           0 :       maxl = atom%state%maxl_calc
     929           0 :       DO l = 0, maxl
     930           0 :          fname = TRIM(fnbody)//"_"//lname(l)//".agr"
     931           0 :          nb = atom%basis%nbas(l)
     932           0 :          nv = atom%state%maxn_calc(l)
     933           0 :          IF (nb > 0 .AND. nv > 0) THEN
     934           0 :             CALL open_file(file_name=fname, file_status="UNKNOWN", file_action="WRITE", unit_number=iw)
     935           0 :             nv = MIN(nv, SIZE(wfn, 2))
     936           0 :             ALLOCATE (wfnr(m, nv))
     937           0 :             wfnr = 0.0_dp
     938           0 :             DO v = 1, nv
     939           0 :                DO b = 1, nb
     940           0 :                   wfnr(:, v) = wfnr(:, v) + wfn(b, v, l)*atom%basis%bf(:, b, l)
     941             :                END DO
     942             :             END DO
     943           0 :             world_coord(1) = 0.0_dp
     944           0 :             world_coord(2) = MINVAL(wfnr) - 0.5_dp
     945           0 :             world_coord(3) = 15.0_dp
     946           0 :             world_coord(4) = MAXVAL(wfnr) + 0.5_dp
     947             :             !
     948           0 :             CALL xm_write_defaults(iw)
     949           0 :             CALL xm_write_frameport(iw)
     950             :             CALL xm_write_frame(iw, world_coord, &
     951             :                                 title="PP Radial Wavefunction", &
     952             :                                 subtitle=lname(l)//"-Quantum Number", &
     953             :                                 xlabel="Radius [Bohr]", &
     954           0 :                                 ylabel="")
     955           0 :             DO i = 0, nv - 1
     956           0 :                legend = "WFN "//wnum(i + 1)
     957           0 :                CALL xm_graph_info(iw, i, 2.5_dp, legend)
     958             :             END DO
     959           0 :             ALLOCATE (gdata(m, 2))
     960           0 :             gdata(1:m, 1) = atom%basis%grid%rad(1:m)
     961           0 :             DO i = 0, nv - 1
     962           0 :                gdata(1:m, 2) = wfnr(1:m, i + 1)
     963           0 :                CALL xm_graph_data(iw, i, gdata)
     964             :             END DO
     965           0 :             DEALLOCATE (gdata, wfnr)
     966           0 :             CALL close_file(iw)
     967             :          END IF
     968             :       END DO
     969           0 :    END SUBROUTINE atom_orbitals_grace
     970             : 
     971             : END MODULE atom_output

Generated by: LCOV version 1.15