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