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