LCOV - code coverage report
Current view: top level - src - atom_basis.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 132 135 97.8 %
Date: 2024-03-29 07:50:05 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : MODULE atom_basis
       9             :    USE atom_fit,                        ONLY: atom_fit_basis
      10             :    USE atom_output,                     ONLY: atom_print_basis,&
      11             :                                               atom_print_info,&
      12             :                                               atom_print_method,&
      13             :                                               atom_print_potential
      14             :    USE atom_types,                      ONLY: &
      15             :         atom_basis_type, atom_integrals, atom_optimization_type, atom_orbitals, atom_p_type, &
      16             :         atom_potential_type, atom_state, create_atom_orbs, create_atom_type, init_atom_basis, &
      17             :         init_atom_potential, lmat, read_atom_opt_section, release_atom_basis, &
      18             :         release_atom_potential, release_atom_type, set_atom
      19             :    USE atom_utils,                      ONLY: atom_consistent_method,&
      20             :                                               atom_set_occupation,&
      21             :                                               get_maxl_occ,&
      22             :                                               get_maxn_occ
      23             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24             :                                               cp_logger_type
      25             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      26             :                                               cp_print_key_unit_nr
      27             :    USE input_constants,                 ONLY: do_analytic
      28             :    USE input_section_types,             ONLY: section_vals_get,&
      29             :                                               section_vals_get_subs_vals,&
      30             :                                               section_vals_type,&
      31             :                                               section_vals_val_get
      32             :    USE kinds,                           ONLY: default_string_length,&
      33             :                                               dp
      34             :    USE periodic_table,                  ONLY: nelem,&
      35             :                                               ptable
      36             : #include "./base/base_uses.f90"
      37             : 
      38             :    IMPLICIT NONE
      39             :    PRIVATE
      40             :    PUBLIC  :: atom_basis_opt
      41             : 
      42             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_basis'
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief Optimize the atomic basis set.
      48             : !> \param atom_section  ATOM input section
      49             : !> \par History
      50             : !>    * 04.2009 created starting from the subroutine atom_energy() [Juerg Hutter]
      51             : ! **************************************************************************************************
      52          10 :    SUBROUTINE atom_basis_opt(atom_section)
      53             :       TYPE(section_vals_type), POINTER                   :: atom_section
      54             : 
      55             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'atom_basis_opt'
      56             : 
      57             :       CHARACTER(LEN=2)                                   :: elem
      58             :       CHARACTER(LEN=default_string_length), &
      59          10 :          DIMENSION(:), POINTER                           :: tmpstringlist
      60             :       INTEGER                                            :: do_eric, do_erie, handle, i, im, in, &
      61             :                                                             iunit, iw, k, maxl, mb, method, mo, &
      62             :                                                             n_meth, n_rep, nr_gh, reltyp, zcore, &
      63             :                                                             zval, zz
      64             :       INTEGER, DIMENSION(0:lmat)                         :: maxn
      65          10 :       INTEGER, DIMENSION(:), POINTER                     :: cn
      66             :       LOGICAL                                            :: do_gh, eri_c, eri_e, had_ae, had_pp, &
      67             :                                                             pp_calc
      68             :       REAL(KIND=dp), DIMENSION(0:lmat, 10)               :: pocc
      69             :       TYPE(atom_basis_type), POINTER                     :: ae_basis, pp_basis
      70             :       TYPE(atom_integrals), POINTER                      :: ae_int, pp_int
      71             :       TYPE(atom_optimization_type)                       :: optimization
      72             :       TYPE(atom_orbitals), POINTER                       :: orbitals
      73          10 :       TYPE(atom_p_type), DIMENSION(:, :), POINTER        :: atom_info
      74             :       TYPE(atom_potential_type), POINTER                 :: ae_pot, p_pot
      75             :       TYPE(atom_state), POINTER                          :: state
      76             :       TYPE(cp_logger_type), POINTER                      :: logger
      77             :       TYPE(section_vals_type), POINTER                   :: basis_section, method_section, &
      78             :                                                             opt_section, potential_section, &
      79             :                                                             powell_section, xc_section
      80             : 
      81          10 :       CALL timeset(routineN, handle)
      82             : 
      83             :       ! What atom do we calculate
      84          10 :       CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
      85          10 :       CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
      86          10 :       zz = 0
      87         222 :       DO i = 1, nelem
      88         222 :          IF (ptable(i)%symbol == elem) THEN
      89             :             zz = i
      90             :             EXIT
      91             :          END IF
      92             :       END DO
      93          10 :       IF (zz /= 1) zval = zz
      94             : 
      95             :       ! read and set up inofrmation on the basis sets
      96         370 :       ALLOCATE (ae_basis, pp_basis)
      97          10 :       basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
      98          10 :       NULLIFY (ae_basis%grid)
      99          10 :       CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
     100          10 :       NULLIFY (pp_basis%grid)
     101          10 :       basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
     102          10 :       CALL init_atom_basis(pp_basis, basis_section, zval, "PP")
     103             : 
     104             :       ! print general and basis set information
     105          10 :       logger => cp_get_default_logger()
     106          10 :       iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
     107          10 :       IF (iw > 0) CALL atom_print_info(zval, "Atomic Basis Optimization", iw)
     108          10 :       CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
     109             : 
     110             :       ! read and setup information on the pseudopotential
     111          10 :       NULLIFY (potential_section)
     112          10 :       potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
     113      105050 :       ALLOCATE (ae_pot, p_pot)
     114          10 :       CALL init_atom_potential(p_pot, potential_section, zval)
     115          10 :       CALL init_atom_potential(ae_pot, potential_section, -1)
     116             : 
     117             :       ! if the ERI's are calculated analytically, we have to precalculate them
     118          10 :       eri_c = .FALSE.
     119          10 :       CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
     120          10 :       IF (do_eric == do_analytic) eri_c = .TRUE.
     121          10 :       eri_e = .FALSE.
     122          10 :       CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
     123          10 :       IF (do_erie == do_analytic) eri_e = .TRUE.
     124          10 :       CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
     125          10 :       CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
     126             : 
     127             :       ! information on the states to be calculated
     128          10 :       CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
     129          10 :       maxn = 0
     130          10 :       CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
     131          20 :       DO in = 1, MIN(SIZE(cn), 4)
     132          20 :          maxn(in - 1) = cn(in)
     133             :       END DO
     134          70 :       DO in = 0, lmat
     135          60 :          maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
     136          70 :          maxn(in) = MIN(maxn(in), pp_basis%nbas(in))
     137             :       END DO
     138             : 
     139             :       ! read optimization section
     140          10 :       opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
     141          10 :       CALL read_atom_opt_section(optimization, opt_section)
     142             : 
     143          10 :       had_ae = .FALSE.
     144          10 :       had_pp = .FALSE.
     145             : 
     146             :       ! Check for the total number of electron configurations to be calculated
     147          10 :       CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
     148             :       ! Check for the total number of method types to be calculated
     149          10 :       method_section => section_vals_get_subs_vals(atom_section, "METHOD")
     150          10 :       CALL section_vals_get(method_section, n_repetition=n_meth)
     151             : 
     152             :       ! integrals
     153        4250 :       ALLOCATE (ae_int, pp_int)
     154             : 
     155          60 :       ALLOCATE (atom_info(n_rep, n_meth))
     156             : 
     157          20 :       DO in = 1, n_rep
     158          30 :          DO im = 1, n_meth
     159             : 
     160          10 :             NULLIFY (atom_info(in, im)%atom)
     161          10 :             CALL create_atom_type(atom_info(in, im)%atom)
     162             : 
     163          10 :             atom_info(in, im)%atom%optimization = optimization
     164             : 
     165          10 :             atom_info(in, im)%atom%z = zval
     166          10 :             xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
     167          10 :             atom_info(in, im)%atom%xc_section => xc_section
     168             : 
     169        3630 :             ALLOCATE (state)
     170             : 
     171             :             ! get the electronic configuration
     172             :             CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
     173          10 :                                       c_vals=tmpstringlist)
     174             : 
     175             :             ! set occupations
     176          10 :             CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
     177          10 :             state%maxl_occ = get_maxl_occ(state%occ)
     178          70 :             state%maxn_occ = get_maxn_occ(state%occ)
     179             : 
     180             :             ! set number of states to be calculated
     181          10 :             state%maxl_calc = MAX(maxl, state%maxl_occ)
     182          10 :             state%maxl_calc = MIN(lmat, state%maxl_calc)
     183          70 :             state%maxn_calc = 0
     184          30 :             DO k = 0, state%maxl_calc
     185          30 :                state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k))
     186             :             END DO
     187             : 
     188             :             ! is there a pseudo potential
     189          22 :             pp_calc = ANY(INDEX(tmpstringlist(1:), "CORE") /= 0)
     190          10 :             IF (pp_calc) THEN
     191             :                ! get and set the core occupations
     192           6 :                CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
     193           6 :                CALL atom_set_occupation(tmpstringlist, state%core, pocc)
     194         426 :                zcore = zval - NINT(SUM(state%core))
     195           6 :                CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
     196           6 :                had_pp = .TRUE.
     197           6 :                CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, potential=p_pot)
     198          78 :                state%maxn_calc(:) = MIN(state%maxn_calc(:), pp_basis%nbas(:))
     199          42 :                CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
     200             :             ELSE
     201         284 :                state%core = 0._dp
     202           4 :                CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
     203           4 :                had_ae = .TRUE.
     204           4 :                CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, potential=ae_pot)
     205          52 :                state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
     206          28 :                CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
     207             :             END IF
     208             : 
     209          10 :             CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_val=im)
     210          10 :             CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
     211          10 :             CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)
     212          10 :             CALL set_atom(atom_info(in, im)%atom, state=state)
     213             :             CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
     214          10 :                           exchange_integral_type=do_erie)
     215          10 :             atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
     216          10 :             atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
     217             : 
     218          10 :             IF (atom_consistent_method(method, state%multiplicity)) THEN
     219          10 :                iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
     220          10 :                CALL atom_print_method(atom_info(in, im)%atom, iw)
     221          10 :                CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
     222          10 :                iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
     223          10 :                IF (pp_calc) THEN
     224           6 :                   IF (iw > 0) CALL atom_print_potential(p_pot, iw)
     225             :                ELSE
     226           4 :                   IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
     227             :                END IF
     228          10 :                CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
     229             :             ELSE
     230           0 :                CPABORT("METHOD_TYPE and MULTIPLICITY are incompatible")
     231             :             END IF
     232             : 
     233          10 :             NULLIFY (orbitals)
     234          70 :             mo = MAXVAL(state%maxn_calc)
     235          70 :             mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
     236          10 :             CALL create_atom_orbs(orbitals, mb, mo)
     237          30 :             CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
     238             : 
     239             :          END DO
     240             :       END DO
     241             : 
     242             :       ! Start the Optimization
     243          10 :       powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
     244          10 :       iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
     245          10 :       iunit = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_BASIS", extension=".log")
     246          10 :       IF (had_ae) THEN
     247           4 :          pp_calc = .FALSE.
     248           4 :          CALL atom_fit_basis(atom_info, ae_basis, pp_calc, iunit, powell_section)
     249             :       END IF
     250          10 :       IF (had_pp) THEN
     251           6 :          pp_calc = .TRUE.
     252           6 :          CALL atom_fit_basis(atom_info, pp_basis, pp_calc, iunit, powell_section)
     253             :       END IF
     254          10 :       CALL cp_print_key_finished_output(iunit, logger, atom_section, "PRINT%FIT_BASIS")
     255          10 :       CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
     256          10 :       iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
     257          10 :       IF (iw > 0) THEN
     258           0 :          CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
     259           0 :          CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
     260             :       END IF
     261          10 :       CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
     262             : 
     263          10 :       CALL release_atom_basis(ae_basis)
     264          10 :       CALL release_atom_basis(pp_basis)
     265             : 
     266          10 :       CALL release_atom_potential(p_pot)
     267          10 :       CALL release_atom_potential(ae_pot)
     268             : 
     269          20 :       DO in = 1, n_rep
     270          30 :          DO im = 1, n_meth
     271          20 :             CALL release_atom_type(atom_info(in, im)%atom)
     272             :          END DO
     273             :       END DO
     274          10 :       DEALLOCATE (atom_info)
     275             : 
     276          10 :       DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
     277             : 
     278          10 :       CALL timestop(handle)
     279             : 
     280          40 :    END SUBROUTINE atom_basis_opt
     281             : 
     282             : END MODULE atom_basis

Generated by: LCOV version 1.15