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