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