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