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