1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief routines that fit parameters for /from atomic calculations 8! ************************************************************************************************** 9MODULE atom_fit 10 USE atom_electronic_structure, ONLY: calculate_atom 11 USE atom_operators, ONLY: atom_int_release,& 12 atom_int_setup,& 13 atom_ppint_release,& 14 atom_ppint_setup,& 15 atom_relint_release,& 16 atom_relint_setup 17 USE atom_output, ONLY: atom_print_basis,& 18 atom_print_basis_file,& 19 atom_write_pseudo_param 20 USE atom_types, ONLY: & 21 GTO_BASIS, STO_BASIS, atom_basis_type, atom_gthpot_type, atom_integrals, atom_p_type, & 22 atom_potential_type, atom_type, create_opgrid, create_opmat, lmat, opgrid_type, & 23 opmat_type, release_opgrid, release_opmat, set_atom 24 USE atom_utils, ONLY: & 25 atom_completeness, atom_condnumber, atom_consistent_method, atom_core_density, & 26 atom_denmat, atom_density, atom_orbital_charge, atom_orbital_max, atom_orbital_nodes, & 27 atom_wfnr0, get_maxn_occ, integrate_grid 28 USE cp_files, ONLY: close_file,& 29 open_file 30 USE input_constants, ONLY: do_analytic,& 31 do_rhf_atom,& 32 do_rks_atom,& 33 do_rohf_atom,& 34 do_uhf_atom,& 35 do_uks_atom 36 USE input_section_types, ONLY: section_vals_type,& 37 section_vals_val_get 38 USE kinds, ONLY: dp 39 USE lapack, ONLY: lapack_sgesv 40 USE mathconstants, ONLY: fac,& 41 fourpi,& 42 pi 43 USE periodic_table, ONLY: ptable 44 USE physcon, ONLY: bohr,& 45 evolt 46 USE powell, ONLY: opt_state_type,& 47 powell_optimize 48#include "./base/base_uses.f90" 49 50 IMPLICIT NONE 51 52 PRIVATE 53 54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_fit' 55 56 PUBLIC :: atom_fit_density, atom_fit_basis, atom_fit_pseudo, atom_fit_kgpot 57 58 TYPE wfn_init 59 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn 60 END TYPE wfn_init 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief Fit the atomic electron density using a geometrical Gaussian basis set. 66!> \param atom information about the atomic kind 67!> \param num_gto number of Gaussian basis functions 68!> \param norder ... 69!> \param iunit output file unit 70!> \param powell_section POWELL input section 71!> \param results (output) array of num_gto+2 real numbers in the following format: 72!> starting exponent, factor of geometrical series, expansion coefficients(1:num_gto) 73! ************************************************************************************************** 74 SUBROUTINE atom_fit_density(atom, num_gto, norder, iunit, powell_section, results) 75 TYPE(atom_type), POINTER :: atom 76 INTEGER, INTENT(IN) :: num_gto, norder, iunit 77 TYPE(section_vals_type), OPTIONAL, POINTER :: powell_section 78 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: results 79 80 CHARACTER(len=*), PARAMETER :: routineN = 'atom_fit_density', & 81 routineP = moduleN//':'//routineN 82 83 INTEGER :: n10 84 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: co 85 REAL(KIND=dp), DIMENSION(2) :: x 86 TYPE(opgrid_type), POINTER :: density 87 TYPE(opt_state_type) :: ostate 88 89 ALLOCATE (co(num_gto)) 90 co = 0._dp 91 NULLIFY (density) 92 CALL create_opgrid(density, atom%basis%grid) 93 CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, & 94 atom%state%maxl_occ, atom%state%maxn_occ) 95 CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, & 96 typ="RHO") 97 density%op = fourpi*density%op 98 IF (norder /= 0) THEN 99 density%op = density%op*atom%basis%grid%rad**norder 100 END IF 101 102 ostate%nf = 0 103 ostate%nvar = 2 104 x(1) = 0.10_dp !starting point of geometric series 105 x(2) = 2.00_dp !factor of series 106 IF (PRESENT(powell_section)) THEN 107 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) 108 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) 109 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) 110 ELSE 111 ostate%rhoend = 1.e-8_dp 112 ostate%rhobeg = 5.e-2_dp 113 ostate%maxfun = 1000 114 END IF 115 ostate%iprint = 1 116 ostate%unit = iunit 117 118 ostate%state = 0 119 IF (iunit > 0) THEN 120 WRITE (iunit, '(/," POWELL| Start optimization procedure")') 121 END IF 122 n10 = ostate%maxfun/10 123 124 DO 125 126 IF (ostate%state == 2) THEN 127 CALL density_fit(density, atom, num_gto, x(1), x(2), co, ostate%f) 128 END IF 129 130 IF (ostate%state == -1) EXIT 131 132 CALL powell_optimize(ostate%nvar, x, ostate) 133 134 IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN 135 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls")') & 136 INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp) 137 END IF 138 139 END DO 140 141 ostate%state = 8 142 CALL powell_optimize(ostate%nvar, x, ostate) 143 144 CALL release_opgrid(density) 145 146 IF (iunit > 0) THEN 147 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf 148 WRITE (iunit, '(" POWELL| Final value of function",T61,G20.10)') ostate%fopt 149 WRITE (iunit, '(" Optimized representation of density using a Geometrical Gaussian basis")') 150 WRITE (iunit, '(A,I3,/,T10,A,F10.6,T48,A,F10.6)') " Number of Gaussians:", num_gto, & 151 "Starting exponent:", x(1), "Proportionality factor:", x(2) 152 WRITE (iunit, '(A)') " Expansion coefficients" 153 WRITE (iunit, '(4F20.10)') co(1:num_gto) 154 END IF 155 156 IF (PRESENT(results)) THEN 157 CPASSERT(SIZE(results) >= num_gto + 2) 158 results(1) = x(1) 159 results(2) = x(2) 160 results(3:2 + num_gto) = co(1:num_gto) 161 END IF 162 163 DEALLOCATE (co) 164 165 END SUBROUTINE atom_fit_density 166! ************************************************************************************************** 167!> \brief ... 168!> \param atom_info ... 169!> \param basis ... 170!> \param pptype ... 171!> \param iunit output file unit 172!> \param powell_section POWELL input section 173! ************************************************************************************************** 174 SUBROUTINE atom_fit_basis(atom_info, basis, pptype, iunit, powell_section) 175 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info 176 TYPE(atom_basis_type), POINTER :: basis 177 LOGICAL, INTENT(IN) :: pptype 178 INTEGER, INTENT(IN) :: iunit 179 TYPE(section_vals_type), POINTER :: powell_section 180 181 CHARACTER(len=*), PARAMETER :: routineN = 'atom_fit_basis', routineP = moduleN//':'//routineN 182 183 INTEGER :: i, j, k, l, ll, m, n, n10, nl, nr, zval 184 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: xtob 185 LOGICAL :: explicit, mult, penalty 186 REAL(KIND=dp) :: al, crad, ear, fopt, pf, rk 187 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x 188 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: wem 189 REAL(KIND=dp), DIMENSION(:), POINTER :: w 190 TYPE(opt_state_type) :: ostate 191 192 penalty = .FALSE. 193 SELECT CASE (basis%basis_type) 194 CASE DEFAULT 195 CPABORT("") 196 CASE (GTO_BASIS) 197 IF (basis%geometrical) THEN 198 ostate%nvar = 2 199 ALLOCATE (x(2)) 200 x(1) = SQRT(basis%aval) 201 x(2) = SQRT(basis%cval) 202 ELSE 203 ll = MAXVAL(basis%nprim(:)) 204 ALLOCATE (xtob(ll, 0:lmat)) 205 xtob = 0 206 ll = SUM(basis%nprim(:)) 207 ALLOCATE (x(ll)) 208 x = 0._dp 209 ll = 0 210 DO l = 0, lmat 211 DO i = 1, basis%nprim(l) 212 mult = .FALSE. 213 DO k = 1, ll 214 IF (ABS(basis%am(i, l) - x(k)) < 1.e-6_dp) THEN 215 mult = .TRUE. 216 xtob(i, l) = k 217 END IF 218 END DO 219 IF (.NOT. mult) THEN 220 ll = ll + 1 221 x(ll) = basis%am(i, l) 222 xtob(i, l) = ll 223 END IF 224 END DO 225 END DO 226 ostate%nvar = ll 227 DO i = 1, ostate%nvar 228 x(i) = SQRT(LOG(1._dp + x(i))) 229 END DO 230 penalty = .TRUE. 231 END IF 232 CASE (STO_BASIS) 233 ll = MAXVAL(basis%nbas(:)) 234 ALLOCATE (xtob(ll, 0:lmat)) 235 xtob = 0 236 ll = SUM(basis%nbas(:)) 237 ALLOCATE (x(ll)) 238 x = 0._dp 239 ll = 0 240 DO l = 0, lmat 241 DO i = 1, basis%nbas(l) 242 ll = ll + 1 243 x(ll) = basis%as(i, l) 244 xtob(i, l) = ll 245 END DO 246 END DO 247 ostate%nvar = ll 248 DO i = 1, ostate%nvar 249 x(i) = SQRT(LOG(1._dp + x(i))) 250 END DO 251 END SELECT 252 253 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) 254 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) 255 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) 256 257 n = SIZE(atom_info, 1) 258 m = SIZE(atom_info, 2) 259 ALLOCATE (wem(n, m)) 260 wem = 1._dp 261 CALL section_vals_val_get(powell_section, "WEIGHT_ELECTRON_CONFIGURATION", explicit=explicit) 262 IF (explicit) THEN 263 CALL section_vals_val_get(powell_section, "WEIGHT_ELECTRON_CONFIGURATION", r_vals=w) 264 DO i = 1, MIN(SIZE(w), n) 265 wem(i, :) = w(i)*wem(i, :) 266 END DO 267 END IF 268 CALL section_vals_val_get(powell_section, "WEIGHT_METHOD", explicit=explicit) 269 IF (explicit) THEN 270 CALL section_vals_val_get(powell_section, "WEIGHT_METHOD", r_vals=w) 271 DO i = 1, MIN(SIZE(w), m) 272 wem(:, i) = w(i)*wem(:, i) 273 END DO 274 END IF 275 276 DO i = 1, SIZE(atom_info, 1) 277 DO j = 1, SIZE(atom_info, 2) 278 atom_info(i, j)%atom%weight = wem(i, j) 279 END DO 280 END DO 281 DEALLOCATE (wem) 282 283 ostate%nf = 0 284 ostate%iprint = 1 285 ostate%unit = iunit 286 287 ostate%state = 0 288 IF (iunit > 0) THEN 289 WRITE (iunit, '(/," POWELL| Start optimization procedure")') 290 WRITE (iunit, '(/," POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar 291 END IF 292 n10 = MAX(ostate%maxfun/100, 1) 293 294 fopt = HUGE(0._dp) 295 296 DO 297 298 IF (ostate%state == 2) THEN 299 SELECT CASE (basis%basis_type) 300 CASE DEFAULT 301 CPABORT("") 302 CASE (GTO_BASIS) 303 IF (basis%geometrical) THEN 304 basis%am = 0._dp 305 DO l = 0, lmat 306 DO i = 1, basis%nbas(l) 307 ll = i - 1 + basis%start(l) 308 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll) 309 END DO 310 END DO 311 basis%aval = x(1)*x(1) 312 basis%cval = x(2)*x(2) 313 ELSE 314 DO l = 0, lmat 315 DO i = 1, basis%nprim(l) 316 al = x(xtob(i, l))**2 317 basis%am(i, l) = EXP(al) - 1._dp 318 END DO 319 END DO 320 END IF 321 basis%bf = 0._dp 322 basis%dbf = 0._dp 323 basis%ddbf = 0._dp 324 nr = basis%grid%nr 325 DO l = 0, lmat 326 DO i = 1, basis%nbas(l) 327 al = basis%am(i, l) 328 DO k = 1, nr 329 rk = basis%grid%rad(k) 330 ear = EXP(-al*basis%grid%rad(k)**2) 331 basis%bf(k, i, l) = rk**l*ear 332 basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear 333 basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - & 334 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear 335 END DO 336 END DO 337 END DO 338 CASE (STO_BASIS) 339 DO l = 0, lmat 340 DO i = 1, basis%nbas(l) 341 al = x(xtob(i, l))**2 342 basis%as(i, l) = EXP(al) - 1._dp 343 END DO 344 END DO 345 basis%bf = 0._dp 346 basis%dbf = 0._dp 347 basis%ddbf = 0._dp 348 nr = basis%grid%nr 349 DO l = 0, lmat 350 DO i = 1, basis%nbas(l) 351 al = basis%as(i, l) 352 nl = basis%ns(i, l) 353 pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl)) 354 DO k = 1, nr 355 rk = basis%grid%rad(k) 356 ear = rk**(nl - 1)*EXP(-al*rk) 357 basis%bf(k, i, l) = pf*ear 358 basis%dbf(k, i, l) = pf*(REAL(nl - 1, dp)/rk - al)*ear 359 basis%ddbf(k, i, l) = pf*(REAL((nl - 2)*(nl - 1), dp)/rk/rk & 360 - al*REAL(2*(nl - 1), dp)/rk + al*al)*ear 361 END DO 362 END DO 363 END DO 364 END SELECT 365 CALL basis_fit(atom_info, basis, pptype, ostate%f, 0, penalty) 366 fopt = MIN(fopt, ostate%f) 367 END IF 368 369 IF (ostate%state == -1) EXIT 370 371 CALL powell_optimize(ostate%nvar, x, ostate) 372 373 IF (ostate%nf == 2 .AND. iunit > 0) THEN 374 WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f 375 END IF 376 IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN 377 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') & 378 INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt 379 END IF 380 381 END DO 382 383 ostate%state = 8 384 CALL powell_optimize(ostate%nvar, x, ostate) 385 386 IF (iunit > 0) THEN 387 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf 388 WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt 389 ! x->basis 390 SELECT CASE (basis%basis_type) 391 CASE DEFAULT 392 CPABORT("") 393 CASE (GTO_BASIS) 394 IF (basis%geometrical) THEN 395 basis%am = 0._dp 396 DO l = 0, lmat 397 DO i = 1, basis%nbas(l) 398 ll = i - 1 + basis%start(l) 399 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll) 400 END DO 401 END DO 402 basis%aval = x(1)*x(1) 403 basis%cval = x(2)*x(2) 404 ELSE 405 DO l = 0, lmat 406 DO i = 1, basis%nprim(l) 407 al = x(xtob(i, l))**2 408 basis%am(i, l) = EXP(al) - 1._dp 409 END DO 410 END DO 411 END IF 412 basis%bf = 0._dp 413 basis%dbf = 0._dp 414 basis%ddbf = 0._dp 415 nr = basis%grid%nr 416 DO l = 0, lmat 417 DO i = 1, basis%nbas(l) 418 al = basis%am(i, l) 419 DO k = 1, nr 420 rk = basis%grid%rad(k) 421 ear = EXP(-al*basis%grid%rad(k)**2) 422 basis%bf(k, i, l) = rk**l*ear 423 basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear 424 basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - & 425 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear 426 END DO 427 END DO 428 END DO 429 CASE (STO_BASIS) 430 DO l = 0, lmat 431 DO i = 1, basis%nprim(l) 432 al = x(xtob(i, l))**2 433 basis%as(i, l) = EXP(al) - 1._dp 434 END DO 435 END DO 436 basis%bf = 0._dp 437 basis%dbf = 0._dp 438 basis%ddbf = 0._dp 439 nr = basis%grid%nr 440 DO l = 0, lmat 441 DO i = 1, basis%nbas(l) 442 al = basis%as(i, l) 443 nl = basis%ns(i, l) 444 pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl)) 445 DO k = 1, nr 446 rk = basis%grid%rad(k) 447 ear = rk**(nl - 1)*EXP(-al*rk) 448 basis%bf(k, i, l) = pf*ear 449 basis%dbf(k, i, l) = pf*(REAL(nl - 1, dp)/rk - al)*ear 450 basis%ddbf(k, i, l) = pf*(REAL((nl - 2)*(nl - 1), dp)/rk/rk & 451 - al*REAL(2*(nl - 1), dp)/rk + al*al)*ear 452 END DO 453 END DO 454 END DO 455 END SELECT 456 CALL atom_print_basis(basis, iunit, " Optimized Basis") 457 CALL atom_print_basis_file(basis) 458 END IF 459 460 DEALLOCATE (x) 461 462 IF (ALLOCATED(xtob)) THEN 463 DEALLOCATE (xtob) 464 END IF 465 466 SELECT CASE (basis%basis_type) 467 CASE DEFAULT 468 CPABORT("") 469 CASE (GTO_BASIS) 470 zval = atom_info(1, 1)%atom%z 471 crad = ptable(zval)%covalent_radius*bohr 472 IF (iunit > 0) THEN 473 CALL atom_condnumber(basis, crad, iunit) 474 CALL atom_completeness(basis, zval, iunit) 475 END IF 476 CASE (STO_BASIS) 477 ! no basis test available 478 END SELECT 479 480 END SUBROUTINE atom_fit_basis 481! ************************************************************************************************** 482!> \brief ... 483!> \param atom_info ... 484!> \param atom_refs ... 485!> \param ppot ... 486!> \param iunit output file unit 487!> \param powell_section POWELL input section 488! ************************************************************************************************** 489 SUBROUTINE atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section) 490 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs 491 TYPE(atom_potential_type), POINTER :: ppot 492 INTEGER, INTENT(IN) :: iunit 493 TYPE(section_vals_type), POINTER :: powell_section 494 495 CHARACTER(len=*), PARAMETER :: routineN = 'atom_fit_pseudo', & 496 routineP = moduleN//':'//routineN 497 LOGICAL, PARAMETER :: debug = .FALSE. 498 499 CHARACTER(len=2) :: pc1, pc2, pct 500 INTEGER :: i, i1, i2, iw, j, j1, j2, k, l, m, n, & 501 n10, np, nre, nreinit, ntarget 502 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: xtob 503 INTEGER, DIMENSION(0:lmat) :: ncore 504 LOGICAL :: explicit, noopt_nlcc, preopt_nlcc 505 REAL(KIND=dp) :: afun, charge, de, deig, drho, dx, eig, fopt, oc, pchg, peig, pv, rcm, rcov, & 506 rmax, semicore_level, step_size_scaling, t_ener, t_psir0, t_semi, t_valence, t_virt, & 507 w_ener, w_node, w_psir0, w_semi, w_valence, w_virt, wtot 508 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x, xi 509 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: wem 510 REAL(KIND=dp), ALLOCATABLE, & 511 DIMENSION(:, :, :, :, :) :: dener, pval 512 REAL(KIND=dp), DIMENSION(:), POINTER :: w 513 TYPE(atom_type), POINTER :: atom 514 TYPE(opt_state_type) :: ostate 515 TYPE(wfn_init), DIMENSION(:, :), POINTER :: wfn_guess 516 517 ! weights for the optimization 518 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) 519 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) 520 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) 521 CALL section_vals_val_get(powell_section, "MAX_INIT", i_val=nreinit) 522 CALL section_vals_val_get(powell_section, "STEP_SIZE_SCALING", r_val=step_size_scaling) 523 524 CALL section_vals_val_get(powell_section, "WEIGHT_POT_VALENCE", r_val=w_valence) 525 CALL section_vals_val_get(powell_section, "WEIGHT_POT_VIRTUAL", r_val=w_virt) 526 CALL section_vals_val_get(powell_section, "WEIGHT_POT_SEMICORE", r_val=w_semi) 527 CALL section_vals_val_get(powell_section, "WEIGHT_POT_NODE", r_val=w_node) 528 CALL section_vals_val_get(powell_section, "WEIGHT_DELTA_ENERGY", r_val=w_ener) 529 530 CALL section_vals_val_get(powell_section, "TARGET_PSIR0", r_val=t_psir0) 531 CALL section_vals_val_get(powell_section, "WEIGHT_PSIR0", r_val=w_psir0) 532 CALL section_vals_val_get(powell_section, "RCOV_MULTIPLICATION", r_val=rcm) 533 534 CALL section_vals_val_get(powell_section, "TARGET_POT_VALENCE", r_val=t_valence) 535 CALL section_vals_val_get(powell_section, "TARGET_POT_VIRTUAL", r_val=t_virt) 536 CALL section_vals_val_get(powell_section, "TARGET_POT_SEMICORE", r_val=t_semi) 537 CALL section_vals_val_get(powell_section, "TARGET_DELTA_ENERGY", r_val=t_ener) 538 539 CALL section_vals_val_get(powell_section, "SEMICORE_LEVEL", r_val=semicore_level) 540 541 CALL section_vals_val_get(powell_section, "NOOPT_NLCC", l_val=noopt_nlcc) 542 CALL section_vals_val_get(powell_section, "PREOPT_NLCC", l_val=preopt_nlcc) 543 544 n = SIZE(atom_info, 1) 545 m = SIZE(atom_info, 2) 546 ALLOCATE (wem(n, m)) 547 wem = 1._dp 548 ALLOCATE (pval(8, 10, 0:lmat, m, n)) 549 ALLOCATE (dener(2, m, m, n, n)) 550 dener = 0.0_dp 551 552 ALLOCATE (wfn_guess(n, m)) 553 DO i = 1, n 554 DO j = 1, m 555 atom => atom_info(i, j)%atom 556 NULLIFY (wfn_guess(i, j)%wfn) 557 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 558 i1 = SIZE(atom%orbitals%wfn, 1) 559 i2 = SIZE(atom%orbitals%wfn, 2) 560 ALLOCATE (wfn_guess(i, j)%wfn(i1, i2, 0:lmat)) 561 END IF 562 END DO 563 END DO 564 565 CALL section_vals_val_get(powell_section, "WEIGHT_ELECTRON_CONFIGURATION", explicit=explicit) 566 IF (explicit) THEN 567 CALL section_vals_val_get(powell_section, "WEIGHT_ELECTRON_CONFIGURATION", r_vals=w) 568 DO i = 1, MIN(SIZE(w), n) 569 wem(i, :) = w(i)*wem(i, :) 570 END DO 571 END IF 572 CALL section_vals_val_get(powell_section, "WEIGHT_METHOD", explicit=explicit) 573 IF (explicit) THEN 574 CALL section_vals_val_get(powell_section, "WEIGHT_METHOD", r_vals=w) 575 DO i = 1, MIN(SIZE(w), m) 576 wem(:, i) = w(i)*wem(:, i) 577 END DO 578 END IF 579 580 IF (debug) THEN 581 CALL open_file(file_name="POWELL_RESULT", file_status="UNKNOWN", file_action="WRITE", unit_number=iw) 582 END IF 583 584 IF (ppot%gth_pot%nlcc) THEN 585 CALL opt_nlcc_param(atom_info, atom_refs, ppot%gth_pot, iunit, preopt_nlcc) 586 ELSE 587 noopt_nlcc = .TRUE. 588 preopt_nlcc = .FALSE. 589 END IF 590 591 ALLOCATE (xi(200)) 592 !decide here what to optimize 593 CALL get_pseudo_param(xi, ostate%nvar, ppot%gth_pot, noopt_nlcc) 594 ALLOCATE (x(ostate%nvar)) 595 x(1:ostate%nvar) = xi(1:ostate%nvar) 596 DEALLOCATE (xi) 597 598 ostate%nf = 0 599 ostate%iprint = 1 600 ostate%unit = iunit 601 602 ostate%state = 0 603 IF (iunit > 0) THEN 604 WRITE (iunit, '(/," POWELL| Start optimization procedure")') 605 WRITE (iunit, '(/," POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar 606 END IF 607 n10 = MAX(ostate%maxfun/100, 1) 608 609 rcov = ptable(atom_refs(1, 1)%atom%z)%covalent_radius*bohr*rcm 610 !set all reference values 611 ntarget = 0 612 wtot = 0.0_dp 613 DO i = 1, SIZE(atom_info, 1) 614 DO j = 1, SIZE(atom_info, 2) 615 atom => atom_info(i, j)%atom 616 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 617 dener(2, j, j, i, i) = atom_refs(i, j)%atom%energy%etot 618 IF (atom%state%multiplicity == -1) THEN 619 ! no spin polarization 620 atom%weight = wem(i, j) 621 ncore = get_maxn_occ(atom%state%core) 622 DO l = 0, lmat 623 np = atom%state%maxn_calc(l) 624 DO k = 1, np 625 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfn(:, ncore(l) + k, l), & 626 rcov, l, atom_refs(i, j)%atom%basis) 627 atom%orbitals%rcmax(k, l, 1) = MAX(rcov, rmax) 628 CALL atom_orbital_charge(charge, atom_refs(i, j)%atom%orbitals%wfn(:, ncore(l) + k, l), & 629 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis) 630 atom%orbitals%refene(k, l, 1) = atom_refs(i, j)%atom%orbitals%ener(ncore(l) + k, l) 631 atom%orbitals%refchg(k, l, 1) = charge 632 IF (k > atom%state%maxn_occ(l)) THEN 633 IF (k <= atom%state%maxn_occ(l) + 1) THEN 634 atom%orbitals%wrefene(k, l, 1) = w_virt 635 atom%orbitals%wrefchg(k, l, 1) = w_virt/100._dp 636 atom%orbitals%crefene(k, l, 1) = t_virt 637 atom%orbitals%reftype(k, l, 1) = "U1" 638 ntarget = ntarget + 2 639 wtot = wtot + atom%weight*(w_virt + w_virt/100._dp) 640 ELSE 641 atom%orbitals%wrefene(k, l, 1) = w_virt/100._dp 642 atom%orbitals%wrefchg(k, l, 1) = 0._dp 643 atom%orbitals%crefene(k, l, 1) = t_virt*10._dp 644 atom%orbitals%reftype(k, l, 1) = "U2" 645 ntarget = ntarget + 1 646 wtot = wtot + atom%weight*w_virt/100._dp 647 END IF 648 ELSEIF (k < atom%state%maxn_occ(l)) THEN 649 atom%orbitals%wrefene(k, l, 1) = w_semi 650 atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp 651 atom%orbitals%crefene(k, l, 1) = t_semi 652 atom%orbitals%reftype(k, l, 1) = "SC" 653 ntarget = ntarget + 2 654 wtot = wtot + atom%weight*(w_semi + w_semi/100._dp) 655 ELSE 656 IF (ABS(atom%state%occupation(l, k) - REAL(4*l + 2, KIND=dp)) < 0.01_dp .AND. & 657 ABS(atom%orbitals%refene(k, l, 1)) > semicore_level) THEN 658 ! full shell semicore 659 atom%orbitals%wrefene(k, l, 1) = w_semi 660 atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp 661 atom%orbitals%crefene(k, l, 1) = t_semi 662 atom%orbitals%reftype(k, l, 1) = "SC" 663 wtot = wtot + atom%weight*(w_semi + w_semi/100._dp) 664 ELSE 665 atom%orbitals%wrefene(k, l, 1) = w_valence 666 atom%orbitals%wrefchg(k, l, 1) = w_valence/100._dp 667 atom%orbitals%crefene(k, l, 1) = t_valence 668 atom%orbitals%reftype(k, l, 1) = "VA" 669 wtot = wtot + atom%weight*(w_valence + w_valence/100._dp) 670 END IF 671 IF (l == 0) THEN 672 atom%orbitals%tpsir0(k, 1) = t_psir0 673 atom%orbitals%wpsir0(k, 1) = w_psir0 674 wtot = wtot + atom%weight*w_psir0 675 END IF 676 ntarget = ntarget + 2 677 END IF 678 END DO 679 DO k = 1, np 680 atom%orbitals%refnod(k, l, 1) = REAL(k - 1, KIND=dp) 681 ! we only enforce 0-nodes for the first state 682 IF (k == 1 .AND. atom%state%occupation(l, k) /= 0._dp) THEN 683 atom%orbitals%wrefnod(k, l, 1) = w_node 684 wtot = wtot + atom%weight*w_node 685 END IF 686 END DO 687 END DO 688 ELSE 689 ! spin polarization 690 atom%weight = wem(i, j) 691 ncore = get_maxn_occ(atom_info(i, j)%atom%state%core) 692 DO l = 0, lmat 693 np = atom%state%maxn_calc(l) 694 DO k = 1, np 695 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfna(:, ncore(l) + k, l), & 696 rcov, l, atom_refs(i, j)%atom%basis) 697 atom%orbitals%rcmax(k, l, 1) = MAX(rcov, rmax) 698 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfnb(:, ncore(l) + k, l), & 699 rcov, l, atom_refs(i, j)%atom%basis) 700 atom%orbitals%rcmax(k, l, 2) = MAX(rcov, rmax) 701 CALL atom_orbital_charge(charge, atom_refs(i, j)%atom%orbitals%wfna(:, ncore(l) + k, l), & 702 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis) 703 atom%orbitals%refene(k, l, 1) = atom_refs(i, j)%atom%orbitals%enera(ncore(l) + k, l) 704 atom%orbitals%refchg(k, l, 1) = charge 705 CALL atom_orbital_charge(charge, atom_refs(i, j)%atom%orbitals%wfnb(:, ncore(l) + k, l), & 706 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis) 707 atom%orbitals%refene(k, l, 2) = atom_refs(i, j)%atom%orbitals%enerb(ncore(l) + k, l) 708 atom%orbitals%refchg(k, l, 2) = charge 709 ! the following assignments could be further specialized 710 IF (k > atom%state%maxn_occ(l)) THEN 711 IF (k <= atom%state%maxn_occ(l) + 1) THEN 712 atom%orbitals%wrefene(k, l, 1:2) = w_virt 713 atom%orbitals%wrefchg(k, l, 1:2) = w_virt/100._dp 714 atom%orbitals%crefene(k, l, 1:2) = t_virt 715 atom%orbitals%reftype(k, l, 1:2) = "U1" 716 ntarget = ntarget + 4 717 wtot = wtot + atom%weight*2._dp*(w_virt + w_virt/100._dp) 718 ELSE 719 atom%orbitals%wrefene(k, l, 1:2) = w_virt/100._dp 720 atom%orbitals%wrefchg(k, l, 1:2) = 0._dp 721 atom%orbitals%crefene(k, l, 1:2) = t_virt*10.0_dp 722 atom%orbitals%reftype(k, l, 1:2) = "U2" 723 wtot = wtot + atom%weight*2._dp*w_virt/100._dp 724 ntarget = ntarget + 2 725 END IF 726 ELSEIF (k < atom%state%maxn_occ(l)) THEN 727 atom%orbitals%wrefene(k, l, 1:2) = w_semi 728 atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp 729 atom%orbitals%crefene(k, l, 1:2) = t_semi 730 atom%orbitals%reftype(k, l, 1:2) = "SC" 731 ntarget = ntarget + 4 732 wtot = wtot + atom%weight*2._dp*(w_semi + w_semi/100._dp) 733 ELSE 734 IF (ABS(atom%state%occupation(l, k) - REAL(2*l + 1, KIND=dp)) < 0.01_dp .AND. & 735 ABS(atom%orbitals%refene(k, l, 1)) > semicore_level) THEN 736 atom%orbitals%wrefene(k, l, 1:2) = w_semi 737 atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp 738 atom%orbitals%crefene(k, l, 1:2) = t_semi 739 atom%orbitals%reftype(k, l, 1:2) = "SC" 740 wtot = wtot + atom%weight*2._dp*(w_semi + w_semi/100._dp) 741 ELSE 742 atom%orbitals%wrefene(k, l, 1:2) = w_valence 743 atom%orbitals%wrefchg(k, l, 1:2) = w_valence/100._dp 744 atom%orbitals%crefene(k, l, 1:2) = t_valence 745 atom%orbitals%reftype(k, l, 1:2) = "VA" 746 wtot = wtot + atom%weight*2._dp*(w_valence + w_valence/100._dp) 747 END IF 748 ntarget = ntarget + 4 749 IF (l == 0) THEN 750 atom%orbitals%tpsir0(k, 1:2) = t_psir0 751 atom%orbitals%wpsir0(k, 1:2) = w_psir0 752 wtot = wtot + atom%weight*2._dp*w_psir0 753 END IF 754 END IF 755 END DO 756 DO k = 1, np 757 atom%orbitals%refnod(k, l, 1:2) = REAL(k - 1, KIND=dp) 758 ! we only enforce 0-nodes for the first state 759 IF (k == 1 .AND. atom%state%occa(l, k) /= 0._dp) THEN 760 atom%orbitals%wrefnod(k, l, 1) = w_node 761 wtot = wtot + atom%weight*w_node 762 END IF 763 IF (k == 1 .AND. atom%state%occb(l, k) /= 0._dp) THEN 764 atom%orbitals%wrefnod(k, l, 2) = w_node 765 wtot = wtot + atom%weight*w_node 766 END IF 767 END DO 768 END DO 769 END IF 770 CALL calculate_atom(atom, 0) 771 END IF 772 END DO 773 END DO 774 ! energy differences 775 DO j1 = 1, SIZE(atom_info, 2) 776 DO j2 = 1, SIZE(atom_info, 2) 777 DO i1 = 1, SIZE(atom_info, 1) 778 DO i2 = 1, SIZE(atom_info, 1) 779 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) CYCLE 780 dener(2, j1, j2, i1, i2) = dener(2, j1, j1, i1, i1) - dener(2, j2, j2, i2, i2) 781 wtot = wtot + w_ener 782 END DO 783 END DO 784 END DO 785 END DO 786 787 DEALLOCATE (wem) 788 789 ALLOCATE (xi(ostate%nvar)) 790 DO nre = 1, nreinit 791 xi(:) = x(:) 792 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc) 793 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .TRUE.) 794 IF (nre == 1) THEN 795 WRITE (iunit, '(/," POWELL| Initial errors of target values")') 796 afun = ostate%f*wtot 797 DO i = 1, SIZE(atom_info, 1) 798 DO j = 1, SIZE(atom_info, 2) 799 atom => atom_info(i, j)%atom 800 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 801 ! start orbitals 802 wfn_guess(i, j)%wfn = atom%orbitals%wfn 803 ! 804 WRITE (iunit, '(/," Reference configuration ",T31,i5,T50," Method number ",T76,i5)') i, j 805 IF (atom%state%multiplicity == -1) THEN 806 ! no spin polarization 807 WRITE (iunit, '(" L N Occupation Eigenvalue [eV] dE [eV] dCharge ")') 808 DO l = 0, lmat 809 np = atom%state%maxn_calc(l) 810 IF (np > 0) THEN 811 DO k = 1, np 812 oc = atom%state%occupation(l, k) 813 eig = atom%orbitals%ener(k, l) 814 deig = eig - atom%orbitals%refene(k, l, 1) 815 peig = pval(1, k, l, j, i)/afun*100._dp 816 IF (pval(5, k, l, j, i) > 0.5_dp) THEN 817 pc1 = " X" 818 ELSE 819 WRITE (pc1, "(I2)") NINT(peig) 820 END IF 821 CALL atom_orbital_charge(charge, atom%orbitals%wfn(:, k, l), & 822 atom%orbitals%rcmax(k, l, 1), l, atom%basis) 823 drho = charge - atom%orbitals%refchg(k, l, 1) 824 pchg = pval(2, k, l, j, i)/afun*100._dp 825 IF (pval(6, k, l, j, i) > 0.5_dp) THEN 826 pc2 = " X" 827 ELSE 828 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 829 END IF 830 pct = atom%orbitals%reftype(k, l, 1) 831 WRITE (iunit, '(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') & 832 l, k, oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 833 END DO 834 END IF 835 END DO 836 ELSE 837 ! spin polarization 838 WRITE (iunit, '(" L N Spin Occupation Eigenvalue [eV] dE [eV] dCharge ")') 839 DO l = 0, lmat 840 np = atom%state%maxn_calc(l) 841 IF (np > 0) THEN 842 DO k = 1, np 843 oc = atom%state%occa(l, k) 844 eig = atom%orbitals%enera(k, l) 845 deig = eig - atom%orbitals%refene(k, l, 1) 846 peig = pval(1, k, l, j, i)/afun*100._dp 847 IF (pval(5, k, l, j, i) > 0.5_dp) THEN 848 pc1 = " X" 849 ELSE 850 WRITE (pc1, "(I2)") NINT(peig) 851 END IF 852 CALL atom_orbital_charge( & 853 charge, atom%orbitals%wfna(:, k, l), atom%orbitals%rcmax(k, l, 1), l, atom%basis) 854 drho = charge - atom%orbitals%refchg(k, l, 1) 855 pchg = pval(2, k, l, j, i)/afun*100._dp 856 IF (pval(6, k, l, j, i) > 0.5_dp) THEN 857 pc2 = " X" 858 ELSE 859 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 860 END IF 861 pct = atom%orbitals%reftype(k, l, 1) 862 WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') & 863 l, k, "alpha", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 864 oc = atom%state%occb(l, k) 865 eig = atom%orbitals%enerb(k, l) 866 deig = eig - atom%orbitals%refene(k, l, 2) 867 peig = pval(3, k, l, j, i)/afun*100._dp 868 IF (pval(7, k, l, j, i) > 0.5_dp) THEN 869 pc1 = " X" 870 ELSE 871 WRITE (pc1, "(I2)") NINT(peig) 872 END IF 873 CALL atom_orbital_charge( & 874 charge, atom%orbitals%wfnb(:, k, l), atom%orbitals%rcmax(k, l, 2), l, atom%basis) 875 drho = charge - atom%orbitals%refchg(k, l, 2) 876 pchg = pval(4, k, l, j, i)/afun*100._dp 877 IF (pval(8, k, l, j, i) > 0.5_dp) THEN 878 pc2 = " X" 879 ELSE 880 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 881 END IF 882 pct = atom%orbitals%reftype(k, l, 2) 883 WRITE (iunit, '(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') & 884 l, k, " beta", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 885 END DO 886 END IF 887 END DO 888 END IF 889 END IF 890 END DO 891 WRITE (iunit, *) 892 END DO 893 ! energy differences 894 IF (n*m > 1) THEN 895 WRITE (iunit, '(" Method/Econf 1 "," Method/Econf 2"," Delta Energy "," Error Energy "," Target")') 896 DO j1 = 1, SIZE(atom_info, 2) 897 DO j2 = 1, SIZE(atom_info, 2) 898 DO i1 = 1, SIZE(atom_info, 1) 899 DO i2 = 1, SIZE(atom_info, 1) 900 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) CYCLE 901 IF (ABS(dener(1, j1, j1, i1, i1)) < 0.000001_dp) CYCLE 902 IF (ABS(dener(2, j1, j1, i1, i1)) < 0.000001_dp) CYCLE 903 IF (ABS(dener(1, j2, j2, i2, i2)) < 0.000001_dp) CYCLE 904 IF (ABS(dener(2, j2, j2, i2, i2)) < 0.000001_dp) CYCLE 905 de = dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2) 906 WRITE (iunit, '(i6,i6,i10,i6,5X,F16.6,F19.6,F12.6)') & 907 j1, i1, j2, i2, dener(2, j1, j2, i1, i2), de, t_ener 908 END DO 909 END DO 910 END DO 911 END DO 912 WRITE (iunit, *) 913 END IF 914 915 WRITE (iunit, '(" Number of target values reached: ",T69,i5," of ",i3)') & 916 INT(SUM(pval(5:8, :, :, :, :))), ntarget 917 WRITE (iunit, *) 918 919 END IF 920 921 WRITE (iunit, '(" POWELL| Start optimization",I4," of total",I4,T60,"rhobeg = ",F12.8)') & 922 nre, nreinit, ostate%rhobeg 923 fopt = HUGE(0._dp) 924 ostate%state = 0 925 CALL powell_optimize(ostate%nvar, x, ostate) 926 DO 927 928 IF (ostate%state == 2) THEN 929 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc) 930 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .FALSE.) 931 fopt = MIN(fopt, ostate%f) 932 END IF 933 934 IF (ostate%state == -1) EXIT 935 936 CALL powell_optimize(ostate%nvar, x, ostate) 937 938 IF (ostate%nf == 2 .AND. iunit > 0) THEN 939 WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f 940 END IF 941 IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0 .AND. ostate%nf > 2) THEN 942 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') & 943 INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt 944 CALL put_pseudo_param(ostate%xopt, ppot%gth_pot, noopt_nlcc) 945 CALL atom_write_pseudo_param(ppot%gth_pot) 946 END IF 947 948 IF (fopt < 1.e-12_dp) EXIT 949 950 IF (debug) THEN 951 WRITE (iw, *) ostate%nf, ostate%f, x(1:ostate%nvar) 952 END IF 953 954 END DO 955 956 dx = SQRT(SUM((ostate%xopt(:) - xi(:))**2)/REAL(ostate%nvar, KIND=dp)) 957 IF (iunit > 0) THEN 958 WRITE (iunit, '(" POWELL| RMS average of variables",T69,F12.10)') dx 959 END IF 960 ostate%state = 8 961 CALL powell_optimize(ostate%nvar, x, ostate) 962 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc) 963 CALL atom_write_pseudo_param(ppot%gth_pot) 964 ! dx < SQRT(ostate%rhoend) 965 IF ((dx*dx) < ostate%rhoend) EXIT 966 ostate%rhobeg = step_size_scaling*ostate%rhobeg 967 END DO 968 DEALLOCATE (xi) 969 970 IF (iunit > 0) THEN 971 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf 972 WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt 973 974 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc) 975 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .FALSE.) 976 afun = wtot*ostate%f 977 978 WRITE (iunit, '(/," POWELL| Final errors of target values")') 979 DO i = 1, SIZE(atom_info, 1) 980 DO j = 1, SIZE(atom_info, 2) 981 atom => atom_info(i, j)%atom 982 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 983 WRITE (iunit, '(/," Reference configuration ",T31,i5,T50," Method number ",T76,i5)') i, j 984 IF (atom%state%multiplicity == -1) THEN 985 !no spin polarization 986 WRITE (iunit, '(" L N Occupation Eigenvalue [eV] dE [eV] dCharge ")') 987 DO l = 0, lmat 988 np = atom%state%maxn_calc(l) 989 IF (np > 0) THEN 990 DO k = 1, np 991 oc = atom%state%occupation(l, k) 992 eig = atom%orbitals%ener(k, l) 993 deig = eig - atom%orbitals%refene(k, l, 1) 994 peig = pval(1, k, l, j, i)/afun*100._dp 995 IF (pval(5, k, l, j, i) > 0.5_dp) THEN 996 pc1 = " X" 997 ELSE 998 WRITE (pc1, "(I2)") NINT(peig) 999 END IF 1000 CALL atom_orbital_charge( & 1001 charge, atom%orbitals%wfn(:, k, l), atom%orbitals%rcmax(k, l, 1), l, atom%basis) 1002 drho = charge - atom%orbitals%refchg(k, l, 1) 1003 pchg = pval(2, k, l, j, i)/afun*100._dp 1004 IF (pval(6, k, l, j, i) > 0.5_dp) THEN 1005 pc2 = " X" 1006 ELSE 1007 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 1008 END IF 1009 pct = atom%orbitals%reftype(k, l, 1) 1010 WRITE (iunit, '(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') & 1011 l, k, oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 1012 END DO 1013 END IF 1014 END DO 1015 np = atom%state%maxn_calc(0) 1016 DO k = 1, np 1017 CALL atom_wfnr0(pv, atom%orbitals%wfn(:, k, 0), atom%basis) 1018 IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN 1019 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN 1020 pv = 0.0_dp 1021 ELSE 1022 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 1))) 1023 END IF 1024 pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun 1025 ELSE 1026 pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp 1027 END IF 1028 WRITE (iunit, '(" s-states"," N=",I5,T40,"Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') & 1029 k, pv, NINT(pchg) 1030 END DO 1031 ELSE 1032 !spin polarization 1033 WRITE (iunit, '(" L N Spin Occupation Eigenvalue [eV] dE [eV] dCharge ")') 1034 DO l = 0, lmat 1035 np = atom%state%maxn_calc(l) 1036 IF (np > 0) THEN 1037 DO k = 1, np 1038 oc = atom%state%occa(l, k) 1039 eig = atom%orbitals%enera(k, l) 1040 deig = eig - atom%orbitals%refene(k, l, 1) 1041 peig = pval(1, k, l, j, i)/afun*100._dp 1042 IF (pval(5, k, l, j, i) > 0.5_dp) THEN 1043 pc1 = " X" 1044 ELSE 1045 WRITE (pc1, "(I2)") NINT(peig) 1046 END IF 1047 CALL atom_orbital_charge( & 1048 charge, atom%orbitals%wfna(:, k, l), atom%orbitals%rcmax(k, l, 1), l, atom%basis) 1049 drho = charge - atom%orbitals%refchg(k, l, 1) 1050 pchg = pval(2, k, l, j, i)/afun*100._dp 1051 IF (pval(6, k, l, j, i) > 0.5_dp) THEN 1052 pc2 = " X" 1053 ELSE 1054 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 1055 END IF 1056 pct = atom%orbitals%reftype(k, l, 1) 1057 WRITE (iunit, '(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') & 1058 l, k, " alpha", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 1059 oc = atom%state%occb(l, k) 1060 eig = atom%orbitals%enerb(k, l) 1061 deig = eig - atom%orbitals%refene(k, l, 2) 1062 peig = pval(3, k, l, j, i)/afun*100._dp 1063 IF (pval(7, k, l, j, i) > 0.5_dp) THEN 1064 pc1 = " X" 1065 ELSE 1066 WRITE (pc1, "(I2)") NINT(peig) 1067 END IF 1068 CALL atom_orbital_charge( & 1069 charge, atom%orbitals%wfnb(:, k, l), atom%orbitals%rcmax(k, l, 2), l, atom%basis) 1070 drho = charge - atom%orbitals%refchg(k, l, 2) 1071 pchg = pval(4, k, l, j, i)/afun*100._dp 1072 IF (pval(8, k, l, j, i) > 0.5_dp) THEN 1073 pc2 = " X" 1074 ELSE 1075 WRITE (pc2, "(I2)") MIN(NINT(pchg), 99) 1076 END IF 1077 pct = atom%orbitals%reftype(k, l, 2) 1078 WRITE (iunit, '(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') & 1079 l, k, " beta", oc, eig*evolt, pct, deig*evolt, pc1, drho, pc2 1080 END DO 1081 END IF 1082 END DO 1083 np = atom%state%maxn_calc(0) 1084 DO k = 1, np 1085 CALL atom_wfnr0(pv, atom%orbitals%wfna(:, k, 0), atom%basis) 1086 IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN 1087 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN 1088 pv = 0.0_dp 1089 ELSE 1090 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 1))) 1091 END IF 1092 pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun 1093 ELSE 1094 pchg = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp 1095 END IF 1096 WRITE (iunit, '(" s-states"," N=",I5,T35,"Alpha Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') & 1097 k, pv, NINT(pchg) 1098 ! 1099 CALL atom_wfnr0(pv, atom%orbitals%wfnb(:, k, 0), atom%basis) 1100 IF (ABS(atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp) THEN 1101 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 2))) THEN 1102 pv = 0.0_dp 1103 ELSE 1104 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 2))) 1105 END IF 1106 pchg = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv/afun 1107 ELSE 1108 pchg = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv/afun*100._dp 1109 END IF 1110 WRITE (iunit, '(" s-states"," N=",I5,T36,"Beta Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') & 1111 k, pv, NINT(pchg) 1112 END DO 1113 END IF 1114 END IF 1115 END DO 1116 END DO 1117 ! energy differences 1118 IF (n*m > 1) THEN 1119 WRITE (iunit, *) 1120 WRITE (iunit, '(" Method/Econf 1 "," Method/Econf 2"," Delta Energy "," Error Energy "," Target")') 1121 DO j1 = 1, SIZE(atom_info, 2) 1122 DO j2 = 1, SIZE(atom_info, 2) 1123 DO i1 = 1, SIZE(atom_info, 1) 1124 DO i2 = 1, SIZE(atom_info, 1) 1125 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) CYCLE 1126 IF (ABS(dener(1, j1, j1, i1, i1)) < 0.000001_dp) CYCLE 1127 IF (ABS(dener(2, j1, j1, i1, i1)) < 0.000001_dp) CYCLE 1128 IF (ABS(dener(1, j2, j2, i2, i2)) < 0.000001_dp) CYCLE 1129 IF (ABS(dener(2, j2, j2, i2, i2)) < 0.000001_dp) CYCLE 1130 de = dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2) 1131 WRITE (iunit, '(i6,i6,i10,i6,5X,F16.6,F19.6,F12.6)') j1, i1, j2, i2, dener(2, j1, j2, i1, i2), de, t_ener 1132 END DO 1133 END DO 1134 END DO 1135 END DO 1136 WRITE (iunit, *) 1137 END IF 1138 1139 WRITE (iunit, '(/," Number of target values reached: ",T69,i5," of ",i3)') INT(SUM(pval(5:8, :, :, :, :))), ntarget 1140 WRITE (iunit, *) 1141 END IF 1142 1143 DEALLOCATE (x, pval, dener) 1144 1145 DO i = 1, SIZE(wfn_guess, 1) 1146 DO j = 1, SIZE(wfn_guess, 2) 1147 IF (ASSOCIATED(wfn_guess(i, j)%wfn)) THEN 1148 DEALLOCATE (wfn_guess(i, j)%wfn) 1149 END IF 1150 END DO 1151 END DO 1152 DEALLOCATE (wfn_guess) 1153 1154 IF (ALLOCATED(xtob)) THEN 1155 DEALLOCATE (xtob) 1156 END IF 1157 1158 IF (debug) THEN 1159 CALL close_file(unit_number=iw) 1160 END IF 1161 1162 END SUBROUTINE atom_fit_pseudo 1163 1164! ************************************************************************************************** 1165!> \brief Fit NLCC density on core densities 1166!> \param atom_info ... 1167!> \param atom_refs ... 1168!> \param gthpot ... 1169!> \param iunit ... 1170!> \param preopt_nlcc ... 1171! ************************************************************************************************** 1172 SUBROUTINE opt_nlcc_param(atom_info, atom_refs, gthpot, iunit, preopt_nlcc) 1173 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs 1174 TYPE(atom_gthpot_type), INTENT(inout) :: gthpot 1175 INTEGER, INTENT(in) :: iunit 1176 LOGICAL, INTENT(in) :: preopt_nlcc 1177 1178 CHARACTER(len=*), PARAMETER :: routineN = 'opt_nlcc_param', routineP = moduleN//':'//routineN 1179 1180 INTEGER :: i, im, j, k, method 1181 REAL(KIND=dp) :: rcov, zcore, zrc, zrch 1182 TYPE(atom_type), POINTER :: aref, atom 1183 TYPE(opgrid_type), POINTER :: corden, den, den1, den2, density 1184 TYPE(opmat_type), POINTER :: denmat, dma, dmb 1185 1186 CPASSERT(gthpot%nlcc) 1187 1188 IF (iunit > 0) THEN 1189 WRITE (iunit, '(/," Core density information")') 1190 WRITE (iunit, '(A,T37,A,T55,A,T75,A)') " State Density:", "Full", "Rcov/2", "Rcov/4" 1191 END IF 1192 1193 rcov = ptable(atom_refs(1, 1)%atom%z)%covalent_radius*bohr 1194 atom => atom_info(1, 1)%atom 1195 NULLIFY (density) 1196 CALL create_opgrid(density, atom%basis%grid) 1197 density%op = 0.0_dp 1198 im = 0 1199 DO i = 1, SIZE(atom_info, 1) 1200 DO j = 1, SIZE(atom_info, 2) 1201 atom => atom_info(i, j)%atom 1202 aref => atom_refs(i, j)%atom 1203 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 1204 NULLIFY (den, denmat) 1205 CALL create_opgrid(den, atom%basis%grid) 1206 CALL create_opmat(denmat, atom%basis%nbas) 1207 1208 method = atom%method_type 1209 SELECT CASE (method) 1210 CASE (do_rks_atom, do_rhf_atom) 1211 CALL atom_denmat(denmat%op, aref%orbitals%wfn, & 1212 atom%basis%nbas, atom%state%core, & 1213 aref%state%maxl_occ, aref%state%maxn_occ) 1214 CASE (do_uks_atom, do_uhf_atom) 1215 NULLIFY (dma, dmb) 1216 CALL create_opmat(dma, atom%basis%nbas) 1217 CALL create_opmat(dmb, atom%basis%nbas) 1218 CALL atom_denmat(dma%op, aref%orbitals%wfna, & 1219 atom%basis%nbas, atom%state%core, & 1220 aref%state%maxl_occ, aref%state%maxn_occ) 1221 CALL atom_denmat(dmb%op, aref%orbitals%wfnb, & 1222 atom%basis%nbas, atom%state%core, & 1223 aref%state%maxl_occ, aref%state%maxn_occ) 1224 denmat%op = 0.5_dp*(dma%op + dmb%op) 1225 CALL release_opmat(dma) 1226 CALL release_opmat(dmb) 1227 CASE (do_rohf_atom) 1228 CPABORT("") 1229 CASE DEFAULT 1230 CPABORT("") 1231 END SELECT 1232 1233 im = im + 1 1234 CALL atom_density(den%op, denmat%op, atom%basis, aref%state%maxl_occ, typ="RHO") 1235 density%op = density%op + den%op 1236 zcore = integrate_grid(den%op, atom%basis%grid) 1237 zcore = fourpi*zcore 1238 NULLIFY (den1, den2) 1239 CALL create_opgrid(den1, atom%basis%grid) 1240 CALL create_opgrid(den2, atom%basis%grid) 1241 den1%op = 0.0_dp 1242 den2%op = 0.0_dp 1243 DO k = 1, atom%basis%grid%nr 1244 IF (atom%basis%grid%rad(k) < 0.50_dp*rcov) THEN 1245 den1%op(k) = den%op(k) 1246 END IF 1247 IF (atom%basis%grid%rad(k) < 0.25_dp*rcov) THEN 1248 den2%op(k) = den%op(k) 1249 END IF 1250 END DO 1251 zrc = integrate_grid(den1%op, atom%basis%grid) 1252 zrc = fourpi*zrc 1253 zrch = integrate_grid(den2%op, atom%basis%grid) 1254 zrch = fourpi*zrch 1255 CALL release_opgrid(den1) 1256 CALL release_opgrid(den2) 1257 CALL release_opgrid(den) 1258 CALL release_opmat(denmat) 1259 IF (iunit > 0) THEN 1260 WRITE (iunit, '(2I5,T31,F10.4,T51,F10.4,T71,F10.4)') i, j, zcore, zrc, zrch 1261 END IF 1262 END IF 1263 END DO 1264 END DO 1265 density%op = density%op/REAL(im, KIND=dp) 1266 ! 1267 IF (preopt_nlcc) THEN 1268 gthpot%nexp_nlcc = 1 1269 gthpot%nct_nlcc = 0 1270 gthpot%cval_nlcc = 0._dp 1271 gthpot%alpha_nlcc = 0._dp 1272 gthpot%nct_nlcc(1) = 1 1273 gthpot%cval_nlcc(1, 1) = 1.0_dp 1274 gthpot%alpha_nlcc(1) = gthpot%rc 1275 NULLIFY (corden) 1276 CALL create_opgrid(corden, atom%basis%grid) 1277 CALL atom_core_density(corden%op, atom%potential, typ="RHO", rr=atom%basis%grid%rad) 1278 DO k = 1, atom%basis%grid%nr 1279 IF (atom%basis%grid%rad(k) > 0.25_dp*rcov) THEN 1280 corden%op(k) = 0.0_dp 1281 END IF 1282 END DO 1283 zrc = integrate_grid(corden%op, atom%basis%grid) 1284 zrc = fourpi*zrc 1285 gthpot%cval_nlcc(1, 1) = zrch/zrc 1286 CALL release_opgrid(corden) 1287 END IF 1288 CALL release_opgrid(density) 1289 1290 END SUBROUTINE opt_nlcc_param 1291 1292! ************************************************************************************************** 1293!> \brief Low level routine to fit an atomic electron density. 1294!> \param density electron density on an atomic radial grid 'atom%basis%grid' 1295!> \param atom information about the atomic kind 1296!> (only 'atom%basis%grid' is accessed, why not to pass it instead?) 1297!> \param n number of Gaussian basis functions 1298!> \param aval exponent of the first Gaussian basis function in the series 1299!> \param cval factor of geometrical series 1300!> \param co computed expansion coefficients 1301!> \param aerr error in fitted density \f[\sqrt{\sum_{i}^{nr} (density_fitted(i)-density(i))^2 }\f] 1302! ************************************************************************************************** 1303 SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr) 1304 TYPE(opgrid_type), POINTER :: density 1305 TYPE(atom_type), POINTER :: atom 1306 INTEGER, INTENT(IN) :: n 1307 REAL(dp), INTENT(IN) :: aval, cval 1308 REAL(dp), DIMENSION(:), INTENT(INOUT) :: co 1309 REAL(dp), INTENT(OUT) :: aerr 1310 1311 CHARACTER(len=*), PARAMETER :: routineN = 'density_fit', routineP = moduleN//':'//routineN 1312 1313 INTEGER :: i, info, ip, iq, k, nr 1314 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 1315 REAL(dp) :: a, rk, zval 1316 REAL(dp), ALLOCATABLE, DIMENSION(:) :: den, pe, uval 1317 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: bf, smat, tval 1318 1319 ALLOCATE (pe(n)) 1320 nr = atom%basis%grid%nr 1321 ALLOCATE (bf(nr, n), den(nr)) 1322 bf = 0._dp 1323 1324 DO i = 1, n 1325 pe(i) = aval*cval**i 1326 a = pe(i) 1327 DO k = 1, nr 1328 rk = atom%basis%grid%rad(k) 1329 bf(k, i) = EXP(-a*rk**2) 1330 END DO 1331 END DO 1332 1333 ! total charge 1334 zval = integrate_grid(density%op, atom%basis%grid) 1335 1336 ! allocate vectors and matrices for overlaps 1337 ALLOCATE (tval(n + 1, 1), uval(n), smat(n + 1, n + 1)) 1338 DO i = 1, n 1339 uval(i) = (pi/pe(i))**1.5_dp 1340 tval(i, 1) = integrate_grid(density%op, bf(:, i), atom%basis%grid) 1341 END DO 1342 tval(n + 1, 1) = zval 1343 1344 DO iq = 1, n 1345 DO ip = 1, n 1346 smat(ip, iq) = (pi/(pe(ip) + pe(iq)))**1.5_dp 1347 END DO 1348 END DO 1349 smat(1:n, n + 1) = uval(1:n) 1350 smat(n + 1, 1:n) = uval(1:n) 1351 smat(n + 1, n + 1) = 0._dp 1352 1353 ALLOCATE (ipiv(n + 1)) 1354 CALL lapack_sgesv(n + 1, 1, smat, n + 1, ipiv, tval, n + 1, info) 1355 DEALLOCATE (ipiv) 1356 CPASSERT(info == 0) 1357 co(1:n) = tval(1:n, 1) 1358 1359 ! calculate density 1360 den(:) = 0._dp 1361 DO i = 1, n 1362 den(:) = den(:) + co(i)*bf(:, i) 1363 END DO 1364 den(:) = den(:)*fourpi 1365 den(:) = (den(:) - density%op(:))**2 1366 aerr = SQRT(integrate_grid(den, atom%basis%grid)) 1367 1368 DEALLOCATE (pe, bf, den) 1369 1370 DEALLOCATE (tval, uval, smat) 1371 1372 END SUBROUTINE density_fit 1373! ************************************************************************************************** 1374!> \brief Low level routine to fit a basis set. 1375!> \param atom_info ... 1376!> \param basis ... 1377!> \param pptype ... 1378!> \param afun ... 1379!> \param iw ... 1380!> \param penalty ... 1381! ************************************************************************************************** 1382 SUBROUTINE basis_fit(atom_info, basis, pptype, afun, iw, penalty) 1383 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info 1384 TYPE(atom_basis_type), POINTER :: basis 1385 LOGICAL, INTENT(IN) :: pptype 1386 REAL(dp), INTENT(OUT) :: afun 1387 INTEGER, INTENT(IN) :: iw 1388 LOGICAL, INTENT(IN) :: penalty 1389 1390 CHARACTER(len=*), PARAMETER :: routineN = 'basis_fit', routineP = moduleN//':'//routineN 1391 1392 INTEGER :: do_eric, do_erie, i, im, in, l, nm, nn, & 1393 reltyp, zval 1394 LOGICAL :: eri_c, eri_e 1395 REAL(KIND=dp) :: amin 1396 TYPE(atom_integrals), POINTER :: atint 1397 TYPE(atom_potential_type), POINTER :: pot 1398 TYPE(atom_type), POINTER :: atom 1399 1400 ALLOCATE (atint) 1401 1402 nn = SIZE(atom_info, 1) 1403 nm = SIZE(atom_info, 2) 1404 1405 ! calculate integrals 1406 NULLIFY (pot) 1407 zval = 0 1408 eri_c = .FALSE. 1409 eri_e = .FALSE. 1410 DO in = 1, nn 1411 DO im = 1, nm 1412 atom => atom_info(in, im)%atom 1413 IF (pptype .EQV. atom%pp_calc) THEN 1414 pot => atom%potential 1415 zval = atom_info(in, im)%atom%z 1416 reltyp = atom_info(in, im)%atom%relativistic 1417 do_eric = atom_info(in, im)%atom%coulomb_integral_type 1418 do_erie = atom_info(in, im)%atom%exchange_integral_type 1419 IF (do_eric == do_analytic) eri_c = .TRUE. 1420 IF (do_erie == do_analytic) eri_e = .TRUE. 1421 EXIT 1422 END IF 1423 END DO 1424 IF (ASSOCIATED(pot)) EXIT 1425 END DO 1426 ! general integrals 1427 CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e) 1428 ! potential 1429 CALL atom_ppint_setup(atint, basis, potential=pot) 1430 IF (pptype) THEN 1431 NULLIFY (atint%tzora, atint%hdkh) 1432 ELSE 1433 ! relativistic correction terms 1434 CALL atom_relint_setup(atint, basis, reltyp, zcore=REAL(zval, dp)) 1435 END IF 1436 1437 afun = 0._dp 1438 1439 DO in = 1, nn 1440 DO im = 1, nm 1441 atom => atom_info(in, im)%atom 1442 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 1443 IF (pptype .EQV. atom%pp_calc) THEN 1444 CALL set_atom(atom, basis=basis) 1445 CALL set_atom(atom, integrals=atint) 1446 CALL calculate_atom(atom, iw) 1447 afun = afun + atom%energy%etot*atom%weight 1448 END IF 1449 END IF 1450 END DO 1451 END DO 1452 1453 ! penalty 1454 IF (penalty) THEN 1455 DO l = 0, lmat 1456 DO i = 1, basis%nbas(l) - 1 1457 amin = MINVAL(ABS(basis%am(i, l) - basis%am(i + 1:basis%nbas(l), l))) 1458 amin = amin/basis%am(i, l) 1459 afun = afun + 10._dp*EXP(-(20._dp*amin)**4) 1460 END DO 1461 END DO 1462 END IF 1463 1464 CALL atom_int_release(atint) 1465 CALL atom_ppint_release(atint) 1466 CALL atom_relint_release(atint) 1467 1468 DEALLOCATE (atint) 1469 1470 END SUBROUTINE basis_fit 1471! ************************************************************************************************** 1472!> \brief Low level routine to fit a pseudo-potential. 1473!> \param atom_info ... 1474!> \param wfn_guess ... 1475!> \param ppot ... 1476!> \param afun ... 1477!> \param wtot ... 1478!> \param pval ... 1479!> \param dener ... 1480!> \param wen ... 1481!> \param ten ... 1482!> \param init ... 1483! ************************************************************************************************** 1484 SUBROUTINE pseudo_fit(atom_info, wfn_guess, ppot, afun, wtot, pval, dener, wen, ten, init) 1485 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info 1486 TYPE(wfn_init), DIMENSION(:, :), POINTER :: wfn_guess 1487 TYPE(atom_potential_type), POINTER :: ppot 1488 REAL(dp), INTENT(OUT) :: afun 1489 REAL(dp), INTENT(IN) :: wtot 1490 REAL(dp), DIMENSION(:, :, 0:, :, :), INTENT(INOUT) :: pval 1491 REAL(dp), DIMENSION(:, :, :, :, :), INTENT(INOUT) :: dener 1492 REAL(dp), INTENT(IN) :: wen, ten 1493 LOGICAL, INTENT(IN) :: init 1494 1495 CHARACTER(len=*), PARAMETER :: routineN = 'pseudo_fit', routineP = moduleN//':'//routineN 1496 1497 INTEGER :: i, i1, i2, j, j1, j2, k, l, n, node 1498 LOGICAL :: converged, noguess, shift 1499 REAL(KIND=dp) :: charge, de, fde, pv, rcov, rcov1, rcov2, & 1500 tv 1501 TYPE(atom_integrals), POINTER :: pp_int 1502 TYPE(atom_type), POINTER :: atom 1503 1504 afun = 0._dp 1505 pval = 0._dp 1506 dener(1, :, :, :, :) = 0._dp 1507 1508 noguess = .NOT. init 1509 1510 pp_int => atom_info(1, 1)%atom%integrals 1511 CALL atom_ppint_release(pp_int) 1512 CALL atom_ppint_setup(pp_int, atom_info(1, 1)%atom%basis, potential=ppot) 1513 1514 DO i = 1, SIZE(atom_info, 1) 1515 DO j = 1, SIZE(atom_info, 2) 1516 atom => atom_info(i, j)%atom 1517 IF (atom_consistent_method(atom%method_type, atom%state%multiplicity)) THEN 1518 IF (noguess) THEN 1519 atom%orbitals%wfn = wfn_guess(i, j)%wfn 1520 END IF 1521 CALL set_atom(atom, integrals=pp_int, potential=ppot) 1522 CALL calculate_atom(atom, 0, noguess=noguess, converged=converged) 1523 IF (.NOT. converged) THEN 1524 CALL calculate_atom(atom, 0, noguess=.FALSE., converged=shift) 1525 IF (.NOT. shift) THEN 1526 atom%orbitals%ener(:, :) = 1.5_dp*atom%orbitals%ener(:, :) + 0.5_dp 1527 END IF 1528 END IF 1529 dener(1, j, j, i, i) = atom%energy%etot 1530 DO l = 0, atom%state%maxl_calc 1531 n = atom%state%maxn_calc(l) 1532 DO k = 1, n 1533 IF (atom%state%multiplicity == -1) THEN 1534 !no spin polarization 1535 rcov = atom%orbitals%rcmax(k, l, 1) 1536 tv = atom%orbitals%crefene(k, l, 1) 1537 de = ABS(atom%orbitals%ener(k, l) - atom%orbitals%refene(k, l, 1)) 1538 fde = get_error_value(de, tv) 1539 IF (fde < 1.e-8) pval(5, k, l, j, i) = 1._dp 1540 pv = atom%weight*atom%orbitals%wrefene(k, l, 1)*fde 1541 afun = afun + pv 1542 pval(1, k, l, j, i) = pv 1543 IF (atom%orbitals%wrefchg(k, l, 1) > 0._dp) THEN 1544 CALL atom_orbital_charge(charge, atom%orbitals%wfn(:, k, l), rcov, l, atom%basis) 1545 de = ABS(charge - atom%orbitals%refchg(k, l, 1)) 1546 fde = get_error_value(de, 25._dp*tv) 1547 IF (fde < 1.e-8) pval(6, k, l, j, i) = 1._dp 1548 pv = atom%weight*atom%orbitals%wrefchg(k, l, 1)*fde 1549 afun = afun + pv 1550 pval(2, k, l, j, i) = pv 1551 END IF 1552 IF (atom%orbitals%wrefnod(k, l, 1) > 0._dp) THEN 1553 CALL atom_orbital_nodes(node, atom%orbitals%wfn(:, k, l), 2._dp*rcov, l, atom%basis) 1554 afun = afun + atom%weight*atom%orbitals%wrefnod(k, l, 1)* & 1555 ABS(REAL(node, dp) - atom%orbitals%refnod(k, l, 1)) 1556 END IF 1557 IF (l == 0) THEN 1558 IF (atom%orbitals%wpsir0(k, 1) > 0._dp) THEN 1559 CALL atom_wfnr0(pv, atom%orbitals%wfn(:, k, l), atom%basis) 1560 IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN 1561 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN 1562 pv = 0.0_dp 1563 ELSE 1564 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 1))) 1565 END IF 1566 pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv 1567 ELSE 1568 pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp 1569 END IF 1570 afun = afun + pv 1571 END IF 1572 END IF 1573 ELSE 1574 !spin polarization 1575 rcov1 = atom%orbitals%rcmax(k, l, 1) 1576 rcov2 = atom%orbitals%rcmax(k, l, 2) 1577 tv = atom%orbitals%crefene(k, l, 1) 1578 de = ABS(atom%orbitals%enera(k, l) - atom%orbitals%refene(k, l, 1)) 1579 fde = get_error_value(de, tv) 1580 IF (fde < 1.e-8) pval(5, k, l, j, i) = 1._dp 1581 pv = atom%weight*atom%orbitals%wrefene(k, l, 1)*fde 1582 afun = afun + pv 1583 pval(1, k, l, j, i) = pv 1584 tv = atom%orbitals%crefene(k, l, 2) 1585 de = ABS(atom%orbitals%enerb(k, l) - atom%orbitals%refene(k, l, 2)) 1586 fde = get_error_value(de, tv) 1587 IF (fde < 1.e-8) pval(7, k, l, j, i) = 1._dp 1588 pv = atom%weight*atom%orbitals%wrefene(k, l, 2)*fde 1589 afun = afun + pv 1590 pval(3, k, l, j, i) = pv 1591 IF (atom%orbitals%wrefchg(k, l, 1) > 0._dp) THEN 1592 CALL atom_orbital_charge(charge, atom%orbitals%wfna(:, k, l), rcov1, l, atom%basis) 1593 de = ABS(charge - atom%orbitals%refchg(k, l, 1)) 1594 fde = get_error_value(de, 25._dp*tv) 1595 IF (fde < 1.e-8) pval(6, k, l, j, i) = 1._dp 1596 pv = atom%weight*atom%orbitals%wrefchg(k, l, 1)*fde 1597 afun = afun + pv 1598 pval(2, k, l, j, i) = pv 1599 END IF 1600 IF (atom%orbitals%wrefchg(k, l, 2) > 0._dp) THEN 1601 CALL atom_orbital_charge(charge, atom%orbitals%wfnb(:, k, l), rcov2, l, atom%basis) 1602 de = ABS(charge - atom%orbitals%refchg(k, l, 2)) 1603 fde = get_error_value(de, 25._dp*tv) 1604 IF (fde < 1.e-8) pval(8, k, l, j, i) = 1._dp 1605 pv = atom%weight*atom%orbitals%wrefchg(k, l, 2)*fde 1606 afun = afun + pv 1607 pval(4, k, l, j, i) = pv 1608 END IF 1609 IF (atom%orbitals%wrefnod(k, l, 1) > 0._dp) THEN 1610 CALL atom_orbital_nodes(node, atom%orbitals%wfna(:, k, l), 2._dp*rcov1, l, atom%basis) 1611 afun = afun + atom%weight*atom%orbitals%wrefnod(k, l, 1)* & 1612 ABS(REAL(node, dp) - atom%orbitals%refnod(k, l, 1)) 1613 END IF 1614 IF (atom%orbitals%wrefnod(k, l, 2) > 0._dp) THEN 1615 CALL atom_orbital_nodes(node, atom%orbitals%wfnb(:, k, l), 2._dp*rcov2, l, atom%basis) 1616 afun = afun + atom%weight*atom%orbitals%wrefnod(k, l, 2)* & 1617 ABS(REAL(node, dp) - atom%orbitals%refnod(k, l, 2)) 1618 END IF 1619 IF (l == 0) THEN 1620 IF (atom%orbitals%wpsir0(k, 1) > 0._dp) THEN 1621 CALL atom_wfnr0(pv, atom%orbitals%wfna(:, k, l), atom%basis) 1622 IF (ABS(atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp) THEN 1623 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 1))) THEN 1624 pv = 0.0_dp 1625 ELSE 1626 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 1))) 1627 END IF 1628 pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv 1629 ELSE 1630 pv = atom%weight*atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp 1631 END IF 1632 afun = afun + pv 1633 END IF 1634 IF (atom%orbitals%wpsir0(k, 2) > 0._dp) THEN 1635 CALL atom_wfnr0(pv, atom%orbitals%wfnb(:, k, l), atom%basis) 1636 IF (ABS(atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp) THEN 1637 IF (ABS(pv) > ABS(atom%orbitals%tpsir0(k, 2))) THEN 1638 pv = 0.0_dp 1639 ELSE 1640 pv = 10._dp*(ABS(pv) - ABS(atom%orbitals%tpsir0(k, 2))) 1641 END IF 1642 pv = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv 1643 ELSE 1644 pv = atom%weight*atom%orbitals%wpsir0(k, 2)*pv*pv*100._dp 1645 END IF 1646 afun = afun + pv 1647 END IF 1648 END IF 1649 ENDIF 1650 END DO 1651 END DO 1652 END IF 1653 END DO 1654 END DO 1655 1656 ! energy differences 1657 DO j1 = 1, SIZE(atom_info, 2) 1658 DO j2 = 1, SIZE(atom_info, 2) 1659 DO i1 = 1, SIZE(atom_info, 1) 1660 DO i2 = i1 + 1, SIZE(atom_info, 1) 1661 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) CYCLE 1662 dener(1, j1, j2, i1, i2) = dener(1, j1, j1, i1, i1) - dener(1, j2, j2, i2, i2) 1663 de = ABS(dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2)) 1664 fde = get_error_value(de, ten) 1665 afun = afun + wen*fde 1666 END DO 1667 END DO 1668 END DO 1669 END DO 1670 1671 ! scaling 1672 afun = afun/wtot 1673 1674 END SUBROUTINE pseudo_fit 1675! ************************************************************************************************** 1676!> \brief Compute the squared relative error. 1677!> \param fval actual value 1678!> \param ftarget target value 1679!> \return squared relative error 1680! ************************************************************************************************** 1681 FUNCTION get_error_value(fval, ftarget) RESULT(errval) 1682 REAL(KIND=dp), INTENT(in) :: fval, ftarget 1683 REAL(KIND=dp) :: errval 1684 1685 CPASSERT(fval >= 0.0_dp) 1686 1687 IF (fval <= ftarget) THEN 1688 errval = 0.0_dp 1689 ELSE 1690 errval = (fval - ftarget)/MAX(ftarget, 1.e-10_dp) 1691 errval = errval*errval 1692 END IF 1693 1694 END FUNCTION get_error_value 1695! ************************************************************************************************** 1696!> \brief ... 1697!> \param pvec ... 1698!> \param nval ... 1699!> \param gthpot ... 1700!> \param noopt_nlcc ... 1701! ************************************************************************************************** 1702 SUBROUTINE get_pseudo_param(pvec, nval, gthpot, noopt_nlcc) 1703 REAL(KIND=dp), DIMENSION(:), INTENT(out) :: pvec 1704 INTEGER, INTENT(out) :: nval 1705 TYPE(atom_gthpot_type), INTENT(in) :: gthpot 1706 LOGICAL, INTENT(in) :: noopt_nlcc 1707 1708 CHARACTER(len=*), PARAMETER :: routineN = 'get_pseudo_param', & 1709 routineP = moduleN//':'//routineN 1710 1711 INTEGER :: i, ival, j, l, n 1712 1713 IF (gthpot%lsdpot) THEN 1714 pvec = 0 1715 ival = 0 1716 DO j = 1, gthpot%nexp_lsd 1717 ival = ival + 1 1718 pvec(ival) = rcpro(-1, gthpot%alpha_lsd(j)) 1719 DO i = 1, gthpot%nct_lsd(j) 1720 ival = ival + 1 1721 pvec(ival) = gthpot%cval_lsd(i, j) 1722 END DO 1723 END DO 1724 ELSE 1725 pvec = 0 1726 ival = 1 1727 pvec(ival) = rcpro(-1, gthpot%rc) 1728 DO i = 1, gthpot%ncl 1729 ival = ival + 1 1730 pvec(ival) = gthpot%cl(i) 1731 END DO 1732 IF (gthpot%lpotextended) THEN 1733 DO j = 1, gthpot%nexp_lpot 1734 ival = ival + 1 1735 pvec(ival) = rcpro(-1, gthpot%alpha_lpot(j)) 1736 DO i = 1, gthpot%nct_lpot(j) 1737 ival = ival + 1 1738 pvec(ival) = gthpot%cval_lpot(i, j) 1739 END DO 1740 END DO 1741 END IF 1742 IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc)) THEN 1743 DO j = 1, gthpot%nexp_nlcc 1744 ival = ival + 1 1745 pvec(ival) = rcpro(-1, gthpot%alpha_nlcc(j)) 1746 DO i = 1, gthpot%nct_nlcc(j) 1747 ival = ival + 1 1748 pvec(ival) = gthpot%cval_nlcc(i, j) 1749 END DO 1750 END DO 1751 END IF 1752 DO l = 0, lmat 1753 IF (gthpot%nl(l) > 0) THEN 1754 ival = ival + 1 1755 pvec(ival) = rcpro(-1, gthpot%rcnl(l)) 1756 END IF 1757 END DO 1758 DO l = 0, lmat 1759 IF (gthpot%nl(l) > 0) THEN 1760 n = gthpot%nl(l) 1761 DO i = 1, n 1762 DO j = i, n 1763 ival = ival + 1 1764 pvec(ival) = gthpot%hnl(i, j, l) 1765 END DO 1766 END DO 1767 END IF 1768 END DO 1769 END IF 1770 nval = ival 1771 1772 END SUBROUTINE get_pseudo_param 1773 1774! ************************************************************************************************** 1775!> \brief ... 1776!> \param pvec ... 1777!> \param gthpot ... 1778!> \param noopt_nlcc ... 1779! ************************************************************************************************** 1780 SUBROUTINE put_pseudo_param(pvec, gthpot, noopt_nlcc) 1781 REAL(KIND=dp), DIMENSION(:), INTENT(in) :: pvec 1782 TYPE(atom_gthpot_type), INTENT(inout) :: gthpot 1783 LOGICAL, INTENT(in) :: noopt_nlcc 1784 1785 CHARACTER(len=*), PARAMETER :: routineN = 'put_pseudo_param', & 1786 routineP = moduleN//':'//routineN 1787 1788 INTEGER :: i, ival, j, l, n 1789 1790 IF (gthpot%lsdpot) THEN 1791 ival = 0 1792 DO j = 1, gthpot%nexp_lsd 1793 ival = ival + 1 1794 gthpot%alpha_lsd(j) = rcpro(1, pvec(ival)) 1795 DO i = 1, gthpot%nct_lsd(j) 1796 ival = ival + 1 1797 gthpot%cval_lsd(i, j) = pvec(ival) 1798 END DO 1799 END DO 1800 ELSE 1801 ival = 1 1802 gthpot%rc = rcpro(1, pvec(ival)) 1803 DO i = 1, gthpot%ncl 1804 ival = ival + 1 1805 gthpot%cl(i) = pvec(ival) 1806 END DO 1807 IF (gthpot%lpotextended) THEN 1808 DO j = 1, gthpot%nexp_lpot 1809 ival = ival + 1 1810 gthpot%alpha_lpot(j) = rcpro(1, pvec(ival)) 1811 DO i = 1, gthpot%nct_lpot(j) 1812 ival = ival + 1 1813 gthpot%cval_lpot(i, j) = pvec(ival) 1814 END DO 1815 END DO 1816 END IF 1817 IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc)) THEN 1818 DO j = 1, gthpot%nexp_nlcc 1819 ival = ival + 1 1820 gthpot%alpha_nlcc(j) = rcpro(1, pvec(ival)) 1821 DO i = 1, gthpot%nct_nlcc(j) 1822 ival = ival + 1 1823 gthpot%cval_nlcc(i, j) = pvec(ival) 1824 END DO 1825 END DO 1826 END IF 1827 DO l = 0, lmat 1828 IF (gthpot%nl(l) > 0) THEN 1829 ival = ival + 1 1830 gthpot%rcnl(l) = rcpro(1, pvec(ival)) 1831 END IF 1832 END DO 1833 DO l = 0, lmat 1834 IF (gthpot%nl(l) > 0) THEN 1835 n = gthpot%nl(l) 1836 DO i = 1, n 1837 DO j = i, n 1838 ival = ival + 1 1839 gthpot%hnl(i, j, l) = pvec(ival) 1840 END DO 1841 END DO 1842 END IF 1843 END DO 1844 END IF 1845 1846 END SUBROUTINE put_pseudo_param 1847 1848! ************************************************************************************************** 1849!> \brief Transform xval according to the following rules: 1850!> direct (id == +1): yval = 2 tanh^{2}(xval / 10) 1851!> inverse (id == -1): yval = 10 arctanh(\sqrt{xval/2}) 1852!> \param id direction of the transformation 1853!> \param xval value to transform 1854!> \return transformed value 1855! ************************************************************************************************** 1856 FUNCTION rcpro(id, xval) RESULT(yval) 1857 INTEGER, INTENT(IN) :: id 1858 REAL(dp), INTENT(IN) :: xval 1859 REAL(dp) :: yval 1860 1861 CHARACTER(len=*), PARAMETER :: routineN = 'rcpro', routineP = moduleN//':'//routineN 1862 1863 REAL(dp) :: x1, x2 1864 1865 IF (id == 1) THEN 1866 yval = 2.0_dp*TANH(0.1_dp*xval)**2 1867 ELSE IF (id == -1) THEN 1868 x1 = SQRT(xval/2.0_dp) 1869 CPASSERT(x1 <= 1._dp) 1870 x2 = 0.5_dp*LOG((1._dp + x1)/(1._dp - x1)) 1871 yval = x2/0.1_dp 1872 ELSE 1873 CPABORT("wrong id") 1874 END IF 1875 END FUNCTION rcpro 1876 1877! ************************************************************************************************** 1878!> \brief ... 1879!> \param atom ... 1880!> \param num_gau ... 1881!> \param num_pol ... 1882!> \param iunit ... 1883!> \param powell_section ... 1884!> \param results ... 1885! ************************************************************************************************** 1886 SUBROUTINE atom_fit_kgpot(atom, num_gau, num_pol, iunit, powell_section, results) 1887 TYPE(atom_type), POINTER :: atom 1888 INTEGER, INTENT(IN) :: num_gau, num_pol, iunit 1889 TYPE(section_vals_type), OPTIONAL, POINTER :: powell_section 1890 REAL(KIND=dp), DIMENSION(:), OPTIONAL :: results 1891 1892 CHARACTER(len=*), PARAMETER :: routineN = 'atom_fit_kgpot', routineP = moduleN//':'//routineN 1893 REAL(KIND=dp), PARAMETER :: t23 = 2._dp/3._dp 1894 1895 INTEGER :: i, ig, ip, iw, j, n10 1896 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x 1897 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: co 1898 TYPE(opgrid_type), POINTER :: density 1899 TYPE(opt_state_type) :: ostate 1900 1901! at least one parameter to be optimized 1902 1903 CPASSERT(num_pol*num_gau > 0) 1904 1905 ALLOCATE (co(num_pol + 1, num_gau), x(num_pol*num_gau + num_gau)) 1906 co = 0._dp 1907 1908 ! calculate density 1909 NULLIFY (density) 1910 CALL create_opgrid(density, atom%basis%grid) 1911 CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, & 1912 atom%state%maxl_occ, atom%state%maxn_occ) 1913 CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO") 1914 ! target functional 1915 density%op = t23*(0.3_dp*(3.0_dp*pi*pi)**t23)*density%op**t23 1916 1917 ! initiallize parameter 1918 ostate%nf = 0 1919 ostate%nvar = num_pol*num_gau + num_gau 1920 DO i = 1, num_gau 1921 co(1, i) = 0.5_dp + REAL(i - 1, KIND=dp) 1922 co(2, i) = 1.0_dp 1923 DO j = 3, num_pol + 1 1924 co(j, i) = 0.1_dp 1925 END DO 1926 END DO 1927 CALL putvar(x, co, num_pol, num_gau) 1928 1929 IF (PRESENT(powell_section)) THEN 1930 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend) 1931 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg) 1932 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun) 1933 ELSE 1934 ostate%rhoend = 1.e-8_dp 1935 ostate%rhobeg = 5.e-2_dp 1936 ostate%maxfun = 1000 1937 END IF 1938 ostate%iprint = 1 1939 ostate%unit = iunit 1940 1941 ostate%state = 0 1942 IF (iunit > 0) THEN 1943 WRITE (iunit, '(/," ",13("*")," Approximated Nonadditive Kinetic Energy Functional ",14("*"))') 1944 WRITE (iunit, '(" POWELL| Start optimization procedure")') 1945 END IF 1946 n10 = 50 1947 1948 DO 1949 1950 IF (ostate%state == 2) THEN 1951 CALL getvar(x, co, num_pol, num_gau) 1952 CALL kgpot_fit(density, num_gau, num_pol, co, ostate%f) 1953 END IF 1954 1955 IF (ostate%state == -1) EXIT 1956 1957 CALL powell_optimize(ostate%nvar, x, ostate) 1958 1959 IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN 1960 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T66,G15.7)') & 1961 INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), ostate%fopt 1962 END IF 1963 1964 END DO 1965 1966 ostate%state = 8 1967 CALL powell_optimize(ostate%nvar, x, ostate) 1968 CALL getvar(x, co, num_pol, num_gau) 1969 1970 CALL release_opgrid(density) 1971 1972 IF (iunit > 0) THEN 1973 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf 1974 WRITE (iunit, '(" POWELL| Final value of function",T61,G20.10)') ostate%fopt 1975 WRITE (iunit, '(" Optimized local potential of approximated nonadditive kinetic energy functional")') 1976 DO ig = 1, num_gau 1977 WRITE (iunit, '(I2,T15,"Gaussian polynomial expansion",T66,"Rc=",F12.4)') ig, co(1, ig) 1978 WRITE (iunit, '(T15,"Coefficients",T33,4F12.4)') (co(1 + ip, ig), ip=1, num_pol) 1979 END DO 1980 END IF 1981 1982 CALL open_file(file_name="FIT_RESULT", file_status="UNKNOWN", file_action="WRITE", unit_number=iw) 1983 WRITE (iw, *) ptable(atom%z)%symbol 1984 WRITE (iw, *) num_gau, num_pol 1985 DO ig = 1, num_gau 1986 WRITE (iw, '(T10,F12.4,6X,4F12.4)') (co(ip, ig), ip=1, num_pol + 1) 1987 END DO 1988 CALL close_file(unit_number=iw) 1989 1990 IF (PRESENT(results)) THEN 1991 CPASSERT(SIZE(results) >= SIZE(x)) 1992 results = x 1993 END IF 1994 1995 DEALLOCATE (co, x) 1996 1997 END SUBROUTINE atom_fit_kgpot 1998 1999! ************************************************************************************************** 2000!> \brief ... 2001!> \param kgpot ... 2002!> \param ng ... 2003!> \param np ... 2004!> \param cval ... 2005!> \param aerr ... 2006! ************************************************************************************************** 2007 SUBROUTINE kgpot_fit(kgpot, ng, np, cval, aerr) 2008 TYPE(opgrid_type), POINTER :: kgpot 2009 INTEGER, INTENT(IN) :: ng, np 2010 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: cval 2011 REAL(dp), INTENT(OUT) :: aerr 2012 2013 CHARACTER(len=*), PARAMETER :: routineN = 'kgpot_fit', routineP = moduleN//':'//routineN 2014 2015 INTEGER :: i, ig, ip, n 2016 REAL(KIND=dp) :: pc, rc 2017 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: pval 2018 2019 n = kgpot%grid%nr 2020 ALLOCATE (pval(n)) 2021 pval = 0.0_dp 2022 DO i = 1, n 2023 DO ig = 1, ng 2024 rc = kgpot%grid%rad(i)/cval(1, ig) 2025 pc = 0.0_dp 2026 DO ip = 1, np 2027 pc = pc + cval(ip + 1, ig)*rc**(2*ip - 2) 2028 END DO 2029 pval(i) = pval(i) + pc*EXP(-0.5_dp*rc*rc) 2030 END DO 2031 END DO 2032 pval(1:n) = (pval(1:n) - kgpot%op(1:n))**2 2033 aerr = fourpi*SUM(pval(1:n)*kgpot%grid%wr(1:n)) 2034 2035 DEALLOCATE (pval) 2036 2037 END SUBROUTINE kgpot_fit 2038 2039! ************************************************************************************************** 2040!> \brief ... 2041!> \param xvar ... 2042!> \param cvar ... 2043!> \param np ... 2044!> \param ng ... 2045! ************************************************************************************************** 2046 PURE SUBROUTINE getvar(xvar, cvar, np, ng) 2047 REAL(KIND=dp), DIMENSION(:), INTENT(in) :: xvar 2048 REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: cvar 2049 INTEGER, INTENT(IN) :: np, ng 2050 2051 INTEGER :: ig, ii, ip 2052 2053 ii = 0 2054 DO ig = 1, ng 2055 ii = ii + 1 2056 cvar(1, ig) = xvar(ii) 2057 DO ip = 1, np 2058 ii = ii + 1 2059 cvar(ip + 1, ig) = xvar(ii)**2 2060 END DO 2061 END DO 2062 2063 END SUBROUTINE getvar 2064 2065! ************************************************************************************************** 2066!> \brief ... 2067!> \param xvar ... 2068!> \param cvar ... 2069!> \param np ... 2070!> \param ng ... 2071! ************************************************************************************************** 2072 PURE SUBROUTINE putvar(xvar, cvar, np, ng) 2073 REAL(KIND=dp), DIMENSION(:), INTENT(inout) :: xvar 2074 REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: cvar 2075 INTEGER, INTENT(IN) :: np, ng 2076 2077 INTEGER :: ig, ii, ip 2078 2079 ii = 0 2080 DO ig = 1, ng 2081 ii = ii + 1 2082 xvar(ii) = cvar(1, ig) 2083 DO ip = 1, np 2084 ii = ii + 1 2085 xvar(ii) = SQRT(ABS(cvar(ip + 1, ig))) 2086 END DO 2087 END DO 2088 2089 END SUBROUTINE putvar 2090 2091END MODULE atom_fit 2092