(git:0de0cc2)
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 ! **************************************************************************************************
11 MODULE atom_fit
13  USE atom_operators, ONLY: atom_int_release,&
19  USE atom_output, ONLY: atom_print_basis,&
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, &
26  USE atom_utils, ONLY: &
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,&
38  USE input_section_types, ONLY: section_vals_type,&
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,&
50 #include "./base/base_uses.f90"
51 
52  IMPLICIT NONE
53 
54  PRIVATE
55 
56  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_fit'
57 
59 
60  TYPE wfn_init
61  REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn => null()
62  END TYPE wfn_init
63 
64 CONTAINS
65 
66 ! **************************************************************************************************
67 !> \brief Fit the atomic electron density using a geometrical Gaussian basis set.
68 !> \param atom information about the atomic kind
69 !> \param num_gto number of Gaussian basis functions
70 !> \param norder ...
71 !> \param iunit output file unit
72 !> \param powell_section POWELL input section
73 !> \param results (output) array of num_gto+2 real numbers in the following format:
74 !> starting exponent, factor of geometrical series, expansion coefficients(1:num_gto)
75 ! **************************************************************************************************
76  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 
2068 END MODULE atom_fit
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
routines that fit parameters for /from atomic calculations
Definition: atom_fit.F:11
subroutine, public atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
...
Definition: atom_fit.F:488
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
subroutine, public atom_fit_kgpot(atom, num_gau, num_pol, iunit, powell_section, results)
...
Definition: atom_fit.F:1867
subroutine, public atom_fit_basis(atom_info, basis, pptype, iunit, powell_section)
...
Definition: atom_fit.F:175
Calculate the atomic operator matrices.
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Routines that print various information about an atomic kind.
Definition: atom_output.F:11
subroutine, public atom_print_basis(atom_basis, iw, title)
Print atomic basis set.
Definition: atom_output.F:292
subroutine, public atom_print_basis_file(atom_basis, wfn)
Print the optimized atomic basis set into a file.
Definition: atom_output.F:385
subroutine, public atom_write_pseudo_param(gthpot, iunit)
Print GTH pseudo-potential parameters.
Definition: atom_output.F:744
Define the atom type and its sub types.
Definition: atom_types.F:15
integer, parameter, public gto_basis
Definition: atom_types.F:69
integer, parameter, public sto_basis
Definition: atom_types.F:69
subroutine, public release_opmat(opmat)
...
Definition: atom_types.F:1368
subroutine, public release_opgrid(opgrid)
...
Definition: atom_types.F:1408
integer, parameter, public lmat
Definition: atom_types.F:67
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
subroutine, public create_opmat(opmat, n, lmax)
...
Definition: atom_types.F:1340
subroutine, public create_opgrid(opgrid, grid)
...
Definition: atom_types.F:1385
Some basic routines for atomic calculations.
Definition: atom_utils.F:15
subroutine, public atom_condnumber(basis, crad, iw)
Print condition numbers of the given atomic basis set.
Definition: atom_utils.F:2310
pure subroutine, public atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
Calculate a density matrix using atomic orbitals.
Definition: atom_utils.F:328
pure subroutine, public atom_orbital_charge(charge, wfn, rcov, l, basis)
...
Definition: atom_utils.F:654
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
subroutine, public atom_core_density(corden, potential, typ, rr)
...
Definition: atom_utils.F:694
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
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition: atom_utils.F:367
pure subroutine, public atom_orbital_max(rmax, wfn, rcov, l, basis)
...
Definition: atom_utils.F:843
subroutine, public atom_completeness(basis, zv, iw)
Estimate completeness of the given atomic basis set.
Definition: atom_utils.F:2345
pure subroutine, public atom_orbital_nodes(node, wfn, rcov, l, basis)
...
Definition: atom_utils.F:882
pure subroutine, public atom_wfnr0(value, wfn, basis)
...
Definition: atom_utils.F:917
Definition: atom.F:9
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
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
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
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_rhf_atom
integer, parameter, public do_rks_atom
integer, parameter, public do_analytic
integer, parameter, public do_uhf_atom
integer, parameter, public do_uks_atom
integer, parameter, public do_rohf_atom
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Interface to the LAPACK F77 library.
Definition: lapack.F:17
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
Definition: powell.F:9
subroutine, public powell_optimize(n, x, optstate)
...
Definition: powell.F:55