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