23895

Committed 04 Sep 2025 11:02PM UTC coverage: 88.582% (-0.02%) from 88.597%

Build # 23895

Build Type

push

gitlab-ci

Committed by

Commit Message

Merge branch 'cleanup-td-prints' into 'development'

Cleanup of oscillator and rotatory strengths output

See merge request eT-program/eT!1578

Run Details

53979 of 60937 relevant lines covered (88.58%)

3115005.9 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

91.07

/src/io/input/input_tool_class.F90



! eT - a coupled cluster program
! Copyright (C) 2016-2024 the authors of eT
!
! eT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! eT is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <https://www.gnu.org/licenses/>.


module input_tool_class

   !!
   !! Input tool class module
   !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
   !!

   use kinds
   use section_class, only: section
   use global_out,    only: output

   use string_utilities, only: string_starts_with, first_index_of_substring
   use string_utilities, only: is_substring_in_string, is_empty, remove_substring_after_symbol


   type :: input_tool

      type(section),     allocatable, dimension(:), private :: sections

      character(len=30), allocatable, dimension(:), private :: rf_wfs
      character(len=30), allocatable, dimension(:), private :: cc_wfs
      character(len=30), allocatable, dimension(:), private :: ci_wfs

      integer, private :: n_keyword_lines ! Number of lines excluding the geometry section
      integer, private :: n_qm_atom_lines ! Number of QM atoms
      integer, private :: n_mm_atom_lines ! Number of MM atoms

      ! Arrays storing lines in the input
      character(len=200), dimension(:), allocatable, private :: keyword_lines
      character(len=200), dimension(:), allocatable, private :: qm_geometry
      character(len=200), dimension(:), allocatable, private :: mm_geometry

   contains

      procedure, public :: read_input
      procedure, public :: set_input
      procedure, public :: process_input

      procedure, private :: print_input_except_geometry

      procedure, public :: get_n_atoms &
                        => get_n_atoms_input_tool

      procedure, public :: get_n_mm_atoms &
                        => get_n_mm_atoms_input_tool

      procedure, public :: get_n_mm_molecules &
                        => get_n_mm_molecules_input_tool

      procedure, public :: get_geometry &
                        => get_geometry_input_tool

      procedure, public :: get_mm_geometry_fq &
                        => get_mm_geometry_fq_input_tool

      procedure, public :: get_mm_geometry_non_polarizable &
                        => get_mm_geometry_non_polarizable_input_tool

      procedure, public :: get_reference_wavefunction
      procedure, public :: get_cc_wavefunction
      procedure, public :: get_ci_wavefunction

      procedure, public :: requested_reference_calculation
      procedure, public :: requested_cc_calculation
      procedure, public :: requested_ci_calculation

      generic, public :: get_keyword                &
                      => get_32bit_integer_keyword, &
                         get_64bit_integer_keyword, &
                         get_string_keyword,        &
                         get_real_dp_keyword

      generic, public :: get_required_keyword                 &
                      => get_required_string_keyword,         &
                         get_required_32bit_integer_keyword,  &
                         get_required_64bit_integer_keyword,  &
                         get_required_real_dp_keyword

      procedure, public :: get_n_elements_for_keyword

      generic, public :: get_array_for_keyword          &
                      => get_integer_array_for_keyword, &
                         get_real_array_for_keyword,    &
                         get_string_array_for_keyword

      procedure, public :: string_in_array_keyword   => string_in_array_keyword

      procedure, public :: cleanup &
                        => cleanup_input_tool

      procedure, public :: is_section_present &
                        => is_section_present_input_tool

      procedure, public :: is_keyword_present &
                        => is_keyword_present_input_tool

      procedure, public :: is_embedding_on &
                        => is_embedding_on_input_tool

      procedure, public :: place_records_in_memory &
                        => place_records_in_memory_input_tool

      procedure, public :: get_n_state_guesses
      procedure, public :: get_state_guesses

      procedure, public :: get_present_keywords &
                        => get_present_keywords_input_tool

      procedure, private :: process_keyword_lines
      procedure, private :: set_section_keywords
      procedure, private :: check_for_errors

      procedure, private :: get_wavefunction
      procedure, private :: requested_calculation

      procedure, private :: get_section_index
      procedure, private :: print_sections

      procedure, private :: get_32bit_integer_keyword
      procedure, private :: get_64bit_integer_keyword
      procedure, private :: get_string_keyword
      procedure, private :: get_real_dp_keyword
      procedure, private :: get_required_string_keyword
      procedure, private :: get_required_32bit_integer_keyword
      procedure, private :: get_required_64bit_integer_keyword
      procedure, private :: get_required_real_dp_keyword

      procedure, private :: get_n_elements_in_set_or_range
      procedure, private :: get_integers_from_set_or_range

      procedure, private :: get_integer_array_for_keyword
      procedure, private :: get_real_array_for_keyword
      procedure, private :: get_string_array_for_keyword

      procedure, private, nopass :: expect_range
      procedure, private, nopass :: expect_set
      procedure, private, nopass :: check_get_array_for_keyword

      procedure, private :: is_section_recognized
      procedure, private :: section_not_recognized

   end type input_tool

   interface input_tool

      procedure :: new_input_tool

   end interface input_tool


contains


   function new_input_tool() result(this)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, June 2018
      !!
      implicit none

      type(input_tool) :: this

      type(section) :: active_atoms
      type(section) :: active_space
      type(section) :: boson
      type(section) :: calculations
      type(section) :: cc_mean_value
      type(section) :: cc_response
      type(section) :: cc_real_time
      type(section) :: ci_mean_value
      type(section) :: ci_transition_property
      type(section) :: dft_functional
      type(section) :: electric_field
      type(section) :: grid_info
      type(section) :: harmonic_frequencies
      type(section) :: hf_mean_value
      type(section) :: hf_response
      type(section) :: integrals
      type(section) :: lanczos
      type(section) :: memory
      type(section) :: method
      type(section) :: mlcc
      type(section) :: mlhf
      type(section) :: mm
      type(section) :: orbital_localization
      type(section) :: pcm
      type(section) :: print_
      type(section) :: solver_cc_es
      type(section) :: solver_cc_gs
      type(section) :: solver_cc_multipliers
      type(section) :: solver_cc_propagation
      type(section) :: solver_cc_response
      type(section) :: solver_cholesky
      type(section) :: solver_ci
      type(section) :: solver_cpp
      type(section) :: solver_fft_dipole_moment
      type(section) :: solver_fft_electric_field
      type(section) :: solver_geoopt
      type(section) :: solver_scf
      type(section) :: solver_hf_es
      type(section) :: solver_hf_response
      type(section) :: system
      type(section) :: visualization

      allocate(this%rf_wfs, source=[character(len=30) :: &
                                       'dft',            &
                                       'hf',             &
                                       'uhf',            &
                                       'mlhf',           &
                                       'qed-hf',         &
                                       'sc-qed-hf',      &
                                       'cuhf',           &
                                       'rohf'])

      allocate(this%cc_wfs, source=[character(len=30) :: &
                                       'ccs',            &
                                       'mp2',            &
                                       'cc2',            &
                                       'lowmem-cc2',     &
                                       'ccsd',           &
                                       'cc3',            &
                                       'ccsd(t)',        &
                                       'ccsdt',          &
                                       'mlcc2',          &
                                       'mlccsd',         &
                                       'qed-ccsd'])

      allocate(this%ci_wfs, source=[character(len=30) :: 'casci',     &
                                                         'fci',       &
                                                         'qed-casci', &
                                                         'qed-fci'])

      method = section(name_            = 'method', &
                       required         = .false.,  &
                       allowed_keywords = [this%rf_wfs, this%cc_wfs, this%ci_wfs])

      ! Set other sections in alphabetical order

      active_atoms = section(name_            = 'active atoms',       &
                             required         = .false.,              &
                             allowed_keywords = [character(len=30) :: &
                                                 'selection type',    &
                                                 'central atom',      &
                                                 'hf',                &
                                                 'ccs',               &
                                                 'cc2',               &
                                                 'ccsd',              &
                                                 'cc3',               &
                                                 'ccsd(t)',           &
                                                 'inactive basis',    &
                                                 'hf basis',          &
                                                 'ccs basis',         &
                                                 'cc2 basis',         &
                                                 'ccsd basis',        &
                                                 'cc3 basis',         &
                                                 'ccsd(t) basis'])


      active_space = section(name_            = 'active space',         &
                             required         = .false.,                &
                             allowed_keywords = [character(len=30) ::   &
                                                 'canonical',           &
                                                 'freeze core',         &
                                                 'freeze atom cores',   &
                                                 'localized'])


      boson = section(name_            = 'boson',                         &
                      required         = .false.,                         &
                      allowed_keywords = [character(len=30) ::            &
                                         'boson number',                  &
                                         'boson states',                  &
                                         'cbo',                           &
                                         'coherent state',                &
                                         'coupling bilinear',             &
                                         'coupling self',                 &
                                         'coupling',                      &
                                         'frequency',                     &
                                         'interaction type',              &
                                         'modes',                         &
                                         'polarization',                  &
                                         'print dipole eri construction', &
                                         'quadrupole oei',                &
                                         'wavevector'])


      calculations = section(name_            = 'do',                      &
                             required         = .false.,                   &
                             allowed_keywords = [character(len=30) ::      &
                                                'cholesky eri',            &
                                                'ground state',            &
                                                'geometry optimization',   &
                                                'harmonic frequencies',    &
                                                'excited state',           &
                                                'response',                &
                                                'mean value',              &
                                                'real time',               &
                                                'restart'])


      cc_mean_value = section(name_            = 'cc mean value',           &
                              required         = .false.,                   &
                              allowed_keywords = [character(len=30) ::      &
                                                 'dipole',                  &
                                                 'quadrupole',              &
                                                 'molecular gradient'])


      cc_response = section(name_            = 'cc response',         &
                            required         = .false.,               &
                            allowed_keywords = [character(len=30) ::  &
                                               'asymmetric',          &
                                               'damping',             &
                                               'dipole length',       &
                                               'dipole velocity',     &
                                               'dyson orbitals',      &
                                               'eom',                 &
                                               'frequencies',         &
                                               'initial states',      &
                                               'lr',                  &
                                               'permanent moments',   &
                                               'polarizabilities',    &
                                               'response components', &
                                               'rotatory strength',   &
                                               'variational energy',  &
                                               'transition moments',  &
                                               'full-space transition moments'])


      cc_real_time = section(name_            = 'cc real time',       &
                             required         = .false.,              &
                             allowed_keywords = [character(len=30) :: &
                                                'propagation',        &
                                                'fft dipole moment',  &
                                                'fft electric field'])


      ci_mean_value = section(name_            = 'ci mean value',     &
                              required         = .false.,             &
                              allowed_keywords = [character(len=30) ::&
                                                 'dipole',            &
                                                 'quadrupole'])


      ci_transition_property = section(name_            = 'ci transition property',  &
                                       required         = .false.,                   &
                                       allowed_keywords = [character(len=30) ::      &
                                                          'dipole',                  &
                                                          'quadrupole',              &
                                                          'final states',            &
                                                          'initial states'])


      dft_functional = section(name_            = 'dft functional',     &
                               required         = .false.,              &
                               allowed_keywords = [character(len=30) :: &
                                                  'correlation',        &
                                                  'exchange',           &
                                                  'functional',         &
                                                  'hf percentage'])


      electric_field = section(name_            = 'electric field',             &
                               required         = .false.,                      &
                               allowed_keywords = [character(len=30) ::         &
                                                  'envelope',                   &
                                                  'x polarization',             &
                                                  'y polarization',             &
                                                  'z polarization',             &
                                                  'central time',               &
                                                  'width',                      &
                                                  'carrier angular frequency',  &
                                                  'peak strength',              &
                                                  'phase shift',                &
                                                  'repetition',                 &
                                                  'separation'])


      harmonic_frequencies = section(name_            = 'harmonic frequencies',      &
                                     required         = .false.,                     &
                                     allowed_keywords = [character(len=30) ::        &
                                                        'gradient displacement',     &
                                                        'gradient method',           &
                                                        'hessian displacement',      &
                                                        'run geometry optimization', &
                                                        'state',                     &
                                                        'wigner samples',            &
                                                        'wigner seed'])


      grid_info = section(name_            = 'grid info',             &
                          required         = .false.,                 &
                          allowed_keywords = [character(len=30) ::    &
                                             'cube side',             &
                                             'cube offset',           &
                                             'minimum angular order', &
                                             'maximum angular order', &
                                             'partitioning',          &
                                             'radial threshold',      &
                                             'quadrature'])


      hf_mean_value = section(name_            = 'hf mean value',      &
                              required         = .false.,              &
                              allowed_keywords = [character(len=30) :: &
                                                 'dipole',             &
                                                 'molecular gradient', &
                                                 'quadrupole'])

      hf_response = section(name_            = 'hf response',        &
                            required         = .false.,              &
                            allowed_keywords = [character(len=30) :: &
                                              'shielding', &
                                              'dipole velocity',        &
                                              'rotatory strength',      &
                                              'polarizabilities'])


      integrals = section(name_            = 'integrals',          &
                          required         = .false.,              &
                          allowed_keywords = [character(len=30) :: &
                                             'cholesky storage',   &
                                             'eri storage',        &
                                             'mo eri in memory',   &
                                             'ri',                 &
                                             't1 eri in memory'])


      lanczos = section(name_            = 'lanczos',                &
                        required         = .false.,                  &
                        allowed_keywords = [character(len=30) ::     &
                                           'biorthogonalize last',   &
                                           'chain length',           &
                                           'deflation threshold',    &
                                           'lanczos normalization',  &
                                           'max energy',             &
                                           'min energy',             &
                                           'min transition strength',&
                                           'operators',              &
                                           'overlap threshold',      &
                                           'restart states',         &
                                           'skip convergence'])


      memory = section(name_            = 'memory',             &
                       required         = .false.,              &
                       allowed_keywords = [character(len=30) :: &
                                          'available',          &
                                          'unit'])


      mlcc = section(name_            = 'mlcc',                   &
                     required         = .false.,                  &
                     allowed_keywords = [character(len=30) ::     &
                                        'levels',                 &
                                        'cc2 orbitals',           &
                                        'ccsd orbitals',          &
                                        'cholesky threshold',     &
                                        'nto restart',            &
                                        'cnto restart',           &
                                        'orbital restart',        &
                                        'cnto occupied cc2',      &
                                        'cnto virtual cc2',       &
                                        'cnto occupied ccsd',     &
                                        'cnto virtual ccsd',      &
                                        'cnto states',            &
                                        'nto states',             &
                                        'nto occupied cc2',       &
                                        'nto occupied ccsd',      &
                                        'canonical virtual cc2',  &
                                        'canonical virtual ccsd', &
                                        'print ccs calculation',  &
                                        'print cc2 calculation'])


      mlhf = section(name_            = 'multilevel hf',            &
                     required         = .false.,                    &
                     allowed_keywords = [character(len=30) ::       &
                                        'initial hf optimization',  &
                                        'initial hf threshold',     &
                                        'print initial hf',         &
                                        'cholesky threshold',       &
                                        'project on minimal basis', &
                                        'cholesky virtuals',        &
                                        'no mo screening'])


      mm = section(name_            = 'molecular mechanics', &
                   required         = .false.,               &
                   allowed_keywords = [character(len=30) ::  &
                                      'forcefield',          &
                                      'algorithm '])


      orbital_localization = section(name_            = 'orbital localization', &
                                     required         = .false.,                &
                                     allowed_keywords = [character(len=30) ::   &
                                                         'orbitals',            &
                                                         'threshold',           &
                                                         'type'])


      pcm = section(name_            = 'pcm',                &
                    required         = .false.,              &
                    allowed_keywords = [character(len=30) :: &
                                       'solvent',            &
                                       'input',              &
                                       'tesserae area',      &
                                       'solver type'])


      print_ = section(name_            = 'print',               &
                       required         = .false.,               &
                       allowed_keywords = [character(len=30) ::  &
                                          'output print level ', &
                                          'timing print level ', &
                                          'full references',     &
                                          'z-matrix'])


      solver_cpp = section(name_    = 'solver cpp',                    &
                           required = .false.,                         &
                           allowed_keywords = [character(len=30) ::    &
                                              'max iterations',        &
                                              'max reduced dimension', &
                                              'residual minimization', &
                                              'residual threshold',    &
                                              'restart',               &
                                              'storage'])


      solver_cc_es = section(name_            = 'solver cc es',          &
                             required         = .false.,                 &
                             allowed_keywords = [character(len=30) ::    &
                                                'algorithm',             &
                                                'olsen',                 &
                                                'core excitation',       &
                                                'electron attachment',   &
                                                'ionization',            &
                                                'energy threshold ',     &
                                                'crop',                  &
                                                'residual threshold',    &
                                                'max iterations',        &
                                                'restart',               &
                                                'left eigenvectors',     &
                                                'right eigenvectors',    &
                                                'storage',               &
                                                'singlet states',        &
                                                'triplet states',        &
                                                'diis dimension',        &
                                                'max reduced dimension', &
                                                'max micro iterations',  &
                                                'rel micro threshold',   &
                                                'remove core',           &
                                                'state guesses'])


      solver_cc_gs = section(name_            = 'solver cc gs',            &
                             required         = .false.,                   &
                             allowed_keywords = [character(len=30) ::      &
                                                'algorithm',               &
                                                'energy threshold',        &
                                                'residual threshold',      &
                                                'crop',                    &
                                                'micro iteration storage', &
                                                'max micro iterations',    &
                                                'multimodel newton',       &
                                                'rel micro threshold',     &
                                                'storage',                 &
                                                'max iterations',          &
                                                'diis dimension',          &
                                                'restart' ])


      solver_cc_multipliers = section(name_            = 'solver cc multipliers',   &
                                      required         = .false.,                   &
                                      allowed_keywords = [character(len=30) ::      &
                                                         'algorithm',               &
                                                         'threshold',               &
                                                         'storage',                 &
                                                         'crop',                    &
                                                         'micro iteration storage', &
                                                         'max micro iterations',    &
                                                         'multimodel newton',       &
                                                         'rel micro threshold',     &
                                                         'diis dimension',          &
                                                         'restart',                 &
                                                         'max reduced dimension',   &
                                                         'max iterations'])


      solver_cc_propagation = section(name_            = 'solver cc propagation',    &
                                      required         = .false.,                    &
                                      allowed_keywords = [character(len=30) ::       &
                                                         'initial time',             &
                                                         'final time',               &
                                                         'time step',                &
                                                         'steps between output',     &
                                                         'step reduction factor',    &
                                                         'implicit threshold',       &
                                                         'max error',                &
                                                         'min error',                &
                                                         'energy output',            &
                                                         'dipole moment output',     &
                                                         'electric field output',    &
                                                         'parameters output',        &
                                                         'mo density matrix output', &
                                                         'method',                   &
                                                         'integrator'])


      solver_cc_response = section(name_            = 'solver cc response', &
                                   required         = .false.,              &
                                   allowed_keywords = [character(len=30) :: &
                                                      'threshold',          &
                                                      'gradient response threshold', &
                                                      'storage',            &
                                                      'max iterations'])


      solver_cholesky = section(name_            = 'solver cholesky',    &
                                required         = .false.,              &
                                allowed_keywords = [character(len=30) :: &
                                                   'threshold',          &
                                                   'span',               &
                                                   'batches',            &
                                                   'qualified',          &
                                                   'one center',         &
                                                   'diagonal test',      &
                                                   'mo screening'])


      solver_ci = section(name_            = 'solver ci',              &
                          required         = .false.,                  &
                          allowed_keywords = [character(len=30) ::     &
                                              'energy threshold ',     &
                                              'max iterations',        &
                                              'max reduced dimension', &
                                              'residual threshold',    &
                                              'restart',               &
                                              'start guess',           &
                                              'states',                &
                                              'storage'])

      solver_fft_dipole_moment = section(name_            = 'solver fft dipole moment', &
                                         required         = .false.,                    &
                                         allowed_keywords = [character(len=30) ::       &
                                                            'initial time',             &
                                                            'final time',               &
                                                            'time step',                &
                                                            'padding initial time',     &
                                                            'hann window',              &
                                                            'rectangular window'])


      solver_fft_electric_field = section(name_            = 'solver fft electric field', &
                                          required         = .false.,                     &
                                          allowed_keywords = [character(len=30) ::        &
                                                             'initial time',              &
                                                             'final time',                &
                                                             'time step',                 &
                                                             'padding initial time',      &
                                                             'hann window',               &
                                                             'rectangular window'])


      solver_geoopt = section(name_            = 'solver geometry optimization',              &
                              required         = .false.,                      &
                              allowed_keywords = [character(len=30) ::         &
                                                 'algorithm',                  &
                                                 'energy threshold',           &
                                                 'residual threshold',         &
                                                 'residual rms threshold',     &
                                                 'state',                      &
                                                 'max iterations',             &
                                                 'max step',                   &
                                                 'gradient method',            &
                                                 'restart',                    &
                                                 'step size',                  &
                                                 'displacement threshold',     &
                                                 'displacement rms threshold'])


      solver_scf = section(name_            = 'solver scf',                  &
                           required         = .false.,                       &
                           allowed_keywords = [character(len=30) ::          &
                                              'algorithm',                   &
                                              'ao density guess',            &
                                              'coulomb exchange terms',      &
                                              'coulomb threshold',           &
                                              'crop',                        &
                                              'cumulative fock threshold',   &
                                              'diis dimension',              &
                                              'diabatize orbitals',          &
                                              'energy threshold',            &
                                              'exchange threshold',          &
                                              'residual threshold',          &
                                              'gradient response threshold', &
                                              'max iterations',              &
                                              'population analysis',         &
                                              'write orbitals',              &
                                              'restart',                     &
                                              'rohf coupling parameters',    &
                                              'skip',                        &
                                              'storage',                     &
                                              'write molden'])


      solver_hf_es = section(name_            = 'solver hf es',         &
                               required         = .false.,                  &
                               allowed_keywords = [character(len=30) ::     &
                                                  'energy threshold',       &
                                                  'max iterations',         &
                                                  'max reduced dimension',  &
                                                  'residual threshold',     &
                                                  'restart',                &
                                                  'storage',                &
                                                  'singlet states',         &
                                                  'tamm-dancoff'])


      solver_hf_response = section(name_           = 'solver hf response',       &
                                  required         = .false.,                    &
                                  allowed_keywords = [character(len=30) ::       &
                                                        'frequencies',           &
                                                        'max iterations',        &
                                                        'max reduced dimension', &
                                                        'print iterations',      &
                                                        'residual threshold',    &
                                                        'restart',               &
                                                        'storage'])


      system = section(name_            = 'system',               &
                       required         = .false.,                &
                       allowed_keywords = [character(len=30) ::   &
                                          'cartesian gaussians',  &
                                          'write orbitals',       &
                                          'pure gaussians',       &
                                          'charge',               &
                                          'multiplicity'])


      visualization = section(name_            = 'visualization',             &
                              required         = .false.,                     &
                              allowed_keywords = [character(len=30) ::        &
                                                 'file format',               &
                                                 'grid spacing',              &
                                                 'grid buffer',               &
                                                 'grid min',                  &
                                                 'grid max',                  &
                                                 'plot cc density',           &
                                                 'plot cc orbitals',          &
                                                 'plot dyson orbitals',       &
                                                 'plot es densities',         &
                                                 'plot hf orbitals',          &
                                                 'plot hf density',           &
                                                 'plot hf active density',    &
                                                 'plot cntos',                &
                                                 'plot ntos',                 &
                                                 'nto threshold',             &
                                                 'plot transition densities', &
                                                 'states to plot'])

      ! Gather all sections into the file's section array

      this%sections = [active_atoms,              &
                       active_space,              &
                       boson,                     &
                       calculations,              &
                       cc_mean_value,             &
                       cc_response,               &
                       cc_real_time,              &
                       ci_mean_value,             &
                       ci_transition_property,    &
                       dft_functional,            &
                       electric_field,            &
                       grid_info,                 &
                       harmonic_frequencies,      &
                       hf_mean_value,             &
                       hf_response,               &
                       integrals,                 &
                       lanczos,                   &
                       memory,                    &
                       method,                    &
                       mlcc,                      &
                       mlhf,                      &
                       mm,                        &
                       orbital_localization,      &
                       pcm,                       &
                       print_,                    &
                       solver_cc_es,              &
                       solver_cc_gs,              &
                       solver_cc_multipliers,     &
                       solver_cc_propagation,     &
                       solver_cc_response,        &
                       solver_cholesky,           &
                       solver_ci,                 &
                       solver_cpp,                &
                       solver_fft_dipole_moment,  &
                       solver_fft_electric_field, &
                       solver_geoopt,             &
                       solver_scf,                &
                       solver_hf_es,            &
                       solver_hf_response,      &
                       system,                    &
                       visualization]

   end function new_input_tool


   subroutine read_input(this, file_name)
      !!
      !! Written by Eirik F. Kjønstad, May 2022
      !!
      use input_file_parser_class, only: input_file_parser

      implicit none

      class(input_tool), intent(inout) :: this

      character(len=*), intent(in) :: file_name

      type(input_file_parser), allocatable :: file_parser

      file_parser = input_file_parser(file_name)

      call file_parser%read_and_process_input_file()
      call file_parser%get_keyword_lines(this%keyword_lines, this%n_keyword_lines)
      call file_parser%get_qm_geometry_lines(this%qm_geometry, this%n_qm_atom_lines)

      if (file_parser%is_mm_geometry_present()) &
         call file_parser%get_mm_geometry_lines(this%mm_geometry, this%n_mm_atom_lines)

      deallocate(file_parser)

   end subroutine read_input


   subroutine process_input(this)
      !!
      !! Written by Eirik F. Kjønstad, Mar 2024
      !!
      implicit none

      class(input_tool), intent(inout) :: this

      call this%process_keyword_lines()
      call this%print_input_except_geometry()
      call this%check_for_errors()

   end subroutine process_input


   subroutine print_input_except_geometry(this)
      !!
      !! Written by Eirik F. Kjønstad, May 2022
      !!
      implicit none

      class(input_tool), intent(in) :: this

      integer :: n

      call output%printf('m', ':: Input file', fs='(//t3,a)')
      call output%print_separator('m', 13, '=', fs='(t3,a/)')

      call output%printf('m', 'Note: geometry section is excluded from this print.', fs='(t6,a/)')

      do n = 1, this%n_keyword_lines

         call output%printf('m', this%keyword_lines(n), fs='(t6,a)', ll=120)

      enddo

   end subroutine print_input_except_geometry


   subroutine set_input(this, keywords, geometry, mm_geometry)
      !!
      !! Written by Eirik F. Kjønstad, Mar 2024
      !!
      implicit none

      class(input_tool) :: this

      character(len=200), dimension(:), intent(in), optional :: keywords
      character(len=200), dimension(:), intent(in), optional :: geometry
      character(len=200), dimension(:), intent(in), optional :: mm_geometry

      if (present(geometry)) then

         this%qm_geometry = geometry
         this%n_qm_atom_lines = size(geometry)

      else

         allocate(this%qm_geometry, &
                  source=[character(len=200) :: '! no geometry'])

         this%n_qm_atom_lines = 1

      endif

      if (present(keywords)) then

         this%keyword_lines = keywords
         this%n_keyword_lines = size(keywords)

      else

         allocate(this%keyword_lines, &
                  source=[character(len=200) :: '! no input keywords'])

         this%n_keyword_lines = 1

      endif

      if (present(mm_geometry)) then

         this%mm_geometry = mm_geometry
         this%n_mm_atom_lines = size(mm_geometry)

      endif

   end subroutine set_input


   subroutine process_keyword_lines(this)
      !!
      !! Written by Sara Angelico and Alexander C. Paul, 2023
      !!
      implicit none

      class(input_tool), intent(inout) :: this

      integer :: start_, end_, i, j

      character(len=200) :: section, line, section_line

      do i = 1, this%n_keyword_lines

         line = trim(adjustl(this%keyword_lines(i)))

         if (.not. string_starts_with(line, '-')) cycle

         start_ = i + 1
         section = trim(adjustl(line(2:)))

         section = remove_substring_after_symbol(section,"!")

         if (.not. this%is_section_recognized(section)) call this%section_not_recognized(section)

         end_ = this%n_keyword_lines
         do j = start_, this%n_keyword_lines

            section_line = trim(adjustl(this%keyword_lines(j)))

            if (.not. string_starts_with(section_line, '-')) cycle

            end_ = j - 1
            exit

         enddo

         call this%set_section_keywords(section, this%keyword_lines(start_:end_))

      enddo

   end subroutine process_keyword_lines


   subroutine print_sections(this)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      integer :: k

      call output%printf('m', 'The valid input sections are:', fs='(/t3,a/)')

      do k = 1, size(this%sections)

         call output%printf('m', '(a0)', chars=[trim(this%sections(k)%name_)], fs='(t6,a)')

      enddo

   end subroutine print_sections


   function requested_reference_calculation(this) result(requested)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      logical :: requested

      requested = this%requested_calculation(this%rf_wfs)

   end function requested_reference_calculation


   function requested_calculation(this, wfs) result(requested)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      !! Looks for wfs(k) in method and returns true if it finds one such wf
      !! in the input file. Gives an error if it finds more than one of wfs(k).
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=30), dimension(:), intent(in) :: wfs

      logical :: requested

      integer :: n_wfs, k, section_index

      n_wfs = 0

      section_index = this%get_section_index('method')

      do k = 1, size(wfs)

         if (this%sections(section_index)%is_keyword_present(wfs(k))) then

            n_wfs = n_wfs + 1

         endif

      enddo

      requested = .false.

      if (n_wfs == 1) then

         requested = .true.

      elseif (n_wfs > 1) then

         call output%error_msg('Requested more than one method &
                               &of a single kind (e.g., two CC methods).')

      endif

   end function requested_calculation


   function requested_cc_calculation(this) result(requested)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      logical :: requested

      requested = this%requested_calculation(this%cc_wfs)

   end function requested_cc_calculation


   function requested_ci_calculation(this) result(requested)
      !!
      !! Written by Enrico Ronca, 2020
      !!
      implicit none

      class(input_tool), intent(in) :: this

      logical :: requested

      requested = this%requested_calculation(this%ci_wfs)

   end function requested_ci_calculation


   function get_reference_wavefunction(this) result(ref_wf)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=30) :: ref_wf

      ref_wf = this%get_wavefunction(this%rf_wfs)

   end function get_reference_wavefunction


   function get_cc_wavefunction(this) result(cc_wf)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=30) :: cc_wf

      cc_wf = this%get_wavefunction(this%cc_wfs)

   end function get_cc_wavefunction


   function get_ci_wavefunction(this) result(ci_wf)
      !!
      !! Written by Enrico Ronca, Apr 2020
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=30) :: ci_wf

      ci_wf = this%get_wavefunction(this%ci_wfs)

   end function get_ci_wavefunction


   function get_wavefunction(this, wfs) result(wf)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Apr 2019
      !!
      !! Looks for wfs(k) in the 'method' section of the input.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=30), dimension(:), intent(in) :: wfs

      character(len=30) :: wf

      integer :: k, section_index

      section_index = this%get_section_index('method')

      do k = 1, size(wfs)

         if (this%sections(section_index)%is_keyword_present(wfs(k))) then

            wf = wfs(k)
            return

         endif

      enddo

      call output%error_msg('Failed to read wavefunction from input. &
         & Did you forget to specify all wavefunctions for the calculation?')

   end function get_wavefunction


   subroutine get_32bit_integer_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! If specified, reads keyword as an integer into keyword value.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer(i32), intent(inout) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) return

      read(keyword_value_string, *) keyword_value

   end subroutine get_32bit_integer_keyword


   subroutine get_64bit_integer_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! If specified, reads keyword as an integer into keyword value.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer(i64), intent(inout) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) return

      read(keyword_value_string, *) keyword_value

   end subroutine get_64bit_integer_keyword


   subroutine get_required_32bit_integer_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Reads keyword as an integer into keyword value. If the keyword or the
      !! section is not specified, an error occurs because the keyword is "required".
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer(i32), intent(inout) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) then
         call output%error_msg('could not find the required keyword '// trim(keyword) // &
                               ' in section ' // trim(section))

      endif

      read(keyword_value_string, *) keyword_value

   end subroutine get_required_32bit_integer_keyword


   subroutine get_required_64bit_integer_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Reads keyword as an integer into keyword value. If the keyword or the
      !! section is not specified, an error occurs because the keyword is "required".
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer(i64), intent(out) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) then
         call output%error_msg('could not find the required keyword '// trim(keyword) // &
                               ' in section ' // trim(section))

      endif

      read(keyword_value_string, *) keyword_value

   end subroutine get_required_64bit_integer_keyword


   subroutine get_real_dp_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! If specified, reads keyword as a real double precision into keyword value.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      real(dp), intent(inout) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) return

      read(keyword_value_string, *) keyword_value

   end subroutine get_real_dp_keyword


   subroutine get_required_real_dp_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Reads keyword as an double precision into keyword value. If the keyword or the
      !! section is not specified, an error occurs because the keyword is "required".
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      real(dp), intent(out) :: keyword_value

      character(len=200) :: keyword_value_string

      keyword_value_string = ""
      call this%get_string_keyword(keyword, section, keyword_value_string)

      if (is_empty(keyword_value_string)) then
         call output%error_msg('could not find the required keyword '// trim(keyword) // &
                               ' in section ' // trim(section))

      endif

      read(keyword_value_string, *) keyword_value

   end subroutine get_required_real_dp_keyword


   subroutine get_string_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! If specified, reads keyword as a string into keyword value.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      character(len=200), intent(inout) :: keyword_value
      character(len=200) :: keyword_value_string

      integer :: section_index

      section_index = this%get_section_index(section)

      keyword_value_string = this%sections(section_index)%get_string_keyword(keyword)

      ! Make sure not to overwrite defaults set on the outside
      if (is_empty(keyword_value_string)) return

      keyword_value = keyword_value_string

   end subroutine get_string_keyword


   subroutine get_required_string_keyword(this, keyword, section, keyword_value)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! If specified, reads keyword as a string into keyword value.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      character(len=200) :: keyword_value

      keyword_value = ""
      call this%get_string_keyword(keyword, section, keyword_value)

      if (is_empty(keyword_value)) then
         call output%error_msg('could not find the required keyword '// trim(keyword) // &
                               ' in section ' // trim(section))

      endif

   end subroutine get_required_string_keyword


   function get_section_index(this, section) result(index)
      !!
      !! Written by Sara Angelico and Alexander C. Paul, Feb 2023
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: section

      integer :: index, i

      index = 0

      do i = 1, size(this%sections)
         if(trim(this%sections(i)%name_) == trim(section)) then
            index = i
            return
         end if
      enddo

      call this%section_not_recognized(section)

   end function get_section_index


   function is_keyword_present_input_tool(this, keyword, section) result(is_present)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Returns true if the keyword is in the section.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      logical :: is_present

      integer :: section_index

      section_index = this%get_section_index(section)

      is_present = this%sections(section_index)%is_keyword_present(keyword)

   end function is_keyword_present_input_tool


   subroutine get_present_keywords_input_tool(this, section, keywords)
      !!
      !! Written by Leo Stoll, 2025
      !!
      !! For a given section returns the keywords given in the input file.
      !!
      implicit none

      class(input_tool), intent(inout) :: this

      character(len=*), intent(in) :: section

      character(len=200), dimension(:), allocatable, intent(inout) :: keywords

      integer :: section_index

      section_index = this%get_section_index(section)

      keywords = this%sections(section_index)%get_present_keywords()

   end subroutine get_present_keywords_input_tool


   function is_section_present_input_tool(this, section) result(is_present)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Returns true if the section exists, false if it doesn't.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: section

      logical :: is_present

      integer :: section_index

      section_index = this%get_section_index(section)

      is_present = this%sections(section_index)%is_section_present()

   end function is_section_present_input_tool


   function get_n_elements_for_keyword(this, keyword, section) result(n_elements)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Gets the number of elements specified for a keyword in set/range notation
      !!
      !! Ranges should always be given as [a,b].
      !! Lists should always be given as {a, b, c, d}.
      !!
      !! Function is called in preparation for get_array_for_keyword
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer :: n_elements

      character(len=200) :: keyword_value_string

      n_elements = 0

      call this%get_keyword(keyword, section, keyword_value_string)

      n_elements = this%get_n_elements_in_set_or_range(keyword_value_string)

   end function get_n_elements_for_keyword


   subroutine get_integer_array_for_keyword(this, keyword, section, n_elements, array_)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Gets input variable array (array_) for keyword which is specified on input
      !! by either an integer range or list (of length n_elements).
      !!
      !! Ranges should always be given as [a,b].
      !! Lists should always be given as {a, b, c, d}.
      !!
      !! Routine should be called after the get_n_elements_for_keyword is called
      !! in order to determine n_elements so that array_ can be allocated.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section

      integer, intent(in) :: n_elements

      integer, dimension(n_elements) :: array_

      character(len=200) :: keyword_value_string

      call check_get_array_for_keyword(keyword, section, n_elements)

      call this%get_keyword(keyword, section, keyword_value_string)

      call this%get_integers_from_set_or_range(keyword_value_string, n_elements, array_)

   end subroutine get_integer_array_for_keyword


   subroutine get_real_array_for_keyword(this, keyword, section, n_elements, array_)
      !!
      !! Written by Sarai D. Folkestad, Eirik F. Kjønstad, Andreas Skeidsvoll and Alexander C. Paul, 2019-2023
      !!
      !! Gets input variable array (array_) for keyword
      !! which is specified on input by a list
      !!
      !! Lists should always be given as {a, b, c, d}
      !!
      !! Routine should be called after the get_n_elements_for_keyword is called
      !! in order to determine n_elements so that array_ can be allocated.
      !!
      use string_utilities, only: read_comma_separated_list

      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section
      integer,          intent(in) :: n_elements

      real(dp), dimension(n_elements), intent(inout) :: array_

      character(len=200) :: keyword_value_string
      character(len=:), allocatable :: string

      call check_get_array_for_keyword(keyword, section, n_elements)

      call this%get_keyword(keyword, section, keyword_value_string)

      string = trim(adjustl(keyword_value_string))
      call this%expect_set(string)

      call read_comma_separated_list(string(2:len(string)-1), array_, n_elements)

   end subroutine get_real_array_for_keyword


   subroutine get_string_array_for_keyword(this, keyword, section, n_elements, array_)
      !!
      !! Written by Anna Kristina Schnack-Petersen and Alexander C. Paul, 2022-2023
      !!
      !! Gets input variable array (array_) for keyword which is specified on input
      !! by a list of strings (of length n_elements).
      !!
      !! Lists should always be given as {a, b, c, d}
      !!
      !! Routine should be called after the get_n_elements_for_keyword is called
      !! in order to determine n_elements so that array_ can be allocated.
      !!
      use string_utilities, only: split_at_delimiter

      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: section
      integer,          intent(in) :: n_elements

      character(len=200), dimension(n_elements), intent(inout) :: array_

      character(len=200)            :: keyword_value_string
      character(len=:), allocatable :: string

      call check_get_array_for_keyword(keyword, section, n_elements)

      call this%get_keyword(keyword, section, keyword_value_string)

      string = trim(adjustl(keyword_value_string))
      call this%expect_set(string)

      call split_at_delimiter(string(2:len(string)-1), n_elements-1, array_, ',')

   end subroutine get_string_array_for_keyword


   pure function get_n_atoms_input_tool(this) result(n_atoms)
      !!
      !! Get n atoms
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Reads the geometry section of the input and counts the number of atoms
      !!
      implicit none

      class(input_tool), intent(in) :: this

      integer :: n_atoms

      integer :: record

      n_atoms = 0

      do record = 1, this%n_qm_atom_lines

         if (this%qm_geometry(record)(1:6) .ne. 'basis:' &
            .and. (this%qm_geometry(record)(1:6) .ne. 'units:') &
            .and. (this%qm_geometry(record)(1:5) .ne. 'ghost')) &
            n_atoms = n_atoms + 1

      enddo

   end function get_n_atoms_input_tool


   pure function get_n_mm_atoms_input_tool(this) result(n_atoms)
      !!
      !! Get n MM atoms
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      implicit none

      class(input_tool), intent(in) :: this

      integer :: n_atoms

      n_atoms = this%n_mm_atom_lines

   end function  get_n_mm_atoms_input_tool


   function get_n_mm_molecules_input_tool(this) result(n_molecules)
      !!
      !! Written by Sarai D. Folkestad, Sep 2020
      !!
      !! Reads the MM geometry and counts the number of molecules
      !!
      implicit none

      class(input_tool), intent(in) :: this

      integer :: n_molecules

      integer :: cursor, record, current_molecule, previous_molecule

      character(len=200) :: string, imolecule

      previous_molecule = 1
      n_molecules       = 1

      do record = 1, this%n_mm_atom_lines

         string = (trim(adjustl(this%mm_geometry(record))))
         cursor = first_index_of_substring(string,'=')

         string = string(cursor+1:200)

         cursor = first_index_of_substring(string,']')

         imolecule = string(1:cursor-1)
         read(imolecule,'(i4)') current_molecule

         if (record == 1 .and. current_molecule .ne. 1) &
            call output%error_msg('Molecule 1 is missing!')

         if (current_molecule .eq. previous_molecule + 1) then

            n_molecules = n_molecules + 1

         elseif (current_molecule .ne. previous_molecule .and. &
               current_molecule .ne. previous_molecule + 1) then

            call output%error_msg('Molecule (i0) missing even though (i0) &
                                 &should be present in QM/MM calculation.', &
                                 ints=[previous_molecule + 1, n_molecules])

         endif

         previous_molecule = current_molecule

      enddo

   end function get_n_mm_molecules_input_tool


   subroutine get_geometry_input_tool(this, n_atoms, symbols,  &
                                      positions, basis_sets,   &
                                      charge, units_angstrom,  &
                                      is_ghost)
      !!
      !! Written by Sarai D. Folkestad, Eirik F. Kjønstad
      !! Åsmund H. Tveten and Tor S. Haugland, 2019-2021
      !!
      !! Reads the geometry from the output file and sets it in the list of atoms.
      !!
      !! Note: This routine should be called after a call to the function
      !! get_n_atoms_input_tool which finds the number of atoms in the system
      !!
      use string_utilities, only: convert_character_to_uppercase

      implicit none

      class(input_tool), intent(in) :: this

      integer, intent(in) :: n_atoms

      character(len=2), dimension(n_atoms), intent(out)   :: symbols
      character(len=100), dimension(n_atoms), intent(out) :: basis_sets
      logical, dimension(n_atoms), intent(out) :: is_ghost

      real(dp), dimension(3, n_atoms), intent(out) :: positions ! x, y, z
      integer, dimension(n_atoms), intent(out) :: charge

      logical, intent(out) :: units_angstrom

      integer :: current_atom, i

      integer :: start_

      character(len=200) :: string
      character(len=100) :: current_basis
      logical :: is_ghost_atom = .false.

      integer :: cursor

      start_ = 1 ! Specifies the line of the first and required basis

      ! Are units specified?
      ! Note that units can only be specified as the first line of the geometry

      units_angstrom = .true. ! Default units are Angstrom
      charge = 0 ! Default charge is zero

      if (this%qm_geometry(1)(1:6) == 'units:') then

         string = (trim(adjustl(this%qm_geometry(1)(7:200))))

         if (string(1:4) == 'bohr') then

            units_angstrom = .false.

         elseif (string(1:8) /= 'angstrom') then

            call output%error_msg('units of atom coordinates must be either angstrom or bohr')

         endif

         start_ = 2

      endif

      ! Error if next line is not a basis set line

      if(this%qm_geometry(start_)(1:6) /= 'basis:') &
            call output%error_msg('did not find basis in geometry section.')

      current_atom = 0

      do i = start_, this%n_qm_atom_lines

         string = trim(adjustl(this%qm_geometry(i)))

         if (string(1:6) == 'units:') &
            call output%error_msg('Units must be specified as the first line in the geometry section.')

         if(string(1:6) == 'basis:') then

            current_basis = trim(adjustl(string(7:200)))
            current_basis = remove_substring_after_symbol(current_basis,"!")

         elseif (string(1:5) == 'ghost') then

            ! Every atom after 'ghost' keyword are ghosts
            is_ghost_atom = .true.

         else

            current_atom = current_atom + 1

            basis_sets(current_atom) = current_basis
            symbols(current_atom)    = string(1:2)
            is_ghost(current_atom)   = is_ghost_atom

            string = adjustl(string(3:200))

            if (is_substring_in_string(string, 'q')) then

               cursor = first_index_of_substring(string, 'q')
               read(string(1 : cursor - 1), *) positions(:, current_atom)

               string = adjustl(string(cursor+1:200))

               if (.not. string_starts_with(string, '=')) &
                     call output%error_msg('in assignment of charge to atom in the geometry')

               string = adjustl(string(2:200))
               read(string, *) charge(current_atom)
            else

               read(string(1:), *) positions(:, current_atom)

            endif

         endif

      enddo

      ! First character of symbol should be upper case

      do i = 1, n_atoms

         symbols(i)(1:1) = convert_character_to_uppercase(symbols(i)(1:1))

      enddo

   end subroutine get_geometry_input_tool


   subroutine get_mm_geometry_fq_input_tool(this, n_atoms, n_molecules, &
                                          n_atoms_per_molecule, symbols,   &
                                          positions, chi, eta)
      !!
      !! Written by Tommaso Giovanini and Sarai D. Folkestad, 2019-2020
      !!
      !! Reads the geometry of the MM portion from the input file in the case
      !! of an fq calculation and sets it in the list of atoms.
      !!
      !! Note: In order to be run, you need to know the number of MM atoms and molecules
      !!       in the system
      !!
      use string_utilities, only: convert_character_to_uppercase
      use string_utilities, only: last_index_of_substring

      implicit none

      class(input_tool),   intent(in) :: this
      integer,             intent(in) :: n_atoms
      integer,             intent(in) :: n_molecules

      character(len=2), dimension(n_atoms),     intent(out) :: symbols
      real(dp),         dimension(3, n_atoms),  intent(out) :: positions
      real(dp),         dimension(n_atoms),     intent(out) :: chi
      real(dp),         dimension(n_atoms),     intent(out) :: eta
      integer,          dimension(n_molecules), intent(out) :: n_atoms_per_molecule

      integer :: record, cursor, current_atom, i, current_molecule

      character(len=200) :: string, coordinate

      current_atom = 0

      current_molecule     = 0
      n_atoms_per_molecule = 0
      cursor = 0

      do record = 1, this%n_mm_atom_lines

         current_atom = current_atom + 1
         string = trim(adjustl(this%mm_geometry(record)))

         ! Determine symbol

         symbols(current_atom) = string(1:2)
         string = adjustl(string(3:200))

         ! Determine molecule index

         if (is_substring_in_string(string, 'mol')) then
            cursor = last_index_of_substring(string, 'mol')
         else
            call output%error_msg('could not find mol in MM geometry')
         endif

         string = adjustl(string(cursor+1:200))

         if (string(1:1) .ne. '=') call output%error_msg('in mol specification in MM geometry')
         string = adjustl(string(2:200))

         cursor = first_index_of_substring(string, ']')

         read(string(1:cursor-1),*) current_molecule
         string = adjustl(string(cursor+1:200))

         n_atoms_per_molecule(current_molecule) = n_atoms_per_molecule(current_molecule) + 1

         cursor = first_index_of_substring(string, '[')

         coordinate = string(1:cursor-1)
         read(coordinate, *) positions(:, current_atom)

         string = adjustl(string(cursor+1:200))

         ! Determine chi

         if (is_substring_in_string(string, 'chi')) then
            cursor = last_index_of_substring(string, 'chi')
         else
            call output%error_msg('could not find chi in MM geometry')
         endif

         string = adjustl(string(cursor+1:200))

         if (string(1:1) .ne. '=') call output%error_msg('in chi specification in MM geometry')
         string = adjustl(string(2:200))

         cursor = first_index_of_substring(string,',')
         read(string(1:cursor-1), * ) chi(current_atom)
         string = adjustl(string(cursor+1:200))

         ! Determine eta

         if (is_substring_in_string(string, 'eta')) then
            cursor = last_index_of_substring(string, 'eta')
         else
            call output%error_msg('could not find eta in MM geometry')
         endif

         string = trim(adjustl(string(cursor+1:200)))

         if (string(1:1) .ne. '=') call output%error_msg('in chi specification in MM geometry')
         string = adjustl(string(2:200))

         cursor = first_index_of_substring(string,']')
         read(string(1:cursor-1), * ) eta(current_atom)

         if(abs(eta(current_atom)).lt.1.0d-8) then

            call output%error_msg('You have put zero chemical hardness on atom: (i0)', &
                                  ints=[current_atom])

         endif

      enddo

      ! First character of symbol should be upper case

      do i = 1, n_atoms

         symbols(i)(1:1) = convert_character_to_uppercase(symbols(i)(1:1))

      enddo

   end subroutine get_mm_geometry_fq_input_tool


   subroutine get_mm_geometry_non_polarizable_input_tool(this, n_atoms, symbols,   &
                                                         positions, charge)
      !!
      !! Written by Tommaso Giovanini and Sarai D. Folkestad, 2019-2020
      !!
      !! Reads the geometry of the MM portion from the input file and
      !! sets it in the list of atoms.
      !!
      !! Note: In order to be run, you need to know the number of MM atoms
      !!       there are in the system
      !!
      use string_utilities, only: convert_character_to_uppercase
      use string_utilities, only: last_index_of_substring

      implicit none

      class(input_tool), intent(in) :: this

      integer, intent(in) :: n_atoms

      character(len=2), dimension(n_atoms),     intent(out) :: symbols
      real(dp),         dimension(3, n_atoms),  intent(out) :: positions
      real(dp),         dimension(n_atoms),     intent(out) :: charge

      integer :: record, cursor, current_atom, i, current_molecule

      character(len=200) :: string, coordinate

      current_atom      = 0
      current_molecule  = 0
      cursor            = 0

      do record = 1, this%n_mm_atom_lines

         current_atom = current_atom + 1
         string = trim(adjustl(this%mm_geometry(record)))

         ! Determine symbol

         symbols(current_atom) = string(1:2)
         string = adjustl(string(3:200))

         ! Determine molecule index

         if (is_substring_in_string(string, 'mol')) then
            cursor = last_index_of_substring(string, 'mol')
         else
            call output%error_msg('could not find mol in MM geometry')
         endif

         string = adjustl(string(cursor+1:200))

         if (string(1:1) .ne. '=') call output%error_msg('in mol specification in MM geometry')
         string = adjustl(string(2:200))

         cursor = first_index_of_substring(string, ']')
         read(string(1:cursor-1),*) current_molecule
         string = adjustl(string(cursor+1:200))

         ! Determine position

         cursor = first_index_of_substring(string, '[')

         coordinate = string(1:cursor-1)
         read(coordinate, *) positions(:, current_atom)

         string = adjustl(string(cursor+1:200))

         ! Determine charge

         cursor = first_index_of_substring(string, 'q')
         string = adjustl(string(cursor+1:200))

         if (string(1:1) .ne. '=') call output%error_msg('in charge specification in MM geometry')
         string = adjustl(string(2:200))

         cursor = first_index_of_substring(string,']')
         read(string(1:cursor-1), * ) charge(current_atom)

         if(abs(charge(current_atom)).lt.1.0d-8) then

            call output%warning_msg('You put zero charge on atom = (i0)', &
                                     ints=[current_atom], fs='(t6,a)')

         endif

      enddo

      ! First character of symbol should be upper case

      do i = 1, n_atoms

         symbols(i)(1:1) = convert_character_to_uppercase(symbols(i)(1:1))

      enddo

   end subroutine get_mm_geometry_non_polarizable_input_tool


   function string_in_array_keyword(this, keyword, section, string) result(in_list)
      !!
      !! Written by Sarai D. Folkestad and Alexander C. Paul, 2019-2023
      !!
      !! Determines if the 'string' is one of the keyword values in the comma separated
      !! list given with 'keyword' in the given 'section' of the input file.
      !!
      !! Examples of such a keyword is 'levels' from the 'mlcc' section:
      !!
      !!    levels: ccs, cc2, ccsd
      !!
      !! This function assumes that the keyword is requested. If it is not, a test should be performed
      !! before the routine is called
      !!
      use string_utilities, only: string_in_comma_separated_list

      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: section
      character(len=*), intent(in) :: keyword
      character(len=*), intent(in) :: string

      logical :: in_list

      character(len=200) :: keyword_value_string

      call this%get_required_keyword(keyword, section, keyword_value_string)
      in_list = string_in_comma_separated_list(keyword_value_string, string)

   end function string_in_array_keyword


   subroutine cleanup_input_tool(this)
      !!
      !! Written by Eirik F. Kjønstad, May 2022
      !!
      implicit none

      class(input_tool), intent(inout) :: this

      deallocate(this%keyword_lines)
      deallocate(this%qm_geometry)
      if (allocated(this%mm_geometry)) deallocate(this%mm_geometry)

   end subroutine cleanup_input_tool


   function is_embedding_on_input_tool(this) result(embedding)
      !!
      !! Written by Sarai D. Folkestad, Sep 2020
      !!
      !! Returns true if wavefunction is embedded.
      !!
      implicit none

      class(input_tool), intent(in) :: this

      logical :: embedding

      embedding = this%is_section_present('molecular mechanics') .or. &
                  this%is_section_present('pcm') .or. &
                  this%is_keyword_present('cbo', 'boson')

   end function is_embedding_on_input_tool


   subroutine place_records_in_memory_input_tool(this, section, records_in_memory)
      !!
      !! Written by Sarai D. Folkestad, 2021
      !!
      !! Checks if storage for 'section' is in memory or on disk
      !!
      !! records_in_mem is intent(inout) because it should be set
      !! to a default value before a call to this routine
      !!
      implicit none

      class(input_tool), intent(in) :: this
      character(len=*),  intent(in) :: section

      logical, intent(inout) :: records_in_memory

      character(len=200) :: storage

      if (.not. this%is_keyword_present('storage', trim(section))) return

      call this%get_keyword('storage', trim(section), storage)

      ! Determine whether to store records in memory or on file

      if (trim(storage) == 'memory') then

         records_in_memory = .true.

      elseif (trim(storage) == 'disk') then

         records_in_memory = .false.

      else

         call output%error_msg('Could not recognize keyword storage in solver: ', &
                                 chars=[trim(storage)])

      endif

   end subroutine place_records_in_memory_input_tool


   function get_n_sets_and_ranges(string) result(n_elements)
      !!
      !! Written by Sarai D. Folkstad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Gets the number of sets/ranges in a string,
      !!
      !! Ranges should always be given as [a,b].
      !!
      !! Lists/ranges should always be given as {a, b, c, d},
      !! that is, in set notation.
      !!
      implicit none

      character(len=200), intent(inout) :: string

      integer :: n_elements

      integer :: i
      integer :: n_set_start, n_set_end
      integer :: n_range_start, n_range_end

      n_set_start = 0
      n_set_end   = 0

      n_range_start = 0
      n_range_end   = 0

      string = adjustl(string)

      do i = 1, len_trim(string)

         if (string(i:i) == '{') n_set_start = n_set_start + 1
         if (string(i:i) == '}') n_set_end = n_set_end + 1

         if (string(i:i) == '[') n_range_start = n_range_start + 1
         if (string(i:i) == ']') n_range_end = n_range_end + 1

      enddo

      if (n_set_start /= n_set_end) call output%error_msg('found open set in input file')
      if (n_range_start /= n_range_end) call output%error_msg('found open range in input file')

      n_elements = n_set_start + n_range_start

   end function get_n_sets_and_ranges


   function get_n_elements_in_set_or_range(this, string) result(n_elements)
      !!
      !! Written by Sarai D. Folkstad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Gets the number of elements in range or list,
      !!
      !! Ranges should always be given as [a,b], where a and b are integers.
      !! Lists/Sets should always be given as {a, b, c, d}.
      !!
      use string_utilities, only: n_instances_of_character

      implicit none

      class(input_tool), intent(in) :: this

      character(len=200), intent(inout) :: string

      integer :: n_elements

      integer :: last, first, n_characters
      integer :: comma_index

      n_elements = 0

      string = adjustl(string)

      n_characters = len_trim(string)

      if (string_starts_with(string, '[')) then ! range

         call expect_range(string)

         comma_index = index(string, ',')

         read(string(2 : comma_index-1), *) first
         read(string(comma_index+1 : n_characters - 1), *) last

         n_elements = last - first + 1

      else if (string_starts_with(string, '{')) then ! set

         call this%expect_set(string)

         n_elements = 1 ! Assuming that the set contains at least one element (otherwise why give list?)
         n_elements = n_elements + n_instances_of_character(string(2:), ',')

      else ! Did not find list or range

         n_elements = 0

      endif

   end function get_n_elements_in_set_or_range


   subroutine get_integers_from_set_or_range(this, string, n_elements, elements)
      !!
      !! Written by Sarai D. Folkstad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Ranges should always be given as [a,b].
      !!
      !! Lists should always be given as {a, b, c, d}.
      !!
      use string_utilities, only: read_comma_separated_list

      implicit none

      class(input_tool), intent(in) :: this

      character(len=200), intent(inout) :: string

      integer, intent(in) :: n_elements

      integer, dimension(n_elements), intent(out) :: elements

      integer :: first, last, length, n_characters
      integer :: i, comma_index

      string = adjustl(string)

      n_characters = len_trim(string)

      if (string_starts_with(string, '[')) then ! range

         call expect_range(string)

         comma_index = index(string, ',')

         read(string(2 : comma_index-1), *) first
         read(string(comma_index+1 : n_characters - 1), *) last
         length = last - first + 1

         if (n_elements /= length) call output%error_msg('Mismatch in number of elements of range.')

         elements = [(i, i=first, last)]

      else if (string_starts_with(string, '{')) then ! set

         call this%expect_set(string)
         call read_comma_separated_list(string(2:n_characters-1), elements, n_elements)

      else ! Did not find list or range

         call output%error_msg('neither list nor range was found.')

      endif

   end subroutine get_integers_from_set_or_range


   function get_n_state_guesses(this) result(n_guesses)
      !!
      !! Written by Eirik F. Kjønstad, Dec 2021
      !!
      implicit none

      class(input_tool) :: this

      integer :: n_guesses

      character(len=200) :: state_guesses

      call this%get_keyword('state guesses', 'solver cc es', state_guesses)

      n_guesses = get_n_sets_and_ranges(state_guesses)

   end function get_n_state_guesses


   subroutine get_state_guesses(this, occupied, virtual, n_states)
      !!
      !! Written by Sarai D. Folkestad, Dec 2021
      !!
      !! state guesses: {i=1, a=2}, {i=2, a=2}, {i=1, a=1}
      !!
      !! Provided as a (comma separated) list of sets.
      !!
      !! A set comprises of two elements which are
      !! the occupied and virtual indices of the start vector
      !!
      !! 'i =' provides occupied index
      !! 'a =' provides virtual index
      !!
      use string_utilities, only: remove_delimiter_from_string, split_at_delimiter
      use string_utilities, only: n_instances_of_character

      implicit none

      class(input_tool) :: this

      integer, intent(in) :: n_states
      integer, dimension(n_states), intent(out) :: virtual, occupied

      character(len=200) :: state_guesses
      character(len=200), dimension(:), allocatable :: split_string

      integer :: n_elements, cursor_a, cursor_i, i, n_occupied, n_virtual

      call this%get_keyword('state guesses', 'solver cc es', state_guesses)
      state_guesses = adjustl(state_guesses)

      n_occupied = n_instances_of_character(state_guesses, 'i')
      n_virtual = n_instances_of_character(state_guesses, 'a')

      if (n_virtual /= n_occupied) &
         call output%error_msg('in keyword (state guesses)')

      if (n_virtual /= n_states) &
         call output%error_msg('in keyword (state guesses)')

      ! {i=1, a=2}, {i=2, a=2}, {i=1, a=1} -> i=1, a=2} i=2, a=2} i=1, a=1}
      call remove_delimiter_from_string(state_guesses, '{')

      n_elements = n_instances_of_character(state_guesses, '}')

      if (n_elements /= n_states) &
         call output%error_msg('in keyword (state guesses)')

      allocate(split_string(n_elements))

      ! i=1, a=2} i=2, a=2} i=1, a=1} -> string array ("i=1, a=2", "i=2, a=2", "i=1, a=1")
      call split_at_delimiter(state_guesses, n_elements, split_string, '}')

      do i = 1, n_elements

         call remove_delimiter_from_string(split_string(i), '=')
         call remove_delimiter_from_string(split_string(i), ',')

         cursor_i = first_index_of_substring(split_string(i), 'i')
         cursor_a = first_index_of_substring(split_string(i), 'a')

         if (cursor_i < cursor_a) then

            read(split_string(i)(cursor_i + 1 : cursor_a -1), *) occupied(i)
            read(split_string(i)(cursor_a + 1:), *) virtual(i)

         else

            read(split_string(i)(cursor_a + 1 : cursor_i -1), *) virtual(i)
            read(split_string(i)(cursor_i + 1:), *) occupied(i)

         endif

      enddo

      deallocate(split_string)

   end subroutine get_state_guesses


   subroutine check_for_errors(this)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, Mar 2019
      !!
      !! Looks for errors the user may have made when writing the input file.
      !!
      implicit none

      class(input_tool), intent(in) :: this
      integer :: k

      do k = 1, size(this%sections)

         call this%sections(k)%check_section_for_illegal_keywords()

      enddo

   end subroutine check_for_errors


   function is_section_recognized(this, section) result(allowed)
      !!
      !! Written by Sarai D. Folkestad and Eirik F. Kjønstad, 2019-2021
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: section

      logical :: allowed

      integer :: k

      allowed = .false.

      do k = 1, size(this%sections)

         if (this%sections(k)%name_ == trim(section)) then
            allowed = .true.
            return
         end if

      enddo

   end function is_section_recognized


   subroutine section_not_recognized(this, section)
      !!
      !! Written by Sara Angelico and Alexander C. Paul, 2023
      !!
      implicit none

      class(input_tool), intent(in) :: this

      character(len=*), intent(in) :: section

      call output%printf('m', 'Could not recognize section named "(a0)".', &
                        chars=[trim(section)], fs='(/t3,a)')
      call this%print_sections()
      call output%error_msg('Something is wrong in the input file. See above.')

   end subroutine section_not_recognized


   subroutine expect_set(string)
      !!
      !! Written by Alexander C. Paul, Mar 2023
      !!
      implicit none

      character(len=*), intent(in) :: string
      integer :: opening, closing

      opening = index(string, '{')
      if (opening == 0) call output%error_msg('Expected set, but no "{" found in string.')

      closing = index(string, '}')
      if (closing == 0) call output%error_msg('Expected set, but no "}" found in string.')

      if (closing < opening) call output%error_msg('Expected set, but "}" found before "{".')

   end subroutine expect_set


   subroutine expect_range(string)
      !!
      !! Written by Alexander C. Paul, Mar 2023
      !!
      implicit none

      character(len=*), intent(in) :: string
      integer :: opening, closing, comma

      opening = index(string, '[')
      if (opening == 0) call output%error_msg('Expected range, but no "[" found in string.')

      closing = index(string, ']')
      if (closing == 0) call output%error_msg('Expected range, but no "]" found in string.')

      comma = index(string, ',')
      if (comma == 0) call output%error_msg('Expected range, but no "," found in string.')

      if (comma < opening) call output%error_msg('Expected range, but "," found before "[".')
      if (closing < comma) call output%error_msg('Expected range, but "]" found before ",".')

   end subroutine expect_range


   subroutine set_section_keywords(this, section, section_lines)
      !!
      !! Written by Sara Angelico and Alexander C. Paul, Apr 2023
      !!
      !! This routine removes empty lines and comment lines and
      !! sets the keywords for the section
      !!
      implicit none

      class(input_tool),                intent(inout) :: this
      character(len=*),                 intent(in)    :: section
      character(len=200), dimension(:), intent(in)    :: section_lines

      character(len=200), dimension(:), allocatable :: keyword_lines
      character(len=200) :: line
      integer :: section_index, n_lines, l

      allocate(keyword_lines, mold=section_lines)

      n_lines = 0
      do l = 1, size(section_lines)

         line = trim(adjustl(section_lines(l)))

         if (string_starts_with(line, "!")) cycle
         if (is_empty(line)) cycle

         line = remove_substring_after_symbol(line,'!')

         n_lines = n_lines + 1
         keyword_lines(n_lines) = adjustl(line)

      end do

      section_index = this%get_section_index(trim(section))

      call this%sections(section_index)%set_keyword_lines(keyword_lines(1:n_lines))

      deallocate(keyword_lines)

   end subroutine set_section_keywords


   subroutine check_get_array_for_keyword(keyword, section, n_elements)
      !!
      !! Written by Alexander C. Paul, May 2023
      !!
      implicit none

      character(len=*), intent(in) :: keyword, section
      integer,          intent(in) :: n_elements

      if (n_elements .le. 0) then
         call output%error_msg("Get array for keyword called with n_elements=(i0). &
                              &Keyword: (a0) Section: " // section, &
                               ints=[n_elements], chars=[keyword])
      end if

   end subroutine check_get_array_for_keyword


end module input_tool_class

eT-program / eT / 23895

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous