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