LCOV - code coverage report
Current view: top level - src - atom_grb.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 689 702 98.1 %
Date: 2024-03-29 07:50:05 Functions: 8 9 88.9 %

          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_grb
       9             :    USE ai_onecenter,                    ONLY: sg_conf,&
      10             :                                               sg_kinetic,&
      11             :                                               sg_nuclear,&
      12             :                                               sg_overlap
      13             :    USE atom_electronic_structure,       ONLY: calculate_atom
      14             :    USE atom_operators,                  ONLY: atom_int_release,&
      15             :                                               atom_int_setup,&
      16             :                                               atom_ppint_release,&
      17             :                                               atom_ppint_setup,&
      18             :                                               atom_relint_release,&
      19             :                                               atom_relint_setup
      20             :    USE atom_types,                      ONLY: &
      21             :         CGTO_BASIS, GTO_BASIS, atom_basis_type, atom_integrals, atom_orbitals, atom_p_type, &
      22             :         atom_potential_type, atom_state, atom_type, create_atom_orbs, create_atom_type, lmat, &
      23             :         release_atom_basis, release_atom_type, set_atom
      24             :    USE atom_utils,                      ONLY: atom_basis_condnum,&
      25             :                                               atom_density
      26             :    USE cp_files,                        ONLY: close_file,&
      27             :                                               open_file
      28             :    USE input_constants,                 ONLY: barrier_conf,&
      29             :                                               do_analytic,&
      30             :                                               do_rhf_atom,&
      31             :                                               do_rks_atom,&
      32             :                                               do_rohf_atom,&
      33             :                                               do_uhf_atom,&
      34             :                                               do_uks_atom
      35             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      36             :                                               section_vals_type,&
      37             :                                               section_vals_val_get
      38             :    USE kinds,                           ONLY: default_string_length,&
      39             :                                               dp
      40             :    USE lapack,                          ONLY: lapack_ssygv
      41             :    USE mathconstants,                   ONLY: dfac,&
      42             :                                               pi
      43             :    USE orbital_pointers,                ONLY: deallocate_orbital_pointers,&
      44             :                                               init_orbital_pointers
      45             :    USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
      46             :                                               init_spherical_harmonics
      47             :    USE periodic_table,                  ONLY: ptable
      48             :    USE physcon,                         ONLY: bohr
      49             :    USE powell,                          ONLY: opt_state_type,&
      50             :                                               powell_optimize
      51             :    USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
      52             :                                               create_grid_atom
      53             : #include "./base/base_uses.f90"
      54             : 
      55             :    IMPLICIT NONE
      56             : 
      57             :    TYPE basis_p_type
      58             :       TYPE(atom_basis_type), POINTER                :: basis => NULL()
      59             :    END TYPE basis_p_type
      60             : 
      61             :    PRIVATE
      62             :    PUBLIC  :: atom_grb_construction
      63             : 
      64             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_grb'
      65             : 
      66             : CONTAINS
      67             : 
      68             : ! **************************************************************************************************
      69             : !> \brief Construct geometrical response basis set.
      70             : !> \param atom_info    information about the atomic kind. Two-dimensional array of size
      71             : !>                     (electronic-configuration, electronic-structure-method)
      72             : !> \param atom_section ATOM input section
      73             : !> \param iw           output file unit
      74             : !> \par History
      75             : !>    * 11.2016 created [Juerg Hutter]
      76             : ! **************************************************************************************************
      77           2 :    SUBROUTINE atom_grb_construction(atom_info, atom_section, iw)
      78             : 
      79             :       TYPE(atom_p_type), DIMENSION(:, :), POINTER        :: atom_info
      80             :       TYPE(section_vals_type), POINTER                   :: atom_section
      81             :       INTEGER, INTENT(IN)                                :: iw
      82             : 
      83             :       CHARACTER(len=default_string_length)               :: abas, basname
      84             :       CHARACTER(len=default_string_length), DIMENSION(1) :: basline
      85             :       CHARACTER(len=default_string_length), DIMENSION(3) :: headline
      86             :       INTEGER                                            :: i, ider, is, iunit, j, k, l, lhomo, ll, &
      87             :                                                             lval, m, maxl, mb, method, mo, n, &
      88             :                                                             nder, ngp, nhomo, nr, num_gto, &
      89             :                                                             num_pol, quadtype, s1, s2
      90             :       INTEGER, DIMENSION(0:7)                            :: nbas
      91             :       INTEGER, DIMENSION(0:lmat)                         :: next_bas, next_prim
      92           2 :       INTEGER, DIMENSION(:), POINTER                     :: num_bas
      93             :       REAL(KIND=dp) :: al, amin, aval, cnum, crad, cradx, cval, delta, dene, ear, emax, &
      94             :          energy_ex(0:lmat), energy_ref, energy_vb(0:lmat), expzet, fhomo, o, prefac, rconf, rk, &
      95             :          rmax, scon, zeta, zval
      96           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ale, alp, rho
      97           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat
      98           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ebasis, pbasis, qbasis, rbasis
      99           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: wfn
     100           2 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: ovlp
     101             :       TYPE(atom_basis_type), POINTER                     :: basis, basis_grb, basis_ref, basis_vrb
     102             :       TYPE(atom_integrals), POINTER                      :: atint
     103             :       TYPE(atom_orbitals), POINTER                       :: orbitals
     104             :       TYPE(atom_state), POINTER                          :: state
     105             :       TYPE(atom_type), POINTER                           :: atom, atom_ref, atom_test
     106          24 :       TYPE(basis_p_type), DIMENSION(0:10)                :: vbasis
     107             :       TYPE(section_vals_type), POINTER                   :: grb_section, powell_section
     108             : 
     109           2 :       IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T28,A,/," ",79("*"))') "GEOMETRICAL RESPONSE BASIS"
     110             : 
     111          24 :       DO i = 0, 10
     112          24 :          NULLIFY (vbasis(i)%basis)
     113             :       END DO
     114             :       ! make some basic checks
     115           2 :       is = SIZE(atom_info)
     116           2 :       IF (iw > 0 .AND. is > 1) THEN
     117           0 :          WRITE (iw, '(/,A,/)') " WARNING: Only use first electronic structure/method for basis set generation"
     118             :       END IF
     119           2 :       atom_ref => atom_info(1, 1)%atom
     120             : 
     121             :       ! check method
     122           2 :       method = atom_ref%method_type
     123           0 :       SELECT CASE (method)
     124             :       CASE (do_rks_atom, do_rhf_atom)
     125             :          ! restricted methods are okay
     126             :       CASE (do_uks_atom, do_uhf_atom, do_rohf_atom)
     127           0 :          CPABORT("Unrestricted methods not allowed for GRB generation")
     128             :       CASE DEFAULT
     129           2 :          CPABORT("")
     130             :       END SELECT
     131             : 
     132             :       ! input for basis optimization
     133           2 :       grb_section => section_vals_get_subs_vals(atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
     134             : 
     135             :       ! generate an atom type
     136           2 :       NULLIFY (atom)
     137           2 :       CALL create_atom_type(atom)
     138           2 :       CALL copy_atom_basics(atom_ref, atom, state=.TRUE., potential=.TRUE., optimization=.TRUE., xc=.TRUE.)
     139             :       ! set confinement potential
     140           2 :       atom%potential%confinement = .TRUE.
     141           2 :       atom%potential%conf_type = barrier_conf
     142           2 :       atom%potential%acon = 200._dp
     143           2 :       atom%potential%rcon = 4._dp
     144           2 :       CALL section_vals_val_get(grb_section, "CONFINEMENT", r_val=scon)
     145           2 :       atom%potential%scon = scon
     146             :       ! generate main block geometrical exponents
     147           2 :       basis_ref => atom_ref%basis
     148          38 :       ALLOCATE (basis)
     149             :       NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
     150             :       ! get information on quadrature type and number of grid points
     151             :       ! allocate and initialize the atomic grid
     152             :       NULLIFY (basis%grid)
     153           2 :       CALL allocate_grid_atom(basis%grid)
     154           2 :       CALL section_vals_val_get(grb_section, "QUADRATURE", i_val=quadtype)
     155           2 :       CALL section_vals_val_get(grb_section, "GRID_POINTS", i_val=ngp)
     156           2 :       IF (ngp <= 0) &
     157           0 :          CPABORT("# point radial grid < 0")
     158           2 :       CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
     159           2 :       basis%grid%nr = ngp
     160             :       !
     161           2 :       maxl = atom%state%maxl_occ
     162           2 :       basis%basis_type = GTO_BASIS
     163           2 :       CALL section_vals_val_get(grb_section, "NUM_GTO_CORE", i_val=num_gto)
     164          14 :       basis%nbas = 0
     165           6 :       basis%nbas(0:maxl) = num_gto
     166          14 :       basis%nprim = basis%nbas
     167           2 :       CALL section_vals_val_get(grb_section, "GEOMETRICAL_FACTOR", r_val=cval)
     168           2 :       CALL section_vals_val_get(grb_section, "GEO_START_VALUE", r_val=aval)
     169          14 :       m = MAXVAL(basis%nbas)
     170           6 :       ALLOCATE (basis%am(m, 0:lmat))
     171          86 :       basis%am = 0._dp
     172          14 :       DO l = 0, lmat
     173          38 :          DO i = 1, basis%nbas(l)
     174          24 :             ll = i - 1
     175          36 :             basis%am(i, l) = aval*cval**(ll)
     176             :          END DO
     177             :       END DO
     178             : 
     179           2 :       basis%eps_eig = basis_ref%eps_eig
     180           2 :       basis%geometrical = .TRUE.
     181           2 :       basis%aval = aval
     182           2 :       basis%cval = cval
     183          14 :       basis%start = 0
     184             : 
     185             :       ! initialize basis function on a radial grid
     186           2 :       nr = basis%grid%nr
     187          14 :       m = MAXVAL(basis%nbas)
     188          10 :       ALLOCATE (basis%bf(nr, m, 0:lmat))
     189          10 :       ALLOCATE (basis%dbf(nr, m, 0:lmat))
     190          10 :       ALLOCATE (basis%ddbf(nr, m, 0:lmat))
     191       28886 :       basis%bf = 0._dp
     192       28886 :       basis%dbf = 0._dp
     193       28886 :       basis%ddbf = 0._dp
     194          14 :       DO l = 0, lmat
     195          38 :          DO i = 1, basis%nbas(l)
     196          24 :             al = basis%am(i, l)
     197        9636 :             DO k = 1, nr
     198        9600 :                rk = basis%grid%rad(k)
     199        9600 :                ear = EXP(-al*basis%grid%rad(k)**2)
     200        9600 :                basis%bf(k, i, l) = rk**l*ear
     201        9600 :                basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
     202             :                basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
     203        9624 :                                       2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
     204             :             END DO
     205             :          END DO
     206             :       END DO
     207             : 
     208           2 :       NULLIFY (orbitals)
     209          14 :       mo = MAXVAL(atom%state%maxn_calc)
     210          14 :       mb = MAXVAL(basis%nbas)
     211           2 :       CALL create_atom_orbs(orbitals, mb, mo)
     212           2 :       CALL set_atom(atom, orbitals=orbitals)
     213             : 
     214           2 :       powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
     215           2 :       CALL atom_fit_grb(atom, basis, iw, powell_section)
     216           2 :       CALL set_atom(atom, basis=basis)
     217             : 
     218             :       ! generate response contractions
     219           2 :       CALL section_vals_val_get(grb_section, "DELTA_CHARGE", r_val=delta)
     220           2 :       CALL section_vals_val_get(grb_section, "DERIVATIVES", i_val=nder)
     221           2 :       IF (iw > 0) THEN
     222           2 :          WRITE (iw, '(/,A,T76,I5)') " Generate Response Basis Sets with Order ", nder
     223             :       END IF
     224             : 
     225           2 :       state => atom%state
     226             :       ! find HOMO
     227           2 :       lhomo = -1
     228           2 :       nhomo = -1
     229           2 :       emax = -HUGE(1._dp)
     230           6 :       DO l = 0, state%maxl_occ
     231          10 :          DO i = 1, state%maxn_occ(l)
     232           8 :             IF (atom%orbitals%ener(i, l) > emax) THEN
     233           4 :                lhomo = l
     234           4 :                nhomo = i
     235           4 :                emax = atom%orbitals%ener(i, l)
     236           4 :                fhomo = state%occupation(l, i)
     237             :             END IF
     238             :          END DO
     239             :       END DO
     240             : 
     241           2 :       s1 = SIZE(atom%orbitals%wfn, 1)
     242           2 :       s2 = SIZE(atom%orbitals%wfn, 2)
     243          12 :       ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
     244          14 :       s2 = MAXVAL(state%maxn_occ) + nder
     245          18 :       ALLOCATE (rbasis(s1, s2, 0:lmat), qbasis(s1, s2, 0:lmat))
     246         350 :       rbasis = 0._dp
     247         350 :       qbasis = 0._dp
     248             : 
     249             :       ! calculate integrals
     250         426 :       ALLOCATE (atint)
     251           2 :       CALL atom_int_setup(atint, basis, potential=atom%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
     252           2 :       CALL atom_ppint_setup(atint, basis, potential=atom%potential)
     253           2 :       IF (atom%pp_calc) THEN
     254           2 :          NULLIFY (atint%tzora, atint%hdkh)
     255             :       ELSE
     256             :          ! relativistic correction terms
     257           0 :          CALL atom_relint_setup(atint, basis, atom%relativistic, zcore=REAL(atom%z, dp))
     258             :       END IF
     259           2 :       CALL set_atom(atom, integrals=atint)
     260             : 
     261           2 :       CALL calculate_atom(atom, iw=0)
     262          16 :       DO ider = -nder, nder
     263          14 :          dene = REAL(ider, KIND=dp)*delta
     264          14 :          CPASSERT(fhomo > ABS(dene))
     265          14 :          state%occupation(lhomo, nhomo) = fhomo + dene
     266          14 :          CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
     267         686 :          wfn(:, :, :, ider) = atom%orbitals%wfn
     268          16 :          state%occupation(lhomo, nhomo) = fhomo
     269             :       END DO
     270           2 :       IF (iw > 0) THEN
     271           2 :          WRITE (iw, '(A,T76,I5)') " Total number of electronic structure calculations ", 2*nder + 1
     272             :       END IF
     273             : 
     274           2 :       ovlp => atom%integrals%ovlp
     275             : 
     276           6 :       DO l = 0, state%maxl_occ
     277           4 :          IF (iw > 0) THEN
     278           4 :             WRITE (iw, '(A,T76,I5)') " Response derivatives for l quantum number ", l
     279             :          END IF
     280             :          ! occupied states
     281           8 :          DO i = 1, MAX(state%maxn_occ(l), 1)
     282          32 :             rbasis(:, i, l) = wfn(:, i, l, 0)
     283             :          END DO
     284             :          ! differentiation
     285          16 :          DO ider = 1, nder
     286          12 :             i = MAX(state%maxn_occ(l), 1)
     287           4 :             SELECT CASE (ider)
     288             :             CASE (1)
     289          28 :                rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
     290             :             CASE (2)
     291          28 :                rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
     292             :             CASE (3)
     293             :                rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
     294          28 :                                                + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
     295             :             CASE DEFAULT
     296          12 :                CPABORT("")
     297             :             END SELECT
     298             :          END DO
     299             : 
     300             :          ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
     301           4 :          n = state%maxn_occ(l) + nder
     302           4 :          m = atom%basis%nbas(l)
     303          20 :          DO i = 1, n
     304          40 :             DO j = 1, i - 1
     305        2304 :                o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
     306         184 :                rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
     307             :             END DO
     308        1536 :             o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
     309         116 :             rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
     310             :          END DO
     311             : 
     312             :          ! check
     313          16 :          ALLOCATE (amat(n, n))
     314        2320 :          amat(1:n, 1:n) = MATMUL(TRANSPOSE(rbasis(1:m, 1:n, l)), MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
     315          20 :          DO i = 1, n
     316          20 :             amat(i, i) = amat(i, i) - 1._dp
     317             :          END DO
     318          84 :          IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
     319           0 :             IF (iw > 0) WRITE (iw, '(A,G20.10)') " Orthogonality error  ", MAXVAL(ABS(amat))
     320             :          END IF
     321           4 :          DEALLOCATE (amat)
     322             : 
     323             :          ! Quickstep normalization
     324           4 :          expzet = 0.25_dp*REAL(2*l + 3, dp)
     325           4 :          prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
     326          30 :          DO i = 1, m
     327          24 :             zeta = (2._dp*atom%basis%am(i, l))**expzet
     328         124 :             qbasis(i, 1:n, l) = rbasis(i, 1:n, l)*prefac/zeta
     329             :          END DO
     330             : 
     331             :       END DO
     332             : 
     333             :       ! check for condition numbers
     334           2 :       IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Valence Response Basis Sets"
     335           2 :       CALL init_orbital_pointers(lmat)
     336           2 :       CALL init_spherical_harmonics(lmat, 0)
     337          10 :       DO ider = 0, nder
     338           8 :          NULLIFY (basis_vrb)
     339         152 :          ALLOCATE (basis_vrb)
     340             :          NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
     341             :                   basis_vrb%dbf, basis_vrb%ddbf)
     342             :          ! allocate and initialize the atomic grid
     343             :          NULLIFY (basis_vrb%grid)
     344           8 :          CALL allocate_grid_atom(basis_vrb%grid)
     345           8 :          CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
     346           8 :          basis_vrb%grid%nr = ngp
     347             :          !
     348           8 :          basis_vrb%eps_eig = basis_ref%eps_eig
     349           8 :          basis_vrb%geometrical = .FALSE.
     350           8 :          basis_vrb%basis_type = CGTO_BASIS
     351         104 :          basis_vrb%nprim = basis%nprim
     352          56 :          basis_vrb%nbas = 0
     353          24 :          DO l = 0, state%maxl_occ
     354          24 :             basis_vrb%nbas(l) = state%maxn_occ(l) + ider
     355             :          END DO
     356          56 :          m = MAXVAL(basis_vrb%nprim)
     357          56 :          n = MAXVAL(basis_vrb%nbas)
     358          24 :          ALLOCATE (basis_vrb%am(m, 0:lmat))
     359         680 :          basis_vrb%am = basis%am
     360             :          ! contractions
     361          40 :          ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
     362          24 :          DO l = 0, state%maxl_occ
     363          16 :             m = basis_vrb%nprim(l)
     364          16 :             n = basis_vrb%nbas(l)
     365         304 :             basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
     366             :          END DO
     367             : 
     368             :          ! initialize basis function on a radial grid
     369           8 :          nr = basis_vrb%grid%nr
     370          56 :          m = MAXVAL(basis_vrb%nbas)
     371          40 :          ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
     372          40 :          ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
     373          40 :          ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
     374       48176 :          basis_vrb%bf = 0._dp
     375       48176 :          basis_vrb%dbf = 0._dp
     376       48176 :          basis_vrb%ddbf = 0._dp
     377          56 :          DO l = 0, lmat
     378         152 :             DO i = 1, basis_vrb%nprim(l)
     379          96 :                al = basis_vrb%am(i, l)
     380       38544 :                DO k = 1, nr
     381       38400 :                   rk = basis_vrb%grid%rad(k)
     382       38400 :                   ear = EXP(-al*basis_vrb%grid%rad(k)**2)
     383      134496 :                   DO j = 1, basis_vrb%nbas(l)
     384       96000 :                      basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
     385             :                      basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
     386       96000 :                                               + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
     387             :                      basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
     388             :                                                (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
     389      134400 :                                                 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
     390             :                   END DO
     391             :                END DO
     392             :             END DO
     393             :          END DO
     394             : 
     395           8 :          IF (iw > 0) THEN
     396           8 :             CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
     397           8 :             WRITE (iw, '(A,A)') " Basis set     ", TRIM(abas)
     398             :          END IF
     399           8 :          crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
     400           8 :          cradx = crad*1.00_dp
     401           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     402           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     403           8 :          cradx = crad*1.10_dp
     404           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     405           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     406           8 :          cradx = crad*1.20_dp
     407           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     408           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     409          34 :          vbasis(ider)%basis => basis_vrb
     410             :       END DO
     411           2 :       CALL deallocate_orbital_pointers
     412           2 :       CALL deallocate_spherical_harmonics
     413             : 
     414             :       ! get density maximum
     415           6 :       ALLOCATE (rho(basis%grid%nr))
     416           2 :       CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
     417           2 :       CALL atom_density(rho(:), atom%orbitals%pmat, atom%basis, maxl, typ="RHO")
     418         802 :       n = SUM(MAXLOC(rho(:)))
     419           2 :       rmax = basis%grid%rad(n)
     420           2 :       IF (rmax < 0.1_dp) rmax = 1.0_dp
     421           2 :       DEALLOCATE (rho)
     422             : 
     423             :       ! generate polarization sets
     424           2 :       maxl = atom%state%maxl_occ
     425           2 :       CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto)
     426           2 :       num_pol = num_gto
     427           2 :       IF (num_gto > 0) THEN
     428           1 :          IF (iw > 0) THEN
     429           1 :             WRITE (iw, '(/,A)') " Polarization basis set  "
     430             :          END IF
     431           7 :          ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto))
     432         169 :          pbasis = 0.0_dp
     433             :          ! optimize exponents
     434           1 :          lval = maxl + 1
     435           1 :          zval = SQRT(REAL(2*lval + 2, dp))*REAL(lval + 1, dp)/(2._dp*rmax)
     436           1 :          aval = atom%basis%am(1, 0)
     437           1 :          cval = 2.5_dp
     438           1 :          rconf = atom%potential%scon
     439           1 :          CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section)
     440             :          ! calculate contractions
     441           5 :          DO i = 1, num_gto
     442           5 :             alp(i) = aval*cval**(i - 1)
     443             :          END DO
     444           3 :          ALLOCATE (rho(num_gto))
     445           5 :          DO l = maxl + 1, MIN(maxl + num_gto, 7)
     446           4 :             zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
     447           4 :             CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l))
     448           4 :             IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
     449           5 :                " Polarization basis set contraction for lval=", l, "zval=", zval
     450             :          END DO
     451           1 :          DEALLOCATE (rho)
     452             :       END IF
     453             : 
     454             :       ! generate valence expansion sets
     455           2 :       maxl = atom%state%maxl_occ
     456           2 :       CALL section_vals_val_get(grb_section, "NUM_GTO_EXTENDED", i_val=num_gto)
     457           2 :       CALL section_vals_val_get(grb_section, "EXTENSION_BASIS", i_vals=num_bas)
     458           2 :       next_bas(0:lmat) = 0
     459           2 :       IF (num_bas(1) == -1) THEN
     460           0 :          DO l = 0, maxl
     461           0 :             next_bas(l) = maxl - l + 1
     462             :          END DO
     463             :       ELSE
     464           2 :          n = MIN(SIZE(num_bas, 1), 4)
     465           6 :          next_bas(0:n - 1) = num_bas(1:n)
     466             :       END IF
     467           2 :       next_prim = 0
     468          14 :       DO l = 0, lmat
     469          14 :          IF (next_bas(l) > 0) next_prim(l) = num_gto
     470             :       END DO
     471           2 :       IF (iw > 0) THEN
     472           2 :          CALL basis_label(abas, next_prim, next_bas)
     473           2 :          WRITE (iw, '(/,A,A)') " Extension basis set     ", TRIM(abas)
     474             :       END IF
     475          14 :       n = MAXVAL(next_prim)
     476          14 :       m = MAXVAL(next_bas)
     477          11 :       ALLOCATE (ebasis(n, n, 0:lmat), ale(n))
     478           2 :       basis_vrb => vbasis(0)%basis
     479           2 :       amin = atom%basis%aval/atom%basis%cval**1.5_dp
     480           6 :       DO i = 1, n
     481           6 :          ale(i) = amin*atom%basis%cval**(i - 1)
     482             :       END DO
     483         134 :       ebasis = 0._dp
     484           5 :       ALLOCATE (rho(n))
     485           2 :       rconf = 2.0_dp*atom%potential%scon
     486          14 :       DO l = 0, lmat
     487          12 :          IF (next_bas(l) < 1) CYCLE
     488           2 :          zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
     489           2 :          CALL hydrogenic(zval, rconf, l, ale, n, rho, ebasis(:, :, l))
     490           2 :          IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
     491           4 :             " Extension basis set contraction for lval=", l, "zval=", zval
     492             :       END DO
     493           2 :       DEALLOCATE (rho)
     494             :       ! check for condition numbers
     495           2 :       IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Extended Basis Sets"
     496           2 :       CALL init_orbital_pointers(lmat)
     497           2 :       CALL init_spherical_harmonics(lmat, 0)
     498          10 :       DO ider = 0, nder
     499           8 :          NULLIFY (basis_vrb)
     500         152 :          ALLOCATE (basis_vrb)
     501             :          NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
     502             :                   basis_vrb%dbf, basis_vrb%ddbf)
     503             :          ! allocate and initialize the atomic grid
     504             :          NULLIFY (basis_vrb%grid)
     505           8 :          CALL allocate_grid_atom(basis_vrb%grid)
     506           8 :          CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
     507           8 :          basis_vrb%grid%nr = ngp
     508             :          !
     509           8 :          basis_vrb%eps_eig = basis_ref%eps_eig
     510           8 :          basis_vrb%geometrical = .FALSE.
     511           8 :          basis_vrb%basis_type = CGTO_BASIS
     512         104 :          basis_vrb%nprim = basis%nprim + next_prim
     513          56 :          basis_vrb%nbas = 0
     514          24 :          DO l = 0, state%maxl_occ
     515          24 :             basis_vrb%nbas(l) = state%maxn_occ(l) + ider + next_bas(l)
     516             :          END DO
     517          56 :          m = MAXVAL(basis_vrb%nprim)
     518          24 :          ALLOCATE (basis_vrb%am(m, 0:lmat))
     519             :          ! exponents
     520           8 :          m = SIZE(basis%am, 1)
     521         680 :          basis_vrb%am(1:m, :) = basis%am(1:m, :)
     522             :          n = SIZE(ale, 1)
     523          24 :          DO l = 0, state%maxl_occ
     524          56 :             basis_vrb%am(m + 1:m + n, l) = ale(1:n)
     525             :          END DO
     526             :          ! contractions
     527          56 :          m = MAXVAL(basis_vrb%nprim)
     528          56 :          n = MAXVAL(basis_vrb%nbas)
     529          40 :          ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
     530        1664 :          basis_vrb%cm = 0.0_dp
     531          24 :          DO l = 0, state%maxl_occ
     532          16 :             m = basis%nprim(l)
     533          16 :             n = state%maxn_occ(l) + ider
     534         296 :             basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
     535          84 :             basis_vrb%cm(m + 1:m + next_prim(l), n + 1:n + next_bas(l), l) = ebasis(1:next_prim(l), 1:next_bas(l), l)
     536             :          END DO
     537             : 
     538             :          ! initialize basis function on a radial grid
     539           8 :          nr = basis_vrb%grid%nr
     540          56 :          m = MAXVAL(basis_vrb%nbas)
     541          40 :          ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
     542          40 :          ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
     543          40 :          ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
     544       67424 :          basis_vrb%bf = 0._dp
     545       67424 :          basis_vrb%dbf = 0._dp
     546       67424 :          basis_vrb%ddbf = 0._dp
     547          56 :          DO l = 0, lmat
     548         184 :             DO i = 1, basis_vrb%nprim(l)
     549         128 :                al = basis_vrb%am(i, l)
     550       51376 :                DO k = 1, nr
     551       51200 :                   rk = basis_vrb%grid%rad(k)
     552       51200 :                   ear = EXP(-al*basis_vrb%grid%rad(k)**2)
     553      227328 :                   DO j = 1, basis_vrb%nbas(l)
     554      176000 :                      basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
     555             :                      basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
     556      176000 :                                               + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
     557             :                      basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
     558             :                                                (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
     559      227200 :                                                 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
     560             :                   END DO
     561             :                END DO
     562             :             END DO
     563             :          END DO
     564             : 
     565           8 :          IF (iw > 0) THEN
     566           8 :             CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
     567           8 :             WRITE (iw, '(A,A)') " Basis set     ", TRIM(abas)
     568             :          END IF
     569           8 :          crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
     570           8 :          cradx = crad*1.00_dp
     571           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     572           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     573           8 :          cradx = crad*1.10_dp
     574           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     575           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     576           8 :          cradx = crad*1.20_dp
     577           8 :          CALL atom_basis_condnum(basis_vrb, cradx, cnum)
     578           8 :          IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
     579          34 :          vbasis(nder + 1 + ider)%basis => basis_vrb
     580             :       END DO
     581           2 :       CALL deallocate_orbital_pointers
     582           2 :       CALL deallocate_spherical_harmonics
     583             : 
     584             :       ! Tests for energy
     585           2 :       energy_ref = atom_ref%energy%etot
     586           2 :       IF (iw > 0) WRITE (iw, '(/,A,A)') " Basis set tests    "
     587           2 :       IF (iw > 0) WRITE (iw, '(T10,A,T59,F22.9)') " Reference Energy [a.u.]  ", energy_ref
     588          18 :       DO ider = 0, 2*nder + 1
     589             :          ! generate an atom type
     590          16 :          NULLIFY (atom_test)
     591          16 :          CALL create_atom_type(atom_test)
     592             :          CALL copy_atom_basics(atom_ref, atom_test, state=.TRUE., potential=.TRUE., &
     593          16 :                                optimization=.TRUE., xc=.TRUE.)
     594          16 :          basis_grb => vbasis(ider)%basis
     595          16 :          NULLIFY (orbitals)
     596         112 :          mo = MAXVAL(atom_test%state%maxn_calc)
     597         112 :          mb = MAXVAL(basis_grb%nbas)
     598          16 :          CALL create_atom_orbs(orbitals, mb, mo)
     599          16 :          CALL set_atom(atom_test, orbitals=orbitals, basis=basis_grb)
     600             :          ! calculate integrals
     601        3408 :          ALLOCATE (atint)
     602          16 :          CALL atom_int_setup(atint, basis_grb, potential=atom_test%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
     603          16 :          CALL atom_ppint_setup(atint, basis_grb, potential=atom_test%potential)
     604          16 :          IF (atom_test%pp_calc) THEN
     605          16 :             NULLIFY (atint%tzora, atint%hdkh)
     606             :          ELSE
     607             :             ! relativistic correction terms
     608           0 :             CALL atom_relint_setup(atint, basis_grb, atom_test%relativistic, zcore=REAL(atom_test%z, dp))
     609             :          END IF
     610          16 :          CALL set_atom(atom_test, integrals=atint)
     611             :          !
     612          16 :          CALL calculate_atom(atom_test, iw=0)
     613          16 :          IF (ider <= nder) THEN
     614           8 :             energy_vb(ider) = atom_test%energy%etot
     615          16 :             IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (VB)", ider, " Energy [a.u.]  ", &
     616          16 :                energy_ref - energy_vb(ider), energy_vb(ider)
     617             :          ELSE
     618           8 :             i = ider - nder - 1
     619           8 :             energy_ex(i) = atom_test%energy%etot
     620          16 :             IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (EX)", i, " Energy [a.u.]  ", &
     621          16 :                energy_ref - energy_ex(i), energy_ex(i)
     622             :          END IF
     623          16 :          CALL atom_int_release(atint)
     624          16 :          CALL atom_ppint_release(atint)
     625          16 :          CALL atom_relint_release(atint)
     626          16 :          DEALLOCATE (atom_test%state, atom_test%potential, atint)
     627          18 :          CALL release_atom_type(atom_test)
     628             :       END DO
     629             : 
     630             :       ! Quickstep normalization polarization basis
     631          18 :       DO l = 0, 7
     632          16 :          expzet = 0.25_dp*REAL(2*l + 3, dp)
     633          16 :          prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
     634          50 :          DO i = 1, num_pol
     635          32 :             zeta = (2._dp*alp(i))**expzet
     636         176 :             pbasis(i, 1:num_pol, l) = pbasis(i, 1:num_pol, l)*prefac/zeta
     637             :          END DO
     638             :       END DO
     639             :       ! Quickstep normalization extended basis
     640          14 :       DO l = 0, lmat
     641          12 :          expzet = 0.25_dp*REAL(2*l + 3, dp)
     642          12 :          prefac = SQRT(SQRT(pi)/2._dp**(l + 2)*dfac(2*l + 1))
     643          22 :          DO i = 1, next_prim(l)
     644           8 :             zeta = (2._dp*ale(i))**expzet
     645          32 :             ebasis(i, 1:next_bas(l), l) = ebasis(i, 1:next_bas(l), l)*prefac/zeta
     646             :          END DO
     647             :       END DO
     648             : 
     649             :       ! Print basis sets
     650           2 :       CALL section_vals_val_get(grb_section, "NAME_BODY", c_val=basname)
     651           2 :       CALL open_file(file_name="GRB_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iunit)
     652             :       ! header info
     653           8 :       headline = ""
     654           2 :       headline(1) = "#"
     655           2 :       headline(2) = "# Generated with CP2K Atom Code"
     656           2 :       headline(3) = "#"
     657           2 :       CALL grb_print_basis(header=headline, iunit=iunit)
     658             :       ! valence basis
     659           2 :       basline(1) = ""
     660           2 :       WRITE (basline(1), "(T2,A)") ADJUSTL(ptable(atom_ref%z)%symbol)
     661          10 :       DO ider = 0, nder
     662           8 :          basline(1) = ""
     663           8 :          WRITE (basline(1), "(T2,A,T5,A,I1)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-VAL-", ider
     664             :          CALL grb_print_basis(header=basline, nprim=vbasis(ider)%basis%nprim(0), nbas=vbasis(ider)%basis%nbas, &
     665          10 :                               al=vbasis(ider)%basis%am(:, 0), gcc=qbasis, iunit=iunit)
     666             :       END DO
     667             :       ! polarization basis
     668           2 :       maxl = atom_ref%state%maxl_occ
     669           6 :       DO l = maxl + 1, MIN(maxl + num_pol, 7)
     670           4 :          nbas = 0
     671          14 :          DO i = maxl + 1, l
     672          14 :             nbas(i) = l - i + 1
     673             :          END DO
     674           4 :          i = l - maxl
     675           4 :          basline(1) = ""
     676           4 :          WRITE (basline(1), "(T2,A,T5,A,I1)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-POL-", i
     677           6 :          CALL grb_print_basis(header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit)
     678             :       END DO
     679             :       ! extension set
     680          14 :       IF (SUM(next_bas) > 0) THEN
     681           1 :          basline(1) = ""
     682           1 :          WRITE (basline(1), "(T2,A,T5,A)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-EXT"
     683           1 :          CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit)
     684             :       END IF
     685             :       !
     686           2 :       CALL close_file(unit_number=iunit)
     687             : 
     688             :       ! clean up
     689           2 :       IF (ALLOCATED(pbasis)) DEALLOCATE (pbasis)
     690           2 :       IF (ALLOCATED(alp)) DEALLOCATE (alp)
     691           2 :       IF (ALLOCATED(ebasis)) DEALLOCATE (ebasis)
     692           2 :       DEALLOCATE (wfn, rbasis, qbasis, ale)
     693             : 
     694          24 :       DO ider = 0, 10
     695          24 :          IF (ASSOCIATED(vbasis(ider)%basis)) THEN
     696          16 :             CALL release_atom_basis(vbasis(ider)%basis)
     697          16 :             DEALLOCATE (vbasis(ider)%basis)
     698             :          END IF
     699             :       END DO
     700             : 
     701           2 :       CALL atom_int_release(atom%integrals)
     702           2 :       CALL atom_ppint_release(atom%integrals)
     703           2 :       CALL atom_relint_release(atom%integrals)
     704           2 :       CALL release_atom_basis(basis)
     705           2 :       DEALLOCATE (atom%potential, atom%state, atom%integrals, basis)
     706           2 :       CALL release_atom_type(atom)
     707             : 
     708           2 :       IF (iw > 0) WRITE (iw, '(" ",79("*"))')
     709             : 
     710          18 :    END SUBROUTINE atom_grb_construction
     711             : 
     712             : ! **************************************************************************************************
     713             : !> \brief Print geometrical response basis set.
     714             : !> \param header  banner to print on top of the basis set
     715             : !> \param nprim   number of primitive exponents
     716             : !> \param nbas    number of basis functions for the given angular momentum
     717             : !> \param al      list of the primitive exponents
     718             : !> \param gcc     array of contraction coefficients of size
     719             : !>                (index-of-the-primitive-exponent, index-of-the-contraction-set, angular-momentum)
     720             : !> \param iunit   output file unit
     721             : !> \par History
     722             : !>    * 11.2016 created [Juerg Hutter]
     723             : ! **************************************************************************************************
     724          15 :    SUBROUTINE grb_print_basis(header, nprim, nbas, al, gcc, iunit)
     725             :       CHARACTER(len=*), DIMENSION(:), INTENT(IN), &
     726             :          OPTIONAL                                        :: header
     727             :       INTEGER, INTENT(IN), OPTIONAL                      :: nprim
     728             :       INTEGER, DIMENSION(0:), INTENT(IN), OPTIONAL       :: nbas
     729             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: al
     730             :       REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN), &
     731             :          OPTIONAL                                        :: gcc
     732             :       INTEGER, INTENT(IN)                                :: iunit
     733             : 
     734             :       INTEGER                                            :: i, j, l, lmax, lmin, nval
     735             : 
     736          15 :       IF (PRESENT(header)) THEN
     737          34 :          DO i = 1, SIZE(header, 1)
     738          34 :             IF (header(i) .NE. "") THEN
     739          19 :                WRITE (iunit, "(A)") TRIM(header(i))
     740             :             END IF
     741             :          END DO
     742             :       END IF
     743             : 
     744          15 :       IF (PRESENT(nprim)) THEN
     745          13 :          IF (nprim > 0) THEN
     746          13 :             CPASSERT(PRESENT(nbas))
     747          13 :             CPASSERT(PRESENT(al))
     748          13 :             CPASSERT(PRESENT(gcc))
     749             : 
     750          47 :             DO i = LBOUND(nbas, 1), UBOUND(nbas, 1)
     751          21 :                IF (nbas(i) > 0) THEN
     752          13 :                   lmin = i
     753          13 :                   EXIT
     754             :                END IF
     755             :             END DO
     756          63 :             DO i = UBOUND(nbas, 1), LBOUND(nbas, 1), -1
     757          63 :                IF (nbas(i) > 0) THEN
     758          13 :                   lmax = i
     759          13 :                   EXIT
     760             :                END IF
     761             :             END DO
     762             : 
     763          13 :             nval = lmax
     764          13 :             WRITE (iunit, *) "  1"
     765          13 :             WRITE (iunit, "(40I3)") nval, lmin, lmax, nprim, (nbas(l), l=lmin, lmax)
     766          81 :             DO i = nprim, 1, -1
     767          68 :                WRITE (iunit, "(G20.12)", advance="no") al(i)
     768         212 :                DO l = lmin, lmax
     769         544 :                   DO j = 1, nbas(l)
     770         476 :                      WRITE (iunit, "(F16.10)", advance="no") gcc(i, j, l)
     771             :                   END DO
     772             :                END DO
     773          81 :                WRITE (iunit, *)
     774             :             END DO
     775          13 :             WRITE (iunit, *)
     776             :          END IF
     777             :       END IF
     778             : 
     779          15 :    END SUBROUTINE grb_print_basis
     780             : 
     781             : ! **************************************************************************************************
     782             : !> \brief Compose the basis set label:
     783             : !>        (np(0)'s'np(1)'p'...) -> [nb(0)'s'nb(1)'p'...] .
     784             : !> \param label  basis set label
     785             : !> \param np     number of primitive basis functions per angular momentum
     786             : !> \param nb     number of contracted basis functions per angular momentum
     787             : !> \par History
     788             : !>    * 11.2016 created [Juerg Hutter]
     789             : ! **************************************************************************************************
     790          18 :    SUBROUTINE basis_label(label, np, nb)
     791             :       CHARACTER(len=*), INTENT(out)                      :: label
     792             :       INTEGER, DIMENSION(0:), INTENT(in)                 :: np, nb
     793             : 
     794             :       INTEGER                                            :: i, l, lmax
     795             :       CHARACTER(len=1), DIMENSION(0:7), PARAMETER :: &
     796             :          lq = (/"s", "p", "d", "f", "g", "h", "i", "k"/)
     797             : 
     798          18 :       label = ""
     799          36 :       lmax = MIN(UBOUND(np, 1), UBOUND(nb, 1), 7)
     800          18 :       i = 1
     801          18 :       label(i:i) = "("
     802         126 :       DO l = 0, lmax
     803         126 :          IF (np(l) > 0) THEN
     804          34 :             i = i + 1
     805          34 :             IF (np(l) > 9) THEN
     806           8 :                WRITE (label(i:i + 1), "(I2)") np(l)
     807           8 :                i = i + 2
     808             :             ELSE
     809          26 :                WRITE (label(i:i), "(I1)") np(l)
     810          26 :                i = i + 1
     811             :             END IF
     812          34 :             label(i:i) = lq(l)
     813             :          END IF
     814             :       END DO
     815          18 :       i = i + 1
     816          18 :       label(i:i + 6) = ") --> ["
     817          18 :       i = i + 6
     818         126 :       DO l = 0, lmax
     819         126 :          IF (nb(l) > 0) THEN
     820          34 :             i = i + 1
     821          34 :             IF (nb(l) > 9) THEN
     822           0 :                WRITE (label(i:i + 1), "(I2)") nb(l)
     823           0 :                i = i + 2
     824             :             ELSE
     825          34 :                WRITE (label(i:i), "(I1)") nb(l)
     826          34 :                i = i + 1
     827             :             END IF
     828          34 :             label(i:i) = lq(l)
     829             :          END IF
     830             :       END DO
     831          18 :       i = i + 1
     832          18 :       label(i:i) = "]"
     833             : 
     834          18 :    END SUBROUTINE basis_label
     835             : 
     836             : ! **************************************************************************************************
     837             : !> \brief Compute the total energy for the given atomic kind and basis set.
     838             : !> \param atom    information about the atomic kind
     839             : !> \param basis   basis set to fit
     840             : !> \param afun    (output) atomic total energy
     841             : !> \param iw      output file unit
     842             : !> \par History
     843             : !>    * 11.2016 created [Juerg Hutter]
     844             : ! **************************************************************************************************
     845         170 :    SUBROUTINE grb_fit(atom, basis, afun, iw)
     846             :       TYPE(atom_type), POINTER                           :: atom
     847             :       TYPE(atom_basis_type), POINTER                     :: basis
     848             :       REAL(dp), INTENT(OUT)                              :: afun
     849             :       INTEGER, INTENT(IN)                                :: iw
     850             : 
     851             :       INTEGER                                            :: do_eric, do_erie, reltyp, zval
     852             :       LOGICAL                                            :: eri_c, eri_e
     853             :       TYPE(atom_integrals), POINTER                      :: atint
     854             :       TYPE(atom_potential_type), POINTER                 :: pot
     855             : 
     856       36890 :       ALLOCATE (atint)
     857             :       ! calculate integrals
     858         170 :       NULLIFY (pot)
     859         170 :       eri_c = .FALSE.
     860         170 :       eri_e = .FALSE.
     861         170 :       pot => atom%potential
     862         170 :       zval = atom%z
     863         170 :       reltyp = atom%relativistic
     864         170 :       do_eric = atom%coulomb_integral_type
     865         170 :       do_erie = atom%exchange_integral_type
     866         170 :       IF (do_eric == do_analytic) eri_c = .TRUE.
     867         170 :       IF (do_erie == do_analytic) eri_e = .TRUE.
     868             :       ! general integrals
     869         170 :       CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
     870             :       ! potential
     871         170 :       CALL atom_ppint_setup(atint, basis, potential=pot)
     872         170 :       IF (atom%pp_calc) THEN
     873         170 :          NULLIFY (atint%tzora, atint%hdkh)
     874             :       ELSE
     875             :          ! relativistic correction terms
     876           0 :          CALL atom_relint_setup(atint, basis, reltyp, zcore=REAL(zval, dp))
     877             :       END IF
     878         170 :       CALL set_atom(atom, basis=basis)
     879         170 :       CALL set_atom(atom, integrals=atint)
     880         170 :       CALL calculate_atom(atom, iw)
     881         170 :       afun = atom%energy%etot
     882         170 :       CALL atom_int_release(atint)
     883         170 :       CALL atom_ppint_release(atint)
     884         170 :       CALL atom_relint_release(atint)
     885         170 :       DEALLOCATE (atint)
     886         170 :    END SUBROUTINE grb_fit
     887             : 
     888             : ! **************************************************************************************************
     889             : !> \brief Copy basic information about the atomic kind.
     890             : !> \param atom_ref      atom to copy
     891             : !> \param atom          new atom to create
     892             : !> \param state         also copy electronic state and occupation numbers
     893             : !> \param potential     also copy pseudo-potential
     894             : !> \param optimization  also copy optimization procedure
     895             : !> \param xc            also copy the XC input section
     896             : !> \par History
     897             : !>    * 11.2016 created [Juerg Hutter]
     898             : ! **************************************************************************************************
     899          18 :    SUBROUTINE copy_atom_basics(atom_ref, atom, state, potential, optimization, xc)
     900             :       TYPE(atom_type), POINTER                           :: atom_ref, atom
     901             :       LOGICAL, INTENT(IN), OPTIONAL                      :: state, potential, optimization, xc
     902             : 
     903          18 :       atom%z = atom_ref%z
     904          18 :       atom%zcore = atom_ref%zcore
     905          18 :       atom%pp_calc = atom_ref%pp_calc
     906          18 :       atom%method_type = atom_ref%method_type
     907          18 :       atom%relativistic = atom_ref%relativistic
     908          18 :       atom%coulomb_integral_type = atom_ref%coulomb_integral_type
     909          18 :       atom%exchange_integral_type = atom_ref%exchange_integral_type
     910             : 
     911          18 :       NULLIFY (atom%potential, atom%state, atom%xc_section)
     912          18 :       NULLIFY (atom%basis, atom%integrals, atom%orbitals, atom%fmat)
     913             : 
     914          18 :       IF (PRESENT(state)) THEN
     915          18 :          IF (state) THEN
     916        6660 :             ALLOCATE (atom%state)
     917          18 :             atom%state = atom_ref%state
     918             :          END IF
     919             :       END IF
     920             : 
     921          18 :       IF (PRESENT(potential)) THEN
     922          18 :          IF (potential) THEN
     923       94914 :             ALLOCATE (atom%potential)
     924          18 :             atom%potential = atom_ref%potential
     925             :          END IF
     926             :       END IF
     927             : 
     928          18 :       IF (PRESENT(optimization)) THEN
     929          18 :          IF (optimization) THEN
     930          18 :             atom%optimization = atom_ref%optimization
     931             :          END IF
     932             :       END IF
     933             : 
     934          18 :       IF (PRESENT(xc)) THEN
     935          18 :          IF (xc) THEN
     936          18 :             atom%xc_section => atom_ref%xc_section
     937             :          END IF
     938             :       END IF
     939             : 
     940          18 :    END SUBROUTINE copy_atom_basics
     941             : 
     942             : ! **************************************************************************************************
     943             : !> \brief Optimise a geometrical response basis set.
     944             : !> \param atom            information about the atomic kind
     945             : !> \param basis           basis set to fit
     946             : !> \param iunit           output file unit
     947             : !> \param powell_section  POWELL input section
     948             : !> \par History
     949             : !>    * 11.2016 created [Juerg Hutter]
     950             : ! **************************************************************************************************
     951           2 :    SUBROUTINE atom_fit_grb(atom, basis, iunit, powell_section)
     952             :       TYPE(atom_type), POINTER                           :: atom
     953             :       TYPE(atom_basis_type), POINTER                     :: basis
     954             :       INTEGER, INTENT(IN)                                :: iunit
     955             :       TYPE(section_vals_type), POINTER                   :: powell_section
     956             : 
     957             :       INTEGER                                            :: i, k, l, ll, n10, nr
     958             :       REAL(KIND=dp)                                      :: al, cnum, crad, cradx, ear, fopt, rk
     959           2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: x
     960             :       TYPE(opt_state_type)                               :: ostate
     961             : 
     962           0 :       CPASSERT(basis%geometrical)
     963             : 
     964           2 :       CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
     965           2 :       CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
     966           2 :       CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
     967             : 
     968           2 :       ostate%nvar = 2
     969           2 :       ALLOCATE (x(2))
     970           2 :       x(1) = SQRT(basis%aval)
     971           2 :       x(2) = SQRT(basis%cval)
     972             : 
     973           2 :       ostate%nf = 0
     974           2 :       ostate%iprint = 1
     975           2 :       ostate%unit = iunit
     976             : 
     977           2 :       ostate%state = 0
     978           2 :       IF (iunit > 0) THEN
     979           2 :          WRITE (iunit, '(/," POWELL| Start optimization procedure")')
     980           2 :          WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
     981             :       END IF
     982           2 :       n10 = MAX(ostate%maxfun/100, 1)
     983             : 
     984           2 :       fopt = HUGE(0._dp)
     985             : 
     986             :       DO
     987             : 
     988         176 :          IF (ostate%state == 2) THEN
     989        7310 :             basis%am = 0._dp
     990        1190 :             DO l = 0, lmat
     991        3230 :                DO i = 1, basis%nbas(l)
     992        2040 :                   ll = i - 1 + basis%start(l)
     993        3060 :                   basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
     994             :                END DO
     995             :             END DO
     996         170 :             basis%aval = x(1)*x(1)
     997         170 :             basis%cval = x(2)*x(2)
     998     2455310 :             basis%bf = 0._dp
     999     2455310 :             basis%dbf = 0._dp
    1000     2455310 :             basis%ddbf = 0._dp
    1001         170 :             nr = basis%grid%nr
    1002        1190 :             DO l = 0, lmat
    1003        3230 :                DO i = 1, basis%nbas(l)
    1004        2040 :                   al = basis%am(i, l)
    1005      819060 :                   DO k = 1, nr
    1006      816000 :                      rk = basis%grid%rad(k)
    1007      816000 :                      ear = EXP(-al*basis%grid%rad(k)**2)
    1008      816000 :                      basis%bf(k, i, l) = rk**l*ear
    1009      816000 :                      basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
    1010             :                      basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
    1011      818040 :                                             2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
    1012             :                   END DO
    1013             :                END DO
    1014             :             END DO
    1015         170 :             CALL grb_fit(atom, basis, ostate%f, 0)
    1016         170 :             fopt = MIN(fopt, ostate%f)
    1017             :          END IF
    1018             : 
    1019         176 :          IF (ostate%state == -1) EXIT
    1020             : 
    1021         174 :          CALL powell_optimize(ostate%nvar, x, ostate)
    1022             : 
    1023         174 :          IF (ostate%nf == 2 .AND. iunit > 0) THEN
    1024           2 :             WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
    1025             :          END IF
    1026         176 :          IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
    1027             :             WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
    1028           2 :                INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
    1029             :          END IF
    1030             : 
    1031             :       END DO
    1032             : 
    1033           2 :       ostate%state = 8
    1034           2 :       CALL powell_optimize(ostate%nvar, x, ostate)
    1035             : 
    1036           2 :       IF (iunit > 0) THEN
    1037           2 :          WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
    1038           2 :          WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
    1039             :       END IF
    1040             :       ! x->basis
    1041          86 :       basis%am = 0._dp
    1042          14 :       DO l = 0, lmat
    1043          38 :          DO i = 1, basis%nbas(l)
    1044          24 :             ll = i - 1 + basis%start(l)
    1045          36 :             basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
    1046             :          END DO
    1047             :       END DO
    1048           2 :       basis%aval = x(1)*x(1)
    1049           2 :       basis%cval = x(2)*x(2)
    1050       28886 :       basis%bf = 0._dp
    1051       28886 :       basis%dbf = 0._dp
    1052       28886 :       basis%ddbf = 0._dp
    1053           2 :       nr = basis%grid%nr
    1054          14 :       DO l = 0, lmat
    1055          38 :          DO i = 1, basis%nbas(l)
    1056          24 :             al = basis%am(i, l)
    1057        9636 :             DO k = 1, nr
    1058        9600 :                rk = basis%grid%rad(k)
    1059        9600 :                ear = EXP(-al*basis%grid%rad(k)**2)
    1060        9600 :                basis%bf(k, i, l) = rk**l*ear
    1061        9600 :                basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
    1062             :                basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
    1063        9624 :                                       2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
    1064             :             END DO
    1065             :          END DO
    1066             :       END DO
    1067             : 
    1068           2 :       DEALLOCATE (x)
    1069             : 
    1070             :       ! final result
    1071           2 :       IF (iunit > 0) THEN
    1072           2 :          WRITE (iunit, '(/,A)') " Optimized Geometrical GTO basis set"
    1073           2 :          WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", basis%aval, &
    1074           4 :             " Proportionality factor: ", basis%cval
    1075          14 :          DO l = 0, lmat
    1076          14 :             WRITE (iunit, '(T41,A,I2,T76,I5)') " Number of exponents for l=", l, basis%nbas(l)
    1077             :          END DO
    1078             :       END IF
    1079             : 
    1080           2 :       IF (iunit > 0) WRITE (iunit, '(/,A)') " Condition number of uncontracted basis set"
    1081           2 :       crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
    1082           2 :       CALL init_orbital_pointers(lmat)
    1083           2 :       CALL init_spherical_harmonics(lmat, 0)
    1084           2 :       cradx = crad*1.00_dp
    1085           2 :       CALL atom_basis_condnum(basis, cradx, cnum)
    1086           2 :       IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
    1087           2 :       cradx = crad*1.10_dp
    1088           2 :       CALL atom_basis_condnum(basis, cradx, cnum)
    1089           2 :       IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
    1090           2 :       cradx = crad*1.20_dp
    1091           2 :       CALL atom_basis_condnum(basis, cradx, cnum)
    1092           2 :       IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
    1093           2 :       CALL deallocate_orbital_pointers
    1094           2 :       CALL deallocate_spherical_harmonics
    1095             : 
    1096           8 :    END SUBROUTINE atom_fit_grb
    1097             : 
    1098             : ! **************************************************************************************************
    1099             : !> \brief Optimize 'aval' and 'cval' parameters which define the geometrical response basis set.
    1100             : !> \param zval            nuclear charge
    1101             : !> \param rconf           confinement radius
    1102             : !> \param lval            angular momentum
    1103             : !> \param aval            (input/output) exponent of the first Gaussian basis function in the series
    1104             : !> \param cval            (input/output) factor of geometrical series
    1105             : !> \param nbas            number of basis functions
    1106             : !> \param iunit           output file unit
    1107             : !> \param powell_section  POWELL input section
    1108             : !> \par History
    1109             : !>    * 11.2016 created [Juerg Hutter]
    1110             : ! **************************************************************************************************
    1111           1 :    SUBROUTINE atom_fit_pol(zval, rconf, lval, aval, cval, nbas, iunit, powell_section)
    1112             :       REAL(KIND=dp), INTENT(IN)                          :: zval, rconf
    1113             :       INTEGER, INTENT(IN)                                :: lval
    1114             :       REAL(KIND=dp), INTENT(INOUT)                       :: aval, cval
    1115             :       INTEGER, INTENT(IN)                                :: nbas, iunit
    1116             :       TYPE(section_vals_type), POINTER                   :: powell_section
    1117             : 
    1118             :       INTEGER                                            :: i, n10
    1119             :       REAL(KIND=dp)                                      :: fopt, x(2)
    1120             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: am, ener
    1121             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: orb
    1122             :       TYPE(opt_state_type)                               :: ostate
    1123             : 
    1124           8 :       ALLOCATE (am(nbas), ener(nbas), orb(nbas, nbas))
    1125             : 
    1126           1 :       CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
    1127           1 :       CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
    1128           1 :       CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
    1129             : 
    1130           1 :       ostate%nvar = 2
    1131           1 :       x(1) = SQRT(aval)
    1132           1 :       x(2) = SQRT(cval)
    1133             : 
    1134           1 :       ostate%nf = 0
    1135           1 :       ostate%iprint = 1
    1136           1 :       ostate%unit = iunit
    1137             : 
    1138           1 :       ostate%state = 0
    1139           1 :       IF (iunit > 0) THEN
    1140           1 :          WRITE (iunit, '(/," POWELL| Start optimization procedure")')
    1141           1 :          WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
    1142             :       END IF
    1143           1 :       n10 = MAX(ostate%maxfun/100, 1)
    1144             : 
    1145           1 :       fopt = HUGE(0._dp)
    1146             : 
    1147             :       DO
    1148             : 
    1149          72 :          IF (ostate%state == 2) THEN
    1150          69 :             aval = x(1)*x(1)
    1151          69 :             cval = x(2)*x(2)
    1152         345 :             DO i = 1, nbas
    1153         345 :                am(i) = aval*cval**(i - 1)
    1154             :             END DO
    1155          69 :             CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
    1156          69 :             ostate%f = ener(1)
    1157          69 :             fopt = MIN(fopt, ostate%f)
    1158             :          END IF
    1159             : 
    1160          72 :          IF (ostate%state == -1) EXIT
    1161             : 
    1162          71 :          CALL powell_optimize(ostate%nvar, x, ostate)
    1163             : 
    1164          71 :          IF (ostate%nf == 2 .AND. iunit > 0) THEN
    1165           1 :             WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
    1166             :          END IF
    1167          72 :          IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
    1168             :             WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
    1169           1 :                INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
    1170             :          END IF
    1171             : 
    1172             :       END DO
    1173             : 
    1174           1 :       ostate%state = 8
    1175           1 :       CALL powell_optimize(ostate%nvar, x, ostate)
    1176             : 
    1177           1 :       IF (iunit > 0) THEN
    1178           1 :          WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
    1179           1 :          WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
    1180             :       END IF
    1181             :       ! x->basis
    1182           1 :       aval = x(1)*x(1)
    1183           1 :       cval = x(2)*x(2)
    1184             : 
    1185             :       ! final result
    1186           1 :       IF (iunit > 0) THEN
    1187           1 :          WRITE (iunit, '(/,A,T51,A,T76,I5)') " Optimized Polarization basis set", &
    1188           2 :             " Number of exponents:", nbas
    1189           1 :          WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", aval, &
    1190           2 :             " Proportionality factor: ", cval
    1191             :       END IF
    1192             : 
    1193           1 :       DEALLOCATE (am, ener, orb)
    1194             : 
    1195           1 :    END SUBROUTINE atom_fit_pol
    1196             : 
    1197             : ! **************************************************************************************************
    1198             : !> \brief Calculate orbitals of a hydrogen-like atom.
    1199             : !> \param zval   nuclear charge
    1200             : !> \param rconf  confinement radius
    1201             : !> \param lval   angular momentum
    1202             : !> \param am     list of basis functions' exponents
    1203             : !> \param nbas   number of basis functions
    1204             : !> \param ener   orbital energies
    1205             : !> \param orb    expansion coefficients of atomic wavefunctions
    1206             : !> \par History
    1207             : !>    * 11.2016 created [Juerg Hutter]
    1208             : ! **************************************************************************************************
    1209          75 :    SUBROUTINE hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
    1210             :       REAL(KIND=dp), INTENT(IN)                          :: zval, rconf
    1211             :       INTEGER, INTENT(IN)                                :: lval
    1212             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: am
    1213             :       INTEGER, INTENT(IN)                                :: nbas
    1214             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: ener
    1215             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: orb
    1216             : 
    1217             :       INTEGER                                            :: info, k, lwork, n
    1218             :       REAL(KIND=dp)                                      :: cf
    1219          75 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: w, work
    1220          75 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: confmat, hmat, potmat, smat, tmat
    1221             : 
    1222          75 :       n = nbas
    1223        1200 :       ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
    1224             :       ! calclulate overlap matrix
    1225          75 :       CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
    1226             :       ! calclulate kinetic energy matrix
    1227          75 :       CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
    1228             :       ! calclulate core potential matrix
    1229          75 :       CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
    1230             :       ! calclulate confinement potential matrix
    1231          75 :       cf = 0.1_dp
    1232          75 :       k = 10
    1233          75 :       CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
    1234             :       ! Hamiltionian
    1235        1575 :       hmat(1:n, 1:n) = tmat(1:n, 1:n) - zval*potmat(1:n, 1:n) + cf*confmat(1:n, 1:n)
    1236             :       ! solve
    1237          75 :       lwork = 100*n
    1238         375 :       ALLOCATE (w(n), work(lwork))
    1239          75 :       CALL lapack_ssygv(1, "V", "U", n, hmat, n, smat, n, w, work, lwork, info)
    1240          75 :       CPASSERT(info == 0)
    1241        1575 :       orb(1:n, 1:n) = hmat(1:n, 1:n)
    1242         375 :       ener(1:n) = w(1:n)
    1243          75 :       DEALLOCATE (w, work)
    1244          75 :       DEALLOCATE (smat, tmat, potmat, confmat, hmat)
    1245             : 
    1246          75 :    END SUBROUTINE hydrogenic
    1247             : 
    1248          44 : END MODULE atom_grb

Generated by: LCOV version 1.15