LCOV - code coverage report
Current view: top level - src - qs_chargemol.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 80.7 % 404 326
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 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 Write wfx file, works as interface to chargemol and multiwfn
      10              : !> \note  Works only for HF and KS-type wavefunctions, execution is aborted otherwise
      11              : !> \par History
      12              : !>      06.2020 created [Hossam Elgabarty]
      13              : !>      01.2021 modified [Hossam]
      14              : !> \author Hossam Elgabarty
      15              : ! **************************************************************************************************
      16              : MODULE qs_chargemol
      17              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18              :                                               get_atomic_kind,&
      19              :                                               get_atomic_kind_set
      20              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      21              :                                               gto_basis_set_type
      22              :    USE cell_types,                      ONLY: cell_type,&
      23              :                                               get_cell
      24              :    USE cp_control_types,                ONLY: dft_control_type
      25              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      26              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm
      27              :    USE cp_files,                        ONLY: open_file
      28              :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      29              :                                               cp_fm_get_submatrix
      30              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      31              :                                               cp_logger_get_default_io_unit,&
      32              :                                               cp_logger_type
      33              :    USE cp_output_handling,              ONLY: cp_print_key_generate_filename
      34              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      35              :                                               section_vals_type,&
      36              :                                               section_vals_val_get
      37              :    USE kinds,                           ONLY: default_string_length,&
      38              :                                               dp
      39              :    USE kpoint_types,                    ONLY: kpoint_type
      40              :    USE machine,                         ONLY: m_datum
      41              :    USE orbital_pointers,                ONLY: nco,&
      42              :                                               nso
      43              :    USE orbital_symbols,                 ONLY: cgf_symbol,&
      44              :                                               sgf_symbol
      45              :    USE orbital_transformation_matrices, ONLY: orbtramat
      46              :    USE particle_types,                  ONLY: particle_type
      47              :    USE periodic_table,                  ONLY: get_ptable_info
      48              :    USE qs_energy_types,                 ONLY: qs_energy_type
      49              :    USE qs_environment_types,            ONLY: get_qs_env,&
      50              :                                               qs_environment_type
      51              :    USE qs_force_types,                  ONLY: qs_force_type
      52              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      53              :                                               get_qs_kind_set,&
      54              :                                               qs_kind_type
      55              :    USE qs_mo_types,                     ONLY: mo_set_type
      56              :    USE qs_scf_types,                    ONLY: qs_scf_env_type
      57              :    USE qs_subsys_types,                 ONLY: qs_subsys_type
      58              :    USE scf_control_types,               ONLY: scf_control_type
      59              : #include "./base/base_uses.f90"
      60              : 
      61              :    IMPLICIT NONE
      62              :    PRIVATE
      63              : 
      64              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_chargemol'
      65              :    LOGICAL, PARAMETER                   :: debug_this_module = .FALSE.
      66              : 
      67              :    INTEGER, PARAMETER                   :: chargemol_lmax = 5
      68              : 
      69              :    PUBLIC :: write_wfx
      70              : 
      71              : CONTAINS
      72              : 
      73              : ! **************************************************************************************************
      74              : !> \brief ...
      75              : !> \param qs_env ...
      76              : !> \param dft_section ...
      77              : ! **************************************************************************************************
      78            2 :    SUBROUTINE write_wfx(qs_env, dft_section)
      79              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      80              :       TYPE(section_vals_type), POINTER                   :: dft_section
      81              : 
      82              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_wfx'
      83              : 
      84            2 :       CHARACTER(LEN=12), DIMENSION(:), POINTER           :: bcgf_symbol
      85              :       CHARACTER(len=2)                                   :: asym
      86            2 :       CHARACTER(len=20), ALLOCATABLE, DIMENSION(:)       :: atom_symbols
      87              :       CHARACTER(LEN=256)                                 :: datx
      88            2 :       CHARACTER(LEN=6), DIMENSION(:), POINTER            :: bsgf_symbol
      89              :       CHARACTER(len=default_string_length)               :: filename
      90              :       INTEGER :: handle, i, iatom, icgf, ico, ikind, iloop, imo, ipgf, irow, iset, isgf, ishell, &
      91              :          ispin, iwfx, lshell, max_l, ncgf, ncol_global, nelec_alpha, nelec_beta, nkpoint, noccup, &
      92              :          nrow_global, nset, nsgf, nspins, num_atoms, output_unit, roks_high, roks_low, tot_npgf
      93            2 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomic_numbers, kind_of
      94            2 :       INTEGER, DIMENSION(:), POINTER                     :: lmax, npgf, nshell
      95            2 :       INTEGER, DIMENSION(:, :), POINTER                  :: l
      96              :       LOGICAL                                            :: do_kpoints, newl, periodic, unrestricted
      97              :       REAL(KIND=dp)                                      :: zeta
      98            2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atoms_cart, cmatrix, smatrix
      99            2 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: core_charges
     100            2 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
     101            2 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
     102            2 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     103              :       TYPE(cell_type), POINTER                           :: cell
     104              :       TYPE(cp_logger_type), POINTER                      :: logger
     105            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     106              :       TYPE(dft_control_type), POINTER                    :: dft_control
     107              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     108              :       TYPE(kpoint_type), POINTER                         :: kpoint_env
     109            2 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     110            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     111              :       TYPE(qs_energy_type), POINTER                      :: energy
     112            2 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     113            2 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     114              :       TYPE(qs_scf_env_type), POINTER                     :: scf_env
     115              :       TYPE(qs_subsys_type), POINTER                      :: subsys
     116              :       TYPE(scf_control_type), POINTER                    :: scf_control
     117              :       TYPE(section_vals_type), POINTER                   :: chargemol_section
     118              : 
     119              :       ! !--------------------------------------------------------------------------------------------!
     120              : 
     121            2 :       CALL timeset(routineN, handle)
     122              : 
     123            2 :       logger => cp_get_default_logger()
     124            2 :       output_unit = cp_logger_get_default_io_unit(logger)
     125              : 
     126            2 :       IF (output_unit > 0) THEN
     127              :          WRITE (output_unit, '(/,T2,A)') &
     128            1 :             '!-----------------------------------------------------------------------------!'
     129            1 :          WRITE (output_unit, '(T32,A)') "Writing Chargemol wfx file..."
     130              :          WRITE (output_unit, '(T2,A)') &
     131            1 :             '!-----------------------------------------------------------------------------!'
     132              :       END IF
     133              : 
     134              :       ! Collect environment info
     135            2 :       CPASSERT(ASSOCIATED(qs_env))
     136              :       CALL get_qs_env(qs_env, cell=cell, mos=mos, particle_set=particle_set, &
     137              :                       qs_kind_set=qs_kind_set, scf_env=scf_env, scf_control=scf_control, &
     138              :                       dft_control=dft_control, energy=energy, force=force, subsys=subsys, &
     139              :                       atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
     140            2 :                       kpoints=kpoint_env, matrix_s=matrix_s)
     141              : 
     142            2 :       nkpoint = kpoint_env%nkp
     143              : 
     144            2 :       IF (scf_control%use_ot) THEN
     145            0 :          CPABORT("OT is incompatible with .wfx output. Switch off OT.")
     146              :       END IF
     147              : 
     148            2 :       IF (dft_control%roks) THEN
     149              :          CALL cp_abort(__LOCATION__, "The output of chargemol .wfx files is currently "// &
     150            0 :                        "implemented only for (un)restricted HF and Kohn-Sham-like methods.")
     151              :       END IF
     152              : 
     153            2 :       IF (dft_control%lsd) THEN
     154              :          unrestricted = .TRUE.
     155              :       ELSE
     156            2 :          unrestricted = .FALSE.
     157              :       END IF
     158              : 
     159            2 :       chargemol_section => section_vals_get_subs_vals(dft_section, "PRINT%CHARGEMOL")
     160              : 
     161            2 :       CALL section_vals_val_get(chargemol_section, "PERIODIC", l_val=periodic)
     162              : 
     163              :     !!---------------------------------------------------------------------------------!
     164              :       ! output unit
     165              :     !!---------------------------------------------------------------------------------!
     166              : 
     167              :       filename = cp_print_key_generate_filename(logger, chargemol_section, &
     168            2 :                                                 extension=".wfx", my_local=.FALSE.)
     169              : 
     170              :       CALL open_file(filename, file_status="UNKNOWN", &
     171              :                      file_action="WRITE", &
     172            2 :                      unit_number=iwfx)
     173              : 
     174              :     !!---------------------------------------------------------------------------------!
     175              : 
     176            2 :       IF (iwfx > 0) THEN
     177              : 
     178            2 :          nspins = SIZE(mos)
     179            2 :          IF (nspins == 1) THEN
     180            2 :             nelec_alpha = mos(1)%nelectron/2
     181            2 :             nelec_beta = nelec_alpha
     182              :          ELSE
     183            0 :             nelec_alpha = mos(1)%nelectron
     184            0 :             nelec_beta = mos(2)%nelectron
     185              :          END IF
     186              : 
     187            2 :          IF (dft_control%roks) THEN
     188            0 :             IF (mos(1)%homo > mos(2)%homo) THEN
     189              :                roks_high = 1
     190              :                roks_low = 2
     191              :             ELSE
     192            0 :                roks_high = 2
     193            0 :                roks_low = 1
     194              :             END IF
     195              :          END IF
     196              : 
     197              :        !!---------------------------------------------------------------------------------!
     198              :          ! Initial parsing of topology and basis set, check maximum l .le. chargemol_lmax
     199              :        !!---------------------------------------------------------------------------------!
     200            2 :          num_atoms = SIZE(particle_set)
     201            6 :          ALLOCATE (atoms_cart(3, num_atoms))
     202            6 :          ALLOCATE (atom_symbols(num_atoms))
     203            6 :          ALLOCATE (atomic_numbers(num_atoms))
     204            6 :          ALLOCATE (core_charges(num_atoms))
     205              : 
     206            2 :          max_l = 0
     207            4 :          DO iatom = 1, num_atoms
     208            2 :             NULLIFY (orb_basis_set)
     209            8 :             atoms_cart(1:3, iatom) = particle_set(iatom)%r(1:3)
     210              :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
     211            2 :                                  element_symbol=asym, kind_number=ikind)
     212            2 :             atom_symbols(iatom) = asym
     213            2 :             CALL get_ptable_info(asym, number=atomic_numbers(iatom))
     214              :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
     215            2 :                              core_charge=core_charges(iatom))
     216            2 :             IF (ASSOCIATED(orb_basis_set)) THEN
     217              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     218            2 :                                       nset=nset, lmax=lmax)
     219              : 
     220            6 :                DO iset = 1, nset
     221            6 :                   max_l = MAX(max_l, lmax(iset))
     222              :                END DO
     223              :             ELSE
     224            0 :                CPABORT("Unknown basis set type. ")
     225              :             END IF
     226            6 :             IF (max_l > chargemol_lmax) THEN
     227              :                CALL cp_abort(__LOCATION__, "Chargemol supports basis sets with l"// &
     228            0 :                              " up to 5 only (h angular functions).")
     229              :             END IF
     230              :          END DO
     231              : 
     232              :          ! !===========================================================================!
     233              :          ! Header comment and Title
     234              :          ! !===========================================================================!
     235            2 :          WRITE (iwfx, "(A)") "# Chargemol input file generated by CP2K "
     236            2 :          CALL m_datum(datx)
     237            2 :          WRITE (iwfx, "(A,/)") "# Creation date "//TRIM(datx)
     238              : 
     239            2 :          WRITE (iwfx, "(A)") "<Title>"
     240            2 :          WRITE (iwfx, *) TRIM(logger%iter_info%project_name)
     241            2 :          WRITE (iwfx, "(A,/)") "</Title>"
     242              : 
     243              :          ! !===========================================================================!
     244              :          ! Keywords
     245              :          ! TODO: check: always GTO??
     246              :          ! !===========================================================================!
     247            2 :          WRITE (iwfx, "(A)") "<Keywords>"
     248            2 :          WRITE (iwfx, "(A)") "GTO"
     249            2 :          WRITE (iwfx, "(A,/)") "</Keywords>"
     250              : 
     251              :          ! !===========================================================================!
     252              :          ! Model
     253              :          ! !===========================================================================!
     254            2 :          WRITE (iwfx, "(A)") "#<Model>"
     255            2 :          IF (unrestricted) THEN
     256            0 :             WRITE (iwfx, "(A)") "#Unrestricted Kohn-Sham"
     257              :          ELSE
     258            2 :             WRITE (iwfx, "(A)") "#Restricted Kohn-Sham"
     259              :          END IF
     260            2 :          WRITE (iwfx, "(A,/)") "#</Model>"
     261              : 
     262              :          ! !===========================================================================!
     263              :          ! Number of nuclei
     264              :          ! !===========================================================================!
     265            2 :          WRITE (iwfx, "(A)") "<Number of Nuclei>"
     266            2 :          WRITE (iwfx, "(I6)") num_atoms
     267            2 :          WRITE (iwfx, "(A,/)") "</Number of Nuclei>"
     268              : 
     269              :          ! !===========================================================================!
     270              :          ! Number of Occupied MOs
     271              :          ! !===========================================================================!
     272            2 :          WRITE (iwfx, "(A)") "<Number of Occupied Molecular Orbitals>"
     273            2 :          noccup = 0
     274            2 :          IF (dft_control%roks .AND. nspins == 2) THEN
     275            0 :             noccup = MAX(mos(1)%homo, mos(2)%homo)
     276              :          ELSE
     277            4 :             DO ispin = 1, dft_control%nspins
     278            4 :                noccup = noccup + mos(ispin)%homo
     279              :             END DO
     280              :          END IF
     281            2 :          WRITE (iwfx, "(2I6)") noccup
     282            2 :          WRITE (iwfx, "(A,/)") "</Number of Occupied Molecular Orbitals>"
     283              : 
     284              :          ! !===========================================================================!
     285              :          ! Number of Perturbations
     286              :          ! !===========================================================================!
     287            2 :          WRITE (iwfx, "(A)") "<Number of Perturbations>"
     288            2 :          WRITE (iwfx, "(I6)") 0
     289            2 :          WRITE (iwfx, "(A,/)") "</Number of Perturbations>"
     290              : 
     291              :          ! !===========================================================================!
     292              :          ! Total charge
     293              :          ! !===========================================================================!
     294            2 :          WRITE (iwfx, "(A)") "<Net Charge>"
     295            2 :          WRITE (iwfx, "(F8.4)") REAL(dft_control%charge, KIND=dp)
     296            2 :          WRITE (iwfx, "(A,/)") "</Net Charge>"
     297              : 
     298              :          ! !===========================================================================!
     299              :          ! Number of electrons
     300              :          ! !===========================================================================!
     301            2 :          WRITE (iwfx, "(A)") "<Number of Electrons>"
     302            2 :          WRITE (iwfx, "(I6)") scf_env%nelectron
     303            2 :          WRITE (iwfx, "(A,/)") "</Number of Electrons>"
     304              : 
     305              :          !===========================================================================!
     306              :          ! Number of Alpha electrons
     307              :          !===========================================================================!
     308            2 :          WRITE (iwfx, "(A)") "<Number of Alpha Electrons>"
     309            2 :          WRITE (iwfx, "(I6)") nelec_alpha
     310            2 :          WRITE (iwfx, "(A,/)") "</Number of Alpha Electrons>"
     311              : 
     312              :          !===========================================================================!
     313              :          ! Number of Beta electrons
     314              :          !===========================================================================!
     315            2 :          WRITE (iwfx, "(A)") "<Number of Beta Electrons>"
     316            2 :          WRITE (iwfx, "(I6)") nelec_beta
     317            2 :          WRITE (iwfx, "(A,/)") "</Number of Beta Electrons>"
     318              : 
     319              :          !===========================================================================!
     320              :          ! Spin multiplicity
     321              :          !===========================================================================!
     322            2 :          WRITE (iwfx, "(A)") "<Electron Spin Multiplicity>"
     323            2 :          WRITE (iwfx, "(I4)") dft_control%multiplicity
     324            2 :          WRITE (iwfx, "(A,/)") "</Electron Spin Multiplicity>"
     325              : 
     326              :          ! !===========================================================================!
     327              :          ! Number of Core Electrons
     328              :          ! !===========================================================================!
     329            2 :          WRITE (iwfx, "(A)") "<Number of Core Electrons>"
     330            6 :          WRITE (iwfx, "(F16.10)") SUM(atomic_numbers) - SUM(core_charges)
     331            2 :          WRITE (iwfx, "(A,/)") "</Number of Core Electrons>"
     332              : 
     333              :          ! !===========================================================================!
     334              :          ! Nuclear Names
     335              :          ! !===========================================================================!
     336            2 :          WRITE (iwfx, "(A)") "<Nuclear Names>"
     337            4 :          DO iatom = 1, num_atoms
     338            4 :             WRITE (iwfx, "(A4)") atom_symbols(iatom)
     339              :          END DO
     340            2 :          WRITE (iwfx, "(A,/)") "</Nuclear Names>"
     341              : 
     342              :          ! !===========================================================================!
     343              :          ! Atomic Numbers
     344              :          ! !===========================================================================!
     345            2 :          WRITE (iwfx, "(A)") "<Atomic Numbers>"
     346            4 :          DO iatom = 1, num_atoms
     347            4 :             WRITE (iwfx, "(I4)") atomic_numbers(iatom)
     348              :          END DO
     349            2 :          WRITE (iwfx, "(A,/)") "</Atomic Numbers>"
     350              : 
     351              :          ! !===========================================================================!
     352              :          ! Nuclear charges
     353              :          ! !===========================================================================!
     354            2 :          WRITE (iwfx, "(A)") "<Nuclear Charges>"
     355            4 :          DO iatom = 1, num_atoms
     356            4 :             WRITE (iwfx, "(F12.8)") core_charges(iatom)
     357              :          END DO
     358            2 :          WRITE (iwfx, "(A,/)") "</Nuclear Charges>"
     359              : 
     360              :          ! !===========================================================================!
     361              :          ! Nuclear coordinates
     362              :          ! !===========================================================================!
     363            2 :          WRITE (iwfx, "(A)") "<Nuclear Cartesian Coordinates>"
     364            4 :          DO iatom = 1, num_atoms
     365            4 :             WRITE (iwfx, "(3ES26.16E3)") atoms_cart(1:3, iatom)
     366              :          END DO
     367            2 :          WRITE (iwfx, "(A,/)") "</Nuclear Cartesian Coordinates>"
     368              : 
     369              :          ! !===========================================================================!
     370              :          ! Periodic boundary conditions
     371              :          ! !===========================================================================!
     372            2 :          IF (periodic) THEN
     373            2 :             CALL get_cell(cell)
     374            8 :             IF (SUM(cell%perd) == 0) THEN
     375              :                CALL cp_warn(__LOCATION__, "Writing of periodic cell information"// &
     376            0 :                             " requested in input, but system is not periodic, ignoring request.")
     377              :             ELSE
     378            2 :                WRITE (iwfx, "(A)") "<Number of Translation Vectors>"
     379            8 :                WRITE (iwfx, "(I3)") SUM(cell%perd)
     380            2 :                WRITE (iwfx, "(A,/)") "</Number of Translation Vectors>"
     381              : 
     382            2 :                WRITE (iwfx, "(A)") "<Translation Vectors>"
     383            8 :                DO iatom = 1, 3
     384            8 :                   IF (cell%perd(iatom) == 1) THEN
     385            6 :                      WRITE (iwfx, "(3F12.6)") cell%hmat(1:3, iatom)
     386              :                   END IF
     387              :                END DO
     388            2 :                WRITE (iwfx, "(A,/)") "</Translation Vectors>"
     389              :             END IF
     390            2 :             WRITE (iwfx, "(A)") "<Number of Kpoints>"
     391            2 :             WRITE (iwfx, "(I3)") 1
     392            2 :             WRITE (iwfx, "(A,/)") "</Number of Kpoints>"
     393            2 :             WRITE (iwfx, "(A)") "<Kpoint Weights>"
     394            2 :             WRITE (iwfx, "(I3)") 1
     395            2 :             WRITE (iwfx, "(A,/)") "</Kpoint Weights>"
     396            2 :             WRITE (iwfx, "(A)") "<Kpoint Fractional Coordinates>"
     397            2 :             WRITE (iwfx, "(F8.6)") 0.0
     398            2 :             WRITE (iwfx, "(A,/)") "</Kpoint Fractional Coordinates>"
     399              :          END IF
     400              : 
     401              :          ! !===========================================================================!
     402              :          ! Primitive centers
     403              :          ! !===========================================================================!
     404            2 :          WRITE (iwfx, "(A)") "<Primitive Centers>"
     405            2 :          tot_npgf = 0
     406            4 :          DO iatom = 1, num_atoms
     407            2 :             iloop = 1
     408            2 :             NULLIFY (orb_basis_set)
     409            2 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     410            2 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     411            6 :             IF (ASSOCIATED(orb_basis_set)) THEN
     412              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     413              :                                       nset=nset, &
     414              :                                       npgf=npgf, &
     415              :                                       nshell=nshell, &
     416              :                                       l=l, &
     417            2 :                                       gcc=gcc)
     418              : 
     419            6 :                DO iset = 1, nset
     420           16 :                   DO ishell = 1, nshell(iset)
     421           10 :                      lshell = l(ishell, iset)
     422           42 :                      DO ico = 1, nco(lshell)
     423          114 :                         DO ipgf = 1, npgf(iset)
     424           76 :                            tot_npgf = tot_npgf + 1
     425           76 :                            IF (MOD(iloop + 10, 10) /= 0) THEN
     426           70 :                               WRITE (iwfx, "(I4)", advance="no") iatom
     427           70 :                               newl = .TRUE.
     428              :                            ELSE
     429            6 :                               WRITE (iwfx, "(I4)") iatom
     430            6 :                               newl = .FALSE.
     431              :                            END IF
     432          104 :                            iloop = iloop + 1
     433              :                            !END IF
     434              :                         END DO
     435              :                      END DO
     436              :                   END DO
     437              :                END DO
     438            2 :                IF (newl) THEN
     439            2 :                   WRITE (iwfx, *)
     440              :                END IF
     441              :             ELSE
     442            0 :                CPABORT("Unknown basis set type. ")
     443              :             END IF
     444              :          END DO
     445            2 :          WRITE (iwfx, "(A,/)") "</Primitive Centers>"
     446              : 
     447              :          ! !===========================================================================!
     448              :          ! Number of primitives
     449              :          ! !===========================================================================!
     450            2 :          WRITE (iwfx, "(A)") "<Number of Primitives>"
     451            2 :          WRITE (iwfx, "(I8)") tot_npgf
     452            2 :          WRITE (iwfx, "(A,/)") "</Number of Primitives>"
     453              : 
     454              :          ! !===========================================================================!
     455              :          ! Primitive Types
     456              :          ! !===========================================================================!
     457            2 :          WRITE (iwfx, "(A)") "<Primitive Types>"
     458            4 :          DO iatom = 1, num_atoms
     459            2 :             iloop = 1
     460            2 :             NULLIFY (orb_basis_set)
     461            2 :             NULLIFY (bcgf_symbol)
     462            2 :             NULLIFY (gcc)
     463            2 :             NULLIFY (zet)
     464            2 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     465            2 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     466            2 :             IF (ASSOCIATED(orb_basis_set)) THEN
     467              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     468              :                                       nset=nset, &
     469              :                                       nshell=nshell, &
     470              :                                       l=l, &
     471              :                                       npgf=npgf, &
     472              :                                       ncgf=ncgf, &
     473              :                                       zet=zet, &
     474              :                                       cgf_symbol=bcgf_symbol, &
     475              :                                       sgf_symbol=bsgf_symbol, &
     476            2 :                                       gcc=gcc)
     477              : 
     478            2 :                icgf = 1
     479            6 :                DO iset = 1, nset
     480           16 :                   DO ishell = 1, nshell(iset)
     481           10 :                      lshell = l(ishell, iset)
     482           42 :                      DO ico = 1, nco(lshell)
     483          104 :                         DO ipgf = 1, orb_basis_set%npgf(iset)
     484           76 :                            IF (MOD(iloop + 10, 10) /= 0) THEN
     485              :                               WRITE (iwfx, '(I6)', advance="no") &
     486           70 :                                  pgf_type_chargemol(bcgf_symbol(icgf) (3:))
     487           70 :                               newl = .TRUE.
     488              :                            ELSE
     489              :                               WRITE (iwfx, '(I6)') &
     490            6 :                                  pgf_type_chargemol(bcgf_symbol(icgf) (3:))
     491            6 :                               newl = .FALSE.
     492              :                            END IF
     493          104 :                            iloop = iloop + 1
     494              :                            !END IF
     495              :                         END DO
     496           38 :                         icgf = icgf + 1
     497              :                      END DO !ico
     498              :                   END DO ! ishell
     499              :                END DO ! iset
     500              : 
     501              :             ELSE
     502            0 :                CPABORT("Unknown basis set type.")
     503              :             END IF
     504            6 :             IF (newl) THEN
     505            2 :                WRITE (iwfx, *)
     506              :             END IF
     507              :          END DO ! iatom
     508            2 :          WRITE (iwfx, "(A,/)") "</Primitive Types>"
     509              : 
     510              :          ! !===========================================================================!
     511              :          ! Primitive Exponents
     512              :          ! !===========================================================================!
     513              :        !! NOTE: CP2K orders the atomic orbitals as follows:
     514              :        !! Cartesian orbitals: x (slowest loop), y, z (fastest loop)
     515              :        !! Spherical orbitals: m = -l, ..., 0, ..., +l
     516              : 
     517            2 :          WRITE (iwfx, "(A)") "<Primitive Exponents>"
     518            4 :          DO iatom = 1, num_atoms
     519            2 :             NULLIFY (orb_basis_set)
     520            2 :             NULLIFY (gcc)
     521            2 :             NULLIFY (zet)
     522            2 :             iloop = 1
     523              : 
     524            2 :             CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     525            2 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
     526            6 :             IF (ASSOCIATED(orb_basis_set)) THEN
     527              :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     528              :                                       nset=nset, &
     529              :                                       npgf=npgf, &
     530              :                                       nshell=nshell, &
     531              :                                       l=l, &
     532              :                                       zet=zet, &
     533            2 :                                       gcc=gcc)
     534              : 
     535            6 :                DO iset = 1, nset
     536           16 :                   DO ishell = 1, nshell(iset)
     537           10 :                      lshell = l(ishell, iset)
     538           42 :                      DO ico = 1, nco(lshell)
     539          114 :                         DO ipgf = 1, npgf(iset)
     540           76 :                            IF (MOD(iloop + 5, 5) /= 0) THEN
     541           62 :                               WRITE (iwfx, "(ES20.10)", advance="no") zet(ipgf, iset)
     542           62 :                               newl = .TRUE.
     543              :                            ELSE
     544           14 :                               WRITE (iwfx, "(ES20.10)") zet(ipgf, iset)
     545           14 :                               newl = .FALSE.
     546              :                            END IF
     547          104 :                            iloop = iloop + 1
     548              :                            !END IF
     549              :                         END DO
     550              :                      END DO
     551              :                   END DO
     552              :                END DO
     553            2 :                IF (newl) THEN
     554            2 :                   WRITE (iwfx, *)
     555              :                END IF
     556              :             ELSE
     557            0 :                CPABORT("Unknown basis set type. ")
     558              :             END IF
     559              :          END DO
     560            2 :          WRITE (iwfx, "(A,/)") "</Primitive Exponents>"
     561              : 
     562              :          ! !===========================================================================!
     563              :          !  MO occupation numbers
     564              :          !  nonesense for roks at the moment
     565              :          ! !===========================================================================!
     566            2 :          WRITE (iwfx, "(A)") "<Molecular Orbital Occupation Numbers>"
     567            2 :          IF (.NOT. dft_control%roks) THEN
     568            4 :             DO ispin = 1, nspins
     569           12 :                DO imo = 1, mos(ispin)%homo
     570              :                   WRITE (iwfx, "(f8.6)") &
     571           10 :                      mos(ispin)%occupation_numbers(imo)
     572              :                END DO
     573              :             END DO
     574              :          ELSE
     575            0 :             DO imo = 1, mos(roks_low)%homo
     576              :                WRITE (iwfx, "(*(f8.6,1x))") &
     577              :                   mos(roks_low)%occupation_numbers(imo) &
     578            0 :                   + mos(roks_low)%occupation_numbers(imo)
     579              :             END DO
     580            0 :             DO imo = mos(roks_low)%homo + 1, mos(roks_high)%homo
     581              :                WRITE (iwfx, "(f8.6)") &
     582            0 :                   mos(roks_high)%occupation_numbers(imo)
     583              :             END DO
     584              :          END IF
     585            2 :          WRITE (iwfx, "(A,/)") "</Molecular Orbital Occupation Numbers>"
     586              : 
     587              :          ! !===========================================================================!
     588              :          ! MO energies
     589              :          ! !===========================================================================!
     590            2 :          WRITE (iwfx, "(A)") "<Molecular Orbital Energies>"
     591            4 :          DO ispin = 1, nspins
     592           12 :             DO imo = 1, mos(ispin)%homo
     593           10 :                WRITE (iwfx, "(ES20.10)") mos(ispin)%eigenvalues(imo)
     594              :             END DO
     595              :          END DO
     596            2 :          WRITE (iwfx, "(A,/)") "</Molecular Orbital Energies>"
     597              : 
     598              :          ! !===========================================================================!
     599              :          ! MO Spin types
     600              :          ! nonesense for ROKS
     601              :          ! !===========================================================================!
     602            2 :          WRITE (iwfx, "(A)") "<Molecular Orbital Spin Types>"
     603           10 :          DO imo = 1, mos(1)%homo
     604           10 :             IF (nspins == 1) THEN
     605            8 :                WRITE (iwfx, '(A15)') "Alpha and Beta"
     606            0 :             ELSEIF (dft_control%uks) THEN
     607            0 :                WRITE (iwfx, '(A6)') "Alpha"
     608              :             ELSE
     609              :                CALL cp_abort(__LOCATION__, "This wavefunction type is currently"// &
     610            0 :                              " not supported by the chargemol interface.")
     611              :             END IF
     612              :          END DO
     613            2 :          IF (nspins == 2) THEN
     614            0 :             DO imo = 1, mos(2)%homo
     615            0 :                WRITE (iwfx, '(A5)') "Beta"
     616              :             END DO
     617              :          END IF
     618            2 :          WRITE (iwfx, "(A,/)") "</Molecular Orbital Spin Types>"
     619              : 
     620              :          ! !===========================================================================!
     621              :          ! Kpoint index of orbitals
     622              :          ! !===========================================================================!
     623            2 :          IF (periodic) THEN
     624            2 :             WRITE (iwfx, "(A)") "<Kpoint Index for Orbitals>"
     625            4 :             DO ispin = 1, nspins
     626           12 :                DO imo = 1, mos(ispin)%homo
     627           10 :                   WRITE (iwfx, "(I6)") 1
     628              :                END DO
     629              :             END DO
     630            2 :             WRITE (iwfx, "(A,/)") "</Kpoint Index for Orbitals>"
     631              :          END IF
     632              : 
     633              :          ! !===========================================================================!
     634              :          ! MO Primitive Coefficients
     635              :          ! !===========================================================================!
     636            2 :          WRITE (iwfx, "(A)") "<Molecular Orbital Primitive Coefficients>"
     637            2 :          NULLIFY (orb_basis_set)
     638            2 :          NULLIFY (gcc)
     639              : 
     640            2 :          IF (mos(1)%use_mo_coeff_b) THEN
     641            0 :             DO ispin = 1, SIZE(mos)
     642            0 :                IF (.NOT. ASSOCIATED(mos(ispin)%mo_coeff_b)) THEN
     643            0 :                   CPASSERT(.FALSE.)
     644              :                END IF
     645              :                CALL copy_dbcsr_to_fm(mos(ispin)%mo_coeff_b, &
     646            0 :                                      mos(ispin)%mo_coeff)
     647              :             END DO
     648              :          END IF
     649              : 
     650            4 :          DO ispin = 1, nspins
     651              :             CALL cp_fm_get_info(mos(ispin)%mo_coeff, &
     652              :                                 nrow_global=nrow_global, &
     653            2 :                                 ncol_global=ncol_global)
     654              : 
     655            8 :             ALLOCATE (smatrix(nrow_global, ncol_global))
     656          114 :             smatrix = 0.0_dp
     657              : 
     658            2 :             CALL cp_fm_get_submatrix(mos(ispin)%mo_coeff, smatrix)
     659              : 
     660            2 :             CALL get_qs_kind_set(qs_kind_set, ncgf=ncgf, nsgf=nsgf)
     661              : 
     662            8 :             ALLOCATE (cmatrix(ncgf, ncol_global))
     663              : 
     664          122 :             cmatrix = 0.0_dp
     665              : 
     666              :             ! get MO coefficients in terms of contracted cartesian functions
     667            2 :             icgf = 1
     668            2 :             isgf = 1
     669            4 :             DO iatom = 1, num_atoms
     670            2 :                NULLIFY (orb_basis_set)
     671            2 :                CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     672              :                CALL get_qs_kind(qs_kind_set(ikind), &
     673            2 :                                 basis_set=orb_basis_set)
     674            6 :                IF (ASSOCIATED(orb_basis_set)) THEN
     675              :                   CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     676              :                                          nset=nset, &
     677              :                                          nshell=nshell, &
     678            2 :                                          l=l)
     679            6 :                   DO iset = 1, nset
     680           16 :                      DO ishell = 1, nshell(iset)
     681           10 :                         lshell = l(ishell, iset)
     682              :                         CALL dgemm("T", "N", nco(lshell), mos(ispin)%nmo, &
     683              :                                    nso(lshell), 1.0_dp, &
     684              :                                    orbtramat(lshell)%c2s, nso(lshell), &
     685              :                                    smatrix(isgf, 1), nsgf, 0.0_dp, &
     686           10 :                                    cmatrix(icgf, 1), ncgf)
     687              :                         !IF (iatom==1 .and. debug_this_module) THEN
     688              :                         !   print *, lshell
     689              :                         !   print *, "size ", size(orbtramat(lshell)%s2c,1),",",size(orbtramat(lshell)%s2c,2)
     690              :                         !   DO itest =  1, size(orbtramat(lshell)%s2c,1)
     691              :                         !      print *, orbtramat(lshell)%s2c(itest,:)
     692              :                         !   END DO
     693              :                         !   print *, " "
     694              :                         !   DO itest =  1, 5
     695              :                         !      print *, my_d_orbtramat(itest,:)
     696              :                         !   END DO
     697              :                         !END IF
     698           10 :                         icgf = icgf + nco(lshell)
     699           14 :                         isgf = isgf + nso(lshell)
     700              :                      END DO
     701              :                   END DO
     702              :                ELSE
     703            0 :                   CPABORT("Unknown basis set type. ")
     704              :                END IF
     705              :             END DO ! iatom
     706              : 
     707              :             ! Now write MO coefficients in terms of cartesian primitives
     708           10 :             DO imo = 1, mos(ispin)%homo
     709              : 
     710            8 :                WRITE (iwfx, "(A)") "<MO Number>"
     711            8 :                IF (nspins == 1 .OR. ispin == 1) THEN
     712            8 :                   WRITE (iwfx, "(I6)") imo
     713              :                ELSE
     714            0 :                   WRITE (iwfx, "(I6)") imo + mos(1)%homo
     715              :                END IF
     716            8 :                WRITE (iwfx, "(A)") "</MO Number>"
     717              : 
     718            8 :                irow = 1    ! row of cmatrix, each row corresponds to
     719              :                ! a contracted Cartesian function
     720           18 :                DO iatom = 1, num_atoms
     721            8 :                   iloop = 1
     722            8 :                   NULLIFY (orb_basis_set)
     723              :                   CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
     724            8 :                                        kind_number=ikind)
     725              :                   CALL get_qs_kind(qs_kind_set(ikind), &
     726            8 :                                    basis_set=orb_basis_set)
     727            8 :                   IF (ASSOCIATED(orb_basis_set)) THEN
     728              :                      CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     729              :                                             nset=nset, &
     730              :                                             nshell=nshell, &
     731              :                                             gcc=gcc, &
     732            8 :                                             l=l)
     733              : 
     734            8 :                      icgf = 1
     735           24 :                      DO iset = 1, nset
     736           64 :                         DO ishell = 1, nshell(iset)
     737           40 :                            lshell = orb_basis_set%l(ishell, iset)
     738              : 
     739          152 :                            DO ico = orb_basis_set%first_cgf(ishell, iset), &
     740           56 :                               orb_basis_set%last_cgf(ishell, iset)
     741              : 
     742              :                               !loop over primitives
     743          416 :                               DO ipgf = 1, orb_basis_set%npgf(iset)
     744              : 
     745          304 :                                  zeta = orb_basis_set%zet(ipgf, iset)
     746              : 
     747          304 :                                  IF (MOD(iloop + 5, 5) /= 0) THEN
     748              :                                     WRITE (iwfx, '(ES20.10)', advance="no") &
     749              :                                        cmatrix(irow, imo)* &
     750              :                                        orb_basis_set%norm_cgf(ico)* &
     751          248 :                                        orb_basis_set%gcc(ipgf, ishell, iset)
     752          248 :                                     newl = .TRUE.
     753              :                                  ELSE
     754              :                                     WRITE (iwfx, '(ES20.10)') &
     755              :                                        cmatrix(irow, imo)* &
     756              :                                        orb_basis_set%norm_cgf(ico)* &
     757           56 :                                        orb_basis_set%gcc(ipgf, ishell, iset)
     758           56 :                                     newl = .FALSE.
     759              :                                  END IF
     760          416 :                                  iloop = iloop + 1
     761              :                               END DO ! end loop over primitives
     762          152 :                               irow = irow + 1
     763              :                            END DO !ico
     764              :                         END DO ! ishell
     765              :                      END DO ! iset
     766              : 
     767              :                   ELSE
     768            0 :                      CPABORT("Unknown basis set type.")
     769              :                   END IF
     770           24 :                   IF (newl) THEN
     771            8 :                      WRITE (iwfx, *)
     772              :                   END IF
     773              :                END DO ! iatom
     774              :                !Write (iwfx,*)
     775              :             END DO ! imo
     776              : 
     777            2 :             IF (ALLOCATED(cmatrix)) DEALLOCATE (cmatrix)
     778            8 :             IF (ALLOCATED(smatrix)) DEALLOCATE (smatrix)
     779              : 
     780              :          END DO ! ispin
     781            2 :          WRITE (iwfx, "(A,/)") "</Molecular Orbital Primitive Coefficients>"
     782              : 
     783              :          ! !===========================================================================!
     784              :          ! Energy
     785              :          ! !===========================================================================!
     786            2 :          WRITE (iwfx, "(A)") "<Energy>"
     787            2 :          WRITE (iwfx, "(ES26.16E3)") energy%total
     788            2 :          WRITE (iwfx, "(A,/)") "</Energy>"
     789              : 
     790              :          ! !===========================================================================!
     791              :          ! Virial ratio
     792              :          ! !===========================================================================!
     793            2 :          WRITE (iwfx, "(A)") "<Virial Ratio (-V/T)>"
     794            2 :          WRITE (iwfx, "(ES20.12)") - 1*(energy%total - energy%kinetic)/energy%kinetic
     795            2 :          WRITE (iwfx, "(A,/)") "</Virial Ratio (-V/T)>"
     796              : 
     797              :          ! !===========================================================================!
     798              :          ! Force-related quantities
     799              :          ! inactivated for now
     800              :          ! !===========================================================================!
     801              :          IF (ASSOCIATED(force) .AND. debug_this_module) THEN
     802              :             WRITE (iwfx, "(A)") "<Nuclear Cartesian Energy Gradients>"
     803              : 
     804              :             CALL get_qs_env(qs_env, force=force, atomic_kind_set=atomic_kind_set)
     805              :             ALLOCATE (atom_of_kind(num_atoms))
     806              :             ALLOCATE (kind_of(num_atoms))
     807              :             CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     808              :                                      atom_of_kind=atom_of_kind, &
     809              :                                      kind_of=kind_of)
     810              : 
     811              :             DO iatom = 1, num_atoms
     812              :                ikind = kind_of(iatom)
     813              :                i = atom_of_kind(iatom)
     814              :                WRITE (iwfx, "(3ES20.12E3)") force(ikind)%total(1:3, i)
     815              :             END DO
     816              : 
     817              :             WRITE (iwfx, "(A,/)") "<Nuclear Cartesian Energy Gradients>"
     818              : 
     819              :             ! W
     820              :             WRITE (iwfx, "(A)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
     821              :             WRITE (iwfx, "(A)") ""
     822              :             WRITE (iwfx, "(A,/)") "<Nuclear Virial of Energy-Gradient-Based Forces on Nuclei, W>"
     823              : 
     824              :             ! Full Virial ratio
     825              :             WRITE (iwfx, "(A)") "<Full Virial Ratio, -(V - W)/T)>"
     826              :             WRITE (iwfx, "(A)") ""
     827              :             WRITE (iwfx, "(A,/)") "</Full Virial Ratio, -(V - W)/T)>"
     828              : 
     829              :             DEALLOCATE (atom_of_kind)
     830              :             DEALLOCATE (kind_of)
     831              : 
     832              :          END IF
     833              : 
     834              :          ! !===========================================================================!
     835              :          ! Core-density
     836              :          ! !===========================================================================!
     837              :          ! Additional Electron Density Function (EDF)
     838              :          ! WRITE (iwfx, "(A)") "#<Additional Electron Density Function (EDF)>"
     839              :          ! WRITE (iwfx, "(A,/)") "#</Additional Electron Density Function (EDF)>"
     840              : 
     841              :          ! !-----------------------------------!
     842              :          ! Done
     843              :          ! !-----------------------------------!
     844            2 :          DEALLOCATE (atoms_cart)
     845            2 :          DEALLOCATE (atom_symbols)
     846            2 :          DEALLOCATE (atomic_numbers)
     847            2 :          DEALLOCATE (core_charges)
     848              : 
     849              :       ELSE  ! no output unit
     850            0 :          iwfx = -1
     851              :       END IF
     852              : 
     853            2 :       IF (output_unit > 0) THEN
     854              :          WRITE (output_unit, '(/,T2,A)') &
     855            1 :             '!--------------------------- Chargemol wfx written ---------------------------!'
     856              :       END IF
     857              : 
     858            2 :       CALL timestop(handle)
     859              : 
     860            4 :    END SUBROUTINE write_wfx
     861              : 
     862              : ! **************************************************************************************************
     863              : !> \brief ...
     864              : !> \param l_symb ...
     865              : !> \return ...
     866              : ! **************************************************************************************************
     867           76 :    INTEGER FUNCTION pgf_type_chargemol(l_symb) RESULT(label)
     868              :       !INTEGER, INTENT(OUT)                                 :: label
     869              :       CHARACTER(len=*), INTENT(IN)                       :: l_symb
     870              : 
     871              :       SELECT CASE (l_symb)
     872              :       CASE ("s")
     873           16 :          label = 1
     874              :       CASE ("px")
     875           16 :          label = 2
     876              :       CASE ("py")
     877           16 :          label = 3
     878              :       CASE ("pz")
     879           16 :          label = 4
     880              :       CASE ("dx2")
     881            2 :          label = 5
     882              :       CASE ("dy2")
     883            2 :          label = 6
     884              :       CASE ("dz2")
     885            2 :          label = 7
     886              :       CASE ("dxy")
     887            2 :          label = 8
     888              :       CASE ("dxz")
     889            2 :          label = 9
     890              :       CASE ("dyz")
     891            2 :          label = 10
     892              :       CASE ("fx3")
     893            0 :          label = 11
     894              :       CASE ("fy3")
     895            0 :          label = 12
     896              :       CASE ("fz3")
     897            0 :          label = 13
     898              :       CASE ("fx2y")
     899            0 :          label = 14
     900              :       CASE ("fx2z")
     901            0 :          label = 15
     902              :       CASE ("fy2z")
     903            0 :          label = 16
     904              :       CASE ("fxy2")
     905            0 :          label = 17
     906              :       CASE ("fxz2")
     907            0 :          label = 18
     908              :       CASE ("fyz2")
     909            0 :          label = 19
     910              :       CASE ("fxyz")
     911            0 :          label = 20
     912              :       CASE ("gx4")
     913            0 :          label = 21
     914              :       CASE ("gy4")
     915            0 :          label = 22
     916              :       CASE ("gz4")
     917            0 :          label = 23
     918              :       CASE ("gx3y")
     919            0 :          label = 24
     920              :       CASE ("gx3z")
     921            0 :          label = 25
     922              :       CASE ("gxy3")
     923            0 :          label = 26
     924              :       CASE ("gy3z")
     925            0 :          label = 27
     926              :       CASE ("gxz3")
     927            0 :          label = 28
     928              :       CASE ("gyz3")
     929            0 :          label = 29
     930              :       CASE ("gx2y2")
     931            0 :          label = 30
     932              :       CASE ("gx2z2")
     933            0 :          label = 31
     934              :       CASE ("gy2z2")
     935            0 :          label = 32
     936              :       CASE ("gx2yz")
     937            0 :          label = 33
     938              :       CASE ("gxy2z")
     939            0 :          label = 34
     940              :       CASE ("gxyz2")
     941            0 :          label = 35
     942              :       CASE ("hz5")
     943            0 :          label = 36
     944              :       CASE ("hyz4")
     945            0 :          label = 37
     946              :       CASE ("hy2z3")
     947            0 :          label = 38
     948              :       CASE ("hy3z2")
     949            0 :          label = 39
     950              :       CASE ("hy4z")
     951            0 :          label = 40
     952              :       CASE ("hy5")
     953            0 :          label = 41
     954              :       CASE ("hxz4")
     955            0 :          label = 42
     956              :       CASE ("hxyz3")
     957            0 :          label = 43
     958              :       CASE ("hxy2z2")
     959            0 :          label = 44
     960              :       CASE ("hxy3z")
     961            0 :          label = 45
     962              :       CASE ("hxy4")
     963            0 :          label = 46
     964              :       CASE ("hx2z3")
     965            0 :          label = 47
     966              :       CASE ("hx2yz2")
     967            0 :          label = 48
     968              :       CASE ("hx2y2z")
     969            0 :          label = 49
     970              :       CASE ("hx2y3")
     971            0 :          label = 50
     972              :       CASE ("hx3z2")
     973            0 :          label = 51
     974              :       CASE ("hx3yz")
     975            0 :          label = 52
     976              :       CASE ("hx3y2")
     977            0 :          label = 53
     978              :       CASE ("hx4z")
     979            0 :          label = 54
     980              :       CASE ("hx4y")
     981            0 :          label = 55
     982              :       CASE ("hx5")
     983            0 :          label = 56
     984              :       CASE default
     985           76 :          CPABORT("The chargemol interface supports basis functions up to l=5 only")
     986              :          !label = 99
     987              :       END SELECT
     988              : 
     989           76 :    END FUNCTION pgf_type_chargemol
     990              : 
     991              : END MODULE qs_chargemol
     992              : 
        

Generated by: LCOV version 2.0-1