(git:b279b6b)
atom_electronic_structure.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 MODULE atom_electronic_structure
9  USE atom_optimization, ONLY: atom_history_init,&
10  atom_history_release,&
11  atom_history_type,&
12  atom_history_update,&
13  atom_opt_fmat
14  USE atom_output, ONLY: atom_print_energies,&
15  atom_print_iteration,&
16  atom_print_state,&
17  atom_print_zmp_iteration
18  USE atom_types, ONLY: &
19  atom_type, create_opgrid, create_opmat, ecp_pseudo, gth_pseudo, lmat, no_pseudo, &
20  opgrid_type, opmat_type, release_opgrid, release_opmat, set_atom, setup_hf_section, &
21  sgp_pseudo, upf_pseudo
22  USE atom_utils, ONLY: &
23  atom_denmat, atom_density, atom_read_external_density, atom_read_external_vxc, &
24  atom_read_zmp_restart, atom_solve, atom_trace, ceri_contract, coulomb_potential_analytic, &
25  coulomb_potential_numeric, eeri_contract, err_matrix, exchange_numeric, &
26  exchange_semi_analytic, numpot_matrix, slater_density, wigner_slater_functional
27  USE atom_xc, ONLY: calculate_atom_ext_vxc,&
28  calculate_atom_vxc_lda,&
29  calculate_atom_vxc_lsd,&
30  calculate_atom_zmp
31  USE input_constants, ONLY: &
32  do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, do_nonrel_atom, &
33  do_numeric, do_rhf_atom, do_rks_atom, do_rohf_atom, do_sczoramp_atom, do_semi_analytic, &
34  do_uhf_atom, do_uks_atom, do_zoramp_atom
35  USE input_section_types, ONLY: section_vals_get,&
36  section_vals_get_subs_vals,&
37  section_vals_type
38  USE kinds, ONLY: dp
39 #include "./base/base_uses.f90"
40 
41  IMPLICIT NONE
42  PRIVATE
43  PUBLIC :: calculate_atom
44 
45  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_electronic_structure'
46 
47 CONTAINS
48 
49 ! **************************************************************************************************
50 !> \brief General routine to perform electronic structure atomic calculations.
51 !> \param atom information about the atomic kind
52 !> \param iw output file unit
53 !> \param noguess skip initial guess
54 !> \param converged whether SCF iterations have been converged
55 !> \par History
56 !> * 06.2017 disable XC input options [Juerg Hutter]
57 !> * 11.2009 created from the subroutine calculate_atom() [Juerg Hutter]
58 ! **************************************************************************************************
59  SUBROUTINE calculate_atom(atom, iw, noguess, converged)
60  TYPE(atom_type), POINTER :: atom
61  INTEGER, INTENT(IN) :: iw
62  LOGICAL, INTENT(IN), OPTIONAL :: noguess
63  LOGICAL, INTENT(OUT), OPTIONAL :: converged
64 
65  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom'
66 
67  INTEGER :: handle, method
68  LOGICAL :: explicit
69  TYPE(section_vals_type), POINTER :: sub_section
70 
71  CALL timeset(routinen, handle)
72 
73  ! test for not supported methods and functionals
74  IF (ASSOCIATED(atom%xc_section)) THEN
75  !
76  sub_section => section_vals_get_subs_vals(atom%xc_section, "ADIABATIC_RESCALING")
77  CALL section_vals_get(sub_section, explicit=explicit)
78  IF (explicit) CALL cp_abort(__location__, "ADIABATIC_RESCALING not supported in ATOM code")
79  !
80  sub_section => section_vals_get_subs_vals(atom%xc_section, "VDW_POTENTIAL")
81  CALL section_vals_get(sub_section, explicit=explicit)
82  IF (explicit) CALL cp_abort(__location__, "VDW_POTENTIAL not supported in ATOM code")
83  !
84  sub_section => section_vals_get_subs_vals(atom%xc_section, "XC_POTENTIAL")
85  CALL section_vals_get(sub_section, explicit=explicit)
86  IF (explicit) CALL cp_abort(__location__, "XC_POTENTIAL not supported in ATOM code")
87  !
88  sub_section => section_vals_get_subs_vals(atom%xc_section, "WF_CORRELATION")
89  CALL section_vals_get(sub_section, explicit=explicit)
90  IF (explicit) CALL cp_abort(__location__, "WF_CORRELATION methods not supported in ATOM code")
91  !
92  END IF
93 
94  method = atom%method_type
95  SELECT CASE (method)
96  CASE (do_rks_atom, do_rhf_atom)
97  CALL calculate_atom_restricted(atom, iw, noguess, converged)
98  CASE (do_uks_atom, do_uhf_atom)
99  CALL calculate_atom_unrestricted(atom, iw, noguess, converged)
100  CASE (do_rohf_atom)
101  cpabort("")
102  CASE DEFAULT
103  cpabort("")
104  END SELECT
105 
106  CALL timestop(handle)
107 
108  END SUBROUTINE calculate_atom
109 
110 ! **************************************************************************************************
111 !> \brief Perform restricted (closed shell) electronic structure atomic calculations.
112 !> \param atom information about the atomic kind
113 !> \param iw output file unit
114 !> \param noguess skip initial guess
115 !> \param converged whether SCF iterations have been converged
116 !> \par History
117 !> * 02.2016 support UPF files and ECP potentials [Juerg Hutter]
118 !> * 09.2015 Polarized Atomic Orbitals (PAO) method [Ole Schuett]
119 !> * 11.2013 Zhao, Morrison, and Parr (ZMP) potential [Daniele Varsano]
120 !> * 07.2013 scaled ZORA with model potential [Juerg Hutter]
121 !> * 11.2009 split into three subroutines calculate_atom(), calculate_atom_restricted(),
122 !> and calculate_atom_unrestricted() [Juerg Hutter]
123 !> * 12.2008 refactored and renamed to calculate_atom() [Juerg Hutter]
124 !> * 08.2008 created as atom_electronic_structure() [Juerg Hutter]
125 ! **************************************************************************************************
126  SUBROUTINE calculate_atom_restricted(atom, iw, noguess, converged)
127  TYPE(atom_type), POINTER :: atom
128  INTEGER, INTENT(IN) :: iw
129  LOGICAL, INTENT(IN), OPTIONAL :: noguess
130  LOGICAL, INTENT(OUT), OPTIONAL :: converged
131 
132  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom_restricted'
133 
134  INTEGER :: handle, i, iter, l, max_iter, method, &
135  reltyp
136  LOGICAL :: do_hfx, doguess, is_converged, need_vxc, &
137  need_x, need_xc, need_zmp
138  REAL(kind=dp) :: deps, eps_scf, hf_frac
139  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ext_density, ext_vxc, tmp_dens
140  TYPE(atom_history_type) :: history
141  TYPE(opgrid_type), POINTER :: cpot, density
142  TYPE(opmat_type), POINTER :: fmat, hcore, jmat, kmat, xcmat
143  TYPE(section_vals_type), POINTER :: xc_section
144 
145  CALL timeset(routinen, handle)
146 
147  IF (PRESENT(noguess)) THEN
148  doguess = .NOT. noguess
149  ELSE
150  doguess = .true.
151  END IF
152 
153  CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
154 
155  method = atom%method_type
156  max_iter = atom%optimization%max_iter
157  eps_scf = atom%optimization%eps_scf
158 
159  SELECT CASE (method)
160  CASE DEFAULT
161  cpabort("")
162  CASE (do_rks_atom)
163  need_x = do_hfx
164  need_xc = .true.
165  CASE (do_uks_atom)
166  cpabort("")
167  CASE (do_rhf_atom)
168  need_x = .true.
169  need_xc = .false.
170  hf_frac = 1._dp
171  CASE (do_uhf_atom)
172  cpabort("")
173  CASE (do_rohf_atom)
174  need_x = .true.
175  need_xc = .false.
176  hf_frac = 1._dp
177  cpabort("")
178  END SELECT
179 
180  ! ZMP starting to read external density for the zmp calculation
181  need_zmp = atom%do_zmp
182  IF (need_zmp) THEN
183  need_x = .false.
184  need_xc = .false.
185  ALLOCATE (ext_density(atom%basis%grid%nr))
186  ext_density = 0._dp
187  CALL atom_read_external_density(ext_density, atom, iw)
188  END IF
189 
190  ! ZMP starting to read external vxc for the zmp calculation
191  need_vxc = atom%read_vxc
192 
193  IF (need_vxc) THEN
194  need_x = .false.
195  need_xc = .false.
196  need_zmp = .false.
197  ALLOCATE (ext_vxc(atom%basis%grid%nr))
198  ext_vxc = 0._dp
199  CALL atom_read_external_vxc(ext_vxc, atom, iw)
200  END IF
201 
202  ! check for relativistic method
203  reltyp = atom%relativistic
204 
205  IF (iw > 0) CALL atom_print_state(atom%state, iw)
206 
207  NULLIFY (hcore)
208  CALL create_opmat(hcore, atom%basis%nbas)
209  ! Pseudopotentials
210  SELECT CASE (atom%potential%ppot_type)
211  CASE DEFAULT
212  cpabort("")
213  CASE (no_pseudo)
214  SELECT CASE (reltyp)
215  CASE DEFAULT
216  cpabort("")
217  CASE (do_nonrel_atom)
218  hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
219  CASE (do_zoramp_atom, do_sczoramp_atom)
220  hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
221  CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
222  hcore%op = atom%integrals%hdkh
223  END SELECT
224  CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
225  hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
226  END SELECT
227  ! add confinement potential (not included in relativistic transformations)
228  IF (atom%potential%confinement) THEN
229  hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
230  END IF
231 
232  NULLIFY (fmat, jmat, kmat, xcmat)
233  CALL create_opmat(fmat, atom%basis%nbas)
234  CALL create_opmat(jmat, atom%basis%nbas)
235  CALL create_opmat(kmat, atom%basis%nbas)
236  CALL create_opmat(xcmat, atom%basis%nbas)
237 
238  NULLIFY (density, cpot)
239  CALL create_opgrid(density, atom%basis%grid)
240  CALL create_opgrid(cpot, atom%basis%grid)
241 
242  ! ZMP reading the file to restart
243  IF (atom%doread) CALL atom_read_zmp_restart(atom, doguess, iw)
244 
245  IF (doguess) THEN
246  ! initial guess
247  ALLOCATE (tmp_dens(SIZE(density%op)))
248  CALL slater_density(density%op, tmp_dens, atom%z, atom%state, atom%basis%grid)
249  density%op = density%op + tmp_dens
250  DEALLOCATE (tmp_dens)
251  CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
252  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
253  CALL wigner_slater_functional(density%op, cpot%op)
254  CALL numpot_matrix(xcmat%op, cpot%op, atom%basis, 0)
255  fmat%op = hcore%op + jmat%op + xcmat%op
256  CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
257  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
258  END IF
259  CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
260  atom%state%maxl_occ, atom%state%maxn_occ)
261 
262  ! wavefunction history
263  NULLIFY (history%dmat, history%hmat)
264  CALL atom_history_init(history, atom%optimization, fmat%op)
265  is_converged = .false.
266  iter = 0
267  DO !SCF Loop
268 
269  ! Kinetic energy
270  atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
271 
272  ! Band energy
273  atom%energy%eband = 0._dp
274  DO l = 0, lmat
275  DO i = 1, min(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
276  atom%energy%eband = atom%energy%eband + atom%orbitals%ener(i, l)*atom%state%occupation(l, i)
277  END DO
278  END DO
279 
280  ! Pseudopotential energy
281  SELECT CASE (atom%potential%ppot_type)
282  CASE DEFAULT
283  cpabort("")
284  CASE (no_pseudo)
285  atom%energy%eploc = 0._dp
286  atom%energy%epnl = 0._dp
287  CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
288  atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
289  atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
290  END SELECT
291  atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
292 
293  ! Core energy
294  atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
295 
296  ! Confinement energy
297  IF (atom%potential%confinement) THEN
298  atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
299  ELSE
300  atom%energy%econfinement = 0._dp
301  END IF
302 
303  ! Hartree Term
304  jmat%op = 0._dp
305  SELECT CASE (atom%coulomb_integral_type)
306  CASE DEFAULT
307  cpabort("")
308  CASE (do_analytic)
309  CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
310  CASE (do_semi_analytic)
311  CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
312  atom%state%maxl_occ)
313  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
314  CASE (do_numeric)
315  CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
316  CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
317  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
318  END SELECT
319  atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
320 
321  ! Exchange Term
322  IF (need_x) THEN
323  kmat%op = 0._dp
324  SELECT CASE (atom%exchange_integral_type)
325  CASE DEFAULT
326  cpabort("")
327  CASE (do_analytic)
328  CALL eeri_contract(kmat%op, atom%integrals%eeri, atom%orbitals%pmat, atom%integrals%n)
329  CASE (do_semi_analytic)
330  CALL exchange_semi_analytic(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
331  CASE (do_numeric)
332  CALL exchange_numeric(kmat%op, atom%state, atom%state%occupation, atom%orbitals%wfn, atom%basis, atom%hfx_pot)
333  END SELECT
334  atom%energy%eexchange = hf_frac*0.5_dp*atom_trace(kmat%op, atom%orbitals%pmat)
335  kmat%op = hf_frac*kmat%op
336  ELSE
337  kmat%op = 0._dp
338  atom%energy%eexchange = 0._dp
339  END IF
340 
341  ! XC
342  IF (need_xc) THEN
343  xcmat%op = 0._dp
344  CALL calculate_atom_vxc_lda(xcmat, atom, xc_section)
345  ! ZMP added options for the zmp calculations, building external density and vxc potential
346  ELSEIF (need_zmp) THEN
347  xcmat%op = 0._dp
348  CALL calculate_atom_zmp(ext_density=ext_density, atom=atom, lprint=.false., xcmat=xcmat)
349  ELSEIF (need_vxc) THEN
350  xcmat%op = 0._dp
351  CALL calculate_atom_ext_vxc(vxc=ext_vxc, atom=atom, lprint=.false., xcmat=xcmat)
352  ELSE
353  xcmat%op = 0._dp
354  atom%energy%exc = 0._dp
355  END IF
356 
357  ! Zero this contribution
358  atom%energy%elsd = 0._dp
359 
360  ! Total energy
361  atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
362 
363  ! Potential energy
364  atom%energy%epot = atom%energy%etot - atom%energy%ekin
365 
366  ! Total HF/KS matrix
367  fmat%op = hcore%op + jmat%op + kmat%op + xcmat%op
368 
369  ! calculate error matrix
370  CALL err_matrix(jmat%op, deps, fmat%op, atom%orbitals%pmat, atom%integrals%utrans, &
371  atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
372 
373  iter = iter + 1
374 
375  IF (iw > 0) THEN
376  IF (need_zmp) THEN
377  CALL atom_print_zmp_iteration(iter, deps, atom, iw)
378  ELSE
379  CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
380  END IF
381  END IF
382 
383  IF (deps < eps_scf) THEN
384  is_converged = .true.
385  EXIT
386  END IF
387  IF (iter >= max_iter) THEN
388  IF (iw > 0) THEN
389  WRITE (iw, "(A)") " No convergence within maximum number of iterations "
390  END IF
391  EXIT
392  END IF
393 
394  ! update history container and extrapolate KS matrix
395  CALL atom_history_update(history, atom%orbitals%pmat, fmat%op, jmat%op, atom%energy%etot, deps)
396  CALL atom_opt_fmat(fmat%op, history, deps)
397 
398  ! Solve HF/KS equations
399  CALL atom_solve(fmat%op, atom%integrals%utrans, atom%orbitals%wfn, atom%orbitals%ener, &
400  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
401  CALL atom_denmat(atom%orbitals%pmat, atom%orbitals%wfn, atom%basis%nbas, atom%state%occupation, &
402  atom%state%maxl_occ, atom%state%maxn_occ)
403 
404  END DO !SCF Loop
405 
406  IF (iw > 0) THEN
407  CALL atom_print_energies(atom, iw)
408  END IF
409 
410  CALL atom_history_release(history)
411 
412  ! conserve fmat within atom_type
413  CALL set_atom(atom, fmat=fmat)
414  CALL release_opmat(jmat)
415  CALL release_opmat(kmat)
416  CALL release_opmat(xcmat)
417  CALL release_opmat(hcore)
418 
419  CALL release_opgrid(density)
420  CALL release_opgrid(cpot)
421 
422  ! ZMP deallocating ext_density ext_vxc
423  IF (need_zmp) DEALLOCATE (ext_density)
424  IF (need_vxc) DEALLOCATE (ext_vxc)
425 
426  IF (PRESENT(converged)) THEN
427  converged = is_converged
428  END IF
429 
430  CALL timestop(handle)
431 
432  END SUBROUTINE calculate_atom_restricted
433 
434 ! **************************************************************************************************
435 !> \brief Perform unrestricted (spin-polarised) electronic structure atomic calculations.
436 !> \param atom information about the atomic kind
437 !> \param iw output file unit
438 !> \param noguess skip initial guess
439 !> \param converged whether SCF iterations have been converged
440 !> \par History
441 !> * identical to the subroutine calculate_atom_restricted()
442 ! **************************************************************************************************
443  SUBROUTINE calculate_atom_unrestricted(atom, iw, noguess, converged)
444  TYPE(atom_type), POINTER :: atom
445  INTEGER, INTENT(IN) :: iw
446  LOGICAL, INTENT(IN), OPTIONAL :: noguess
447  LOGICAL, INTENT(OUT), OPTIONAL :: converged
448 
449  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_atom_unrestricted'
450 
451  INTEGER :: handle, i, iter, k, l, max_iter, method, &
452  reltyp
453  LOGICAL :: do_hfx, doguess, is_converged, lsdpot, &
454  need_x, need_xc
455  REAL(kind=dp) :: deps, depsa, depsb, eps_scf, hf_frac, &
456  ne, nm
457  TYPE(atom_history_type) :: historya, historyb
458  TYPE(opgrid_type), POINTER :: cpot, density, rhoa, rhob
459  TYPE(opmat_type), POINTER :: fmata, fmatb, hcore, hlsd, jmat, kmata, &
460  kmatb, xcmata, xcmatb
461  TYPE(section_vals_type), POINTER :: xc_section
462 
463  CALL timeset(routinen, handle)
464 
465  IF (PRESENT(noguess)) THEN
466  doguess = .NOT. noguess
467  ELSE
468  doguess = .true.
469  END IF
470 
471  CALL setup_hf_section(hf_frac, do_hfx, atom, xc_section, atom%exchange_integral_type)
472 
473  method = atom%method_type
474  max_iter = atom%optimization%max_iter
475  eps_scf = atom%optimization%eps_scf
476 
477  SELECT CASE (method)
478  CASE DEFAULT
479  cpabort("")
480  CASE (do_rks_atom)
481  cpabort("")
482  CASE (do_uks_atom)
483  need_x = do_hfx
484  need_xc = .true.
485  CASE (do_rhf_atom)
486  cpabort("")
487  CASE (do_uhf_atom)
488  need_x = .true.
489  need_xc = .false.
490  hf_frac = 1._dp
491  CASE (do_rohf_atom)
492  need_x = .true.
493  need_xc = .false.
494  hf_frac = 1._dp
495  cpabort("")
496  END SELECT
497 
498  ! set alpha and beta occupations
499  IF (sum(abs(atom%state%occa) + abs(atom%state%occb)) == 0.0_dp) THEN
500  DO l = 0, 3
501  nm = real((2*l + 1), kind=dp)
502  DO k = 1, 10
503  ne = atom%state%occupation(l, k)
504  IF (ne == 0._dp) THEN !empty shell
505  EXIT !assume there are no holes
506  ELSEIF (ne == 2._dp*nm) THEN !closed shell
507  atom%state%occa(l, k) = nm
508  atom%state%occb(l, k) = nm
509  ELSEIF (atom%state%multiplicity == -2) THEN !High spin case
510  atom%state%occa(l, k) = min(ne, nm)
511  atom%state%occb(l, k) = max(0._dp, ne - nm)
512  ELSE
513  atom%state%occa(l, k) = 0.5_dp*(ne + atom%state%multiplicity - 1._dp)
514  atom%state%occb(l, k) = ne - atom%state%occa(l, k)
515  END IF
516  END DO
517  END DO
518  END IF
519  ! check for relativistic method
520  reltyp = atom%relativistic
521 
522  IF (iw > 0) CALL atom_print_state(atom%state, iw)
523 
524  NULLIFY (hcore, hlsd)
525  CALL create_opmat(hcore, atom%basis%nbas)
526  CALL create_opmat(hlsd, atom%basis%nbas)
527  hlsd%op = 0._dp
528  ! Pseudopotentials
529  lsdpot = .false.
530  SELECT CASE (atom%potential%ppot_type)
531  CASE DEFAULT
532  cpabort("")
533  CASE (no_pseudo)
534  SELECT CASE (reltyp)
535  CASE DEFAULT
536  cpabort("")
537  CASE (do_nonrel_atom)
538  hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core
539  CASE (do_zoramp_atom, do_sczoramp_atom)
540  hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core
541  CASE (do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom)
542  hcore%op = atom%integrals%hdkh
543  END SELECT
544  CASE (gth_pseudo)
545  hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
546  IF (atom%potential%gth_pot%lsdpot) THEN
547  lsdpot = .true.
548  hlsd%op = atom%integrals%clsd
549  END IF
550  CASE (upf_pseudo, sgp_pseudo, ecp_pseudo)
551  hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl
552  END SELECT
553  ! add confinement potential (not included in relativistic transformations)
554  IF (atom%potential%confinement) THEN
555  hcore%op = hcore%op + atom%potential%acon*atom%integrals%conf
556  END IF
557 
558  NULLIFY (fmata, fmatb, jmat, kmata, kmatb, xcmata, xcmatb)
559  CALL create_opmat(fmata, atom%basis%nbas)
560  CALL create_opmat(fmatb, atom%basis%nbas)
561  CALL create_opmat(jmat, atom%basis%nbas)
562  CALL create_opmat(kmata, atom%basis%nbas)
563  CALL create_opmat(kmatb, atom%basis%nbas)
564  CALL create_opmat(xcmata, atom%basis%nbas)
565  CALL create_opmat(xcmatb, atom%basis%nbas)
566 
567  NULLIFY (density, rhoa, rhob, cpot)
568  CALL create_opgrid(density, atom%basis%grid)
569  CALL create_opgrid(rhoa, atom%basis%grid)
570  CALL create_opgrid(rhob, atom%basis%grid)
571  CALL create_opgrid(cpot, atom%basis%grid)
572 
573  IF (doguess) THEN
574  ! initial guess
575  CALL slater_density(rhoa%op, rhob%op, atom%z, atom%state, atom%basis%grid)
576  density%op = rhoa%op + rhob%op
577  CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
578  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
579  ! alpha spin
580  density%op = 2._dp*rhoa%op
581  CALL wigner_slater_functional(density%op, cpot%op)
582  CALL numpot_matrix(xcmata%op, cpot%op, atom%basis, 0)
583  fmata%op = hcore%op + hlsd%op + jmat%op + xcmata%op
584  CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
585  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
586  ! beta spin
587  density%op = 2._dp*rhob%op
588  CALL wigner_slater_functional(density%op, cpot%op)
589  CALL numpot_matrix(xcmatb%op, cpot%op, atom%basis, 0)
590  fmatb%op = hcore%op - hlsd%op + jmat%op + xcmatb%op
591  CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
592  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
593  END IF
594  CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
595  atom%state%maxl_occ, atom%state%maxn_occ)
596  CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
597  atom%state%maxl_occ, atom%state%maxn_occ)
598  atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
599 
600  ! wavefunction history
601  NULLIFY (historya%dmat, historya%hmat)
602  CALL atom_history_init(historya, atom%optimization, fmata%op)
603  NULLIFY (historyb%dmat, historyb%hmat)
604  CALL atom_history_init(historyb, atom%optimization, fmatb%op)
605 
606  is_converged = .false.
607  iter = 0
608  DO !SCF Loop
609 
610  ! Kinetic energy
611  atom%energy%ekin = atom_trace(atom%integrals%kin, atom%orbitals%pmat)
612 
613  ! Band energy
614  atom%energy%eband = 0._dp
615  DO l = 0, 3
616  DO i = 1, min(SIZE(atom%state%occupation, 2), SIZE(atom%orbitals%ener, 1))
617  atom%energy%eband = atom%energy%eband + atom%orbitals%enera(i, l)*atom%state%occa(l, i)
618  atom%energy%eband = atom%energy%eband + atom%orbitals%enerb(i, l)*atom%state%occb(l, i)
619  END DO
620  END DO
621 
622  ! Pseudopotential energy
623  SELECT CASE (atom%potential%ppot_type)
624  CASE DEFAULT
625  cpabort("")
626  CASE (no_pseudo)
627  atom%energy%eploc = 0._dp
628  atom%energy%epnl = 0._dp
629  CASE (gth_pseudo, upf_pseudo, sgp_pseudo, ecp_pseudo)
630  atom%energy%eploc = atom_trace(atom%integrals%core, atom%orbitals%pmat)
631  atom%energy%epnl = atom_trace(atom%integrals%hnl, atom%orbitals%pmat)
632  END SELECT
633  atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl
634 
635  ! Core energy
636  atom%energy%ecore = atom_trace(hcore%op, atom%orbitals%pmat)
637 
638  ! Confinement energy
639  IF (atom%potential%confinement) THEN
640  atom%energy%econfinement = atom_trace(atom%integrals%conf, atom%orbitals%pmat)
641  ELSE
642  atom%energy%econfinement = 0._dp
643  END IF
644 
645  ! Hartree Term
646  jmat%op = 0._dp
647  SELECT CASE (atom%coulomb_integral_type)
648  CASE DEFAULT
649  cpabort("")
650  CASE (do_analytic)
651  CALL ceri_contract(jmat%op, atom%integrals%ceri, atom%orbitals%pmat, atom%integrals%n)
652  CASE (do_semi_analytic)
653  CALL coulomb_potential_analytic(cpot%op, atom%orbitals%pmat, atom%basis, atom%basis%grid, &
654  atom%state%maxl_occ)
655  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
656  CASE (do_numeric)
657  CALL atom_density(density%op, atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
658  CALL coulomb_potential_numeric(cpot%op, density%op, density%grid)
659  CALL numpot_matrix(jmat%op, cpot%op, atom%basis, 0)
660  END SELECT
661  atom%energy%ecoulomb = 0.5_dp*atom_trace(jmat%op, atom%orbitals%pmat)
662 
663  ! Exchange Term
664  IF (need_x) THEN
665  kmata%op = 0._dp
666  kmatb%op = 0._dp
667  SELECT CASE (atom%exchange_integral_type)
668  CASE DEFAULT
669  cpabort("")
670  CASE (do_analytic)
671  CALL eeri_contract(kmata%op, atom%integrals%eeri, atom%orbitals%pmata, atom%integrals%n)
672  CALL eeri_contract(kmatb%op, atom%integrals%eeri, atom%orbitals%pmatb, atom%integrals%n)
673  CASE (do_semi_analytic)
674  CALL exchange_semi_analytic(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
675  CALL exchange_semi_analytic(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
676  CASE (do_numeric)
677  CALL exchange_numeric(kmata%op, atom%state, atom%state%occa, atom%orbitals%wfna, atom%basis, atom%hfx_pot)
678  CALL exchange_numeric(kmatb%op, atom%state, atom%state%occb, atom%orbitals%wfnb, atom%basis, atom%hfx_pot)
679  END SELECT
680  atom%energy%eexchange = hf_frac*(atom_trace(kmata%op, atom%orbitals%pmata) + &
681  atom_trace(kmatb%op, atom%orbitals%pmatb))
682  kmata%op = 2._dp*hf_frac*kmata%op
683  kmatb%op = 2._dp*hf_frac*kmatb%op
684  ELSE
685  kmata%op = 0._dp
686  kmatb%op = 0._dp
687  atom%energy%eexchange = 0._dp
688  END IF
689 
690  ! XC
691  IF (need_xc) THEN
692  xcmata%op = 0._dp
693  xcmatb%op = 0._dp
694  CALL calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
695  ELSE
696  xcmata%op = 0._dp
697  xcmatb%op = 0._dp
698  atom%energy%exc = 0._dp
699  END IF
700 
701  IF (lsdpot) THEN
702  atom%energy%elsd = atom_trace(hlsd%op, atom%orbitals%pmata) - &
703  atom_trace(hlsd%op, atom%orbitals%pmatb)
704  atom%energy%epseudo = atom%energy%epseudo + atom%energy%elsd
705  atom%energy%ecore = atom%energy%ecore + atom%energy%elsd
706  ELSE
707  atom%energy%elsd = 0._dp
708  END IF
709 
710  ! Total energy
711  atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc
712 
713  ! Potential energy
714  atom%energy%epot = atom%energy%etot - atom%energy%ekin
715 
716  ! Total HF/KS matrix
717  fmata%op = hcore%op + hlsd%op + jmat%op + kmata%op + xcmata%op
718  fmatb%op = hcore%op - hlsd%op + jmat%op + kmatb%op + xcmatb%op
719 
720  ! calculate error matrix
721  CALL err_matrix(xcmata%op, depsa, fmata%op, atom%orbitals%pmata, atom%integrals%utrans, &
722  atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
723  CALL err_matrix(xcmatb%op, depsb, fmatb%op, atom%orbitals%pmatb, atom%integrals%utrans, &
724  atom%integrals%uptrans, atom%basis%nbas, atom%integrals%nne)
725  deps = 2._dp*max(depsa, depsb)
726 
727  iter = iter + 1
728 
729  IF (iw > 0) THEN
730  CALL atom_print_iteration(iter, deps, atom%energy%etot, iw)
731  END IF
732 
733  IF (deps < eps_scf) THEN
734  is_converged = .true.
735  EXIT
736  END IF
737  IF (iter >= max_iter) THEN
738  IF (iw > 0) THEN
739  WRITE (iw, "(A)") " No convergence within maximum number of iterations "
740  END IF
741  EXIT
742  END IF
743 
744  ! update history container and extrapolate KS matrix
745  CALL atom_history_update(historya, atom%orbitals%pmata, fmata%op, xcmata%op, atom%energy%etot, deps)
746  CALL atom_history_update(historyb, atom%orbitals%pmatb, fmatb%op, xcmatb%op, atom%energy%etot, deps)
747  CALL atom_opt_fmat(fmata%op, historya, depsa)
748  CALL atom_opt_fmat(fmatb%op, historyb, depsb)
749 
750  ! Solve HF/KS equations
751  CALL atom_solve(fmata%op, atom%integrals%utrans, atom%orbitals%wfna, atom%orbitals%enera, &
752  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
753  CALL atom_denmat(atom%orbitals%pmata, atom%orbitals%wfna, atom%basis%nbas, atom%state%occa, &
754  atom%state%maxl_occ, atom%state%maxn_occ)
755  CALL atom_solve(fmatb%op, atom%integrals%utrans, atom%orbitals%wfnb, atom%orbitals%enerb, &
756  atom%basis%nbas, atom%integrals%nne, atom%state%maxl_calc)
757  CALL atom_denmat(atom%orbitals%pmatb, atom%orbitals%wfnb, atom%basis%nbas, atom%state%occb, &
758  atom%state%maxl_occ, atom%state%maxn_occ)
759  atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb
760 
761  END DO !SCF Loop
762 
763  IF (iw > 0) THEN
764  CALL atom_print_energies(atom, iw)
765  END IF
766 
767  CALL atom_history_release(historya)
768  CALL atom_history_release(historyb)
769 
770  CALL release_opgrid(density)
771  CALL release_opgrid(rhoa)
772  CALL release_opgrid(rhob)
773  CALL release_opgrid(cpot)
774 
775  ! conserve fmata as fmat within atom_type.
776  CALL set_atom(atom, fmat=fmata)
777  CALL release_opmat(fmatb)
778  CALL release_opmat(jmat)
779  CALL release_opmat(kmata)
780  CALL release_opmat(kmatb)
781  CALL release_opmat(xcmata)
782  CALL release_opmat(xcmatb)
783  CALL release_opmat(hlsd)
784  CALL release_opmat(hcore)
785 
786  IF (PRESENT(converged)) THEN
787  converged = is_converged
788  END IF
789 
790  CALL timestop(handle)
791 
792  END SUBROUTINE calculate_atom_unrestricted
793 
794 END MODULE atom_electronic_structure
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_optimization
Optimizer for the atomic code.
Definition: atom_optimization.F:11
atom_optimization::atom_opt_fmat
subroutine, public atom_opt_fmat(fmat, history, err)
Construct a Kohn-Sham matrix for the next iteration based on the historic data.
Definition: atom_optimization.F:172
atom_optimization::atom_history_init
pure subroutine, public atom_history_init(history, optimization, matrix)
Initialise a circular buffer to keep Kohn-Sham and density matrices from previous iteration.
Definition: atom_optimization.F:57
atom_optimization::atom_history_release
pure subroutine, public atom_history_release(history)
Release circular buffer to keep historic matrices.
Definition: atom_optimization.F:133
atom_optimization::atom_history_update
pure subroutine, public atom_history_update(history, pmat, fmat, emat, energy, error)
Add matrices from the current iteration into the circular buffer.
Definition: atom_optimization.F:104
atom_output
Routines that print various information about an atomic kind.
Definition: atom_output.F:11
atom_output::atom_print_energies
subroutine, public atom_print_energies(atom, iw)
Print energy components.
Definition: atom_output.F:175
atom_output::atom_print_iteration
subroutine, public atom_print_iteration(iter, deps, etot, iw)
Print convergence information.
Definition: atom_output.F:272
atom_output::atom_print_state
subroutine, public atom_print_state(state, iw)
Print information about electronic state.
Definition: atom_output.F:83
atom_output::atom_print_zmp_iteration
subroutine, public atom_print_zmp_iteration(iter, deps, atom, iw)
Printing of the atomic iterations when ZMP is active.
Definition: atom_output.F:248
atom_types
Define the atom type and its sub types.
Definition: atom_types.F:15
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::setup_hf_section
subroutine, public setup_hf_section(hf_frac, do_hfx, atom, xc_section, extype)
...
Definition: atom_types.F:1170
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::slater_density
subroutine, public slater_density(density1, density2, zcore, state, grid)
Calculate Slater density on a radial grid.
Definition: atom_utils.F:2094
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_read_external_vxc
subroutine, public atom_read_external_vxc(vxc, atom, iw)
ZMP subroutine to read external v_xc in the atomic code.
Definition: atom_utils.F:608
atom_utils::eeri_contract
pure subroutine, public eeri_contract(kmat, erint, pmat, nsize)
Contract exchange Electron Repulsion Integrals.
Definition: atom_utils.F:1974
atom_utils::exchange_semi_analytic
subroutine, public exchange_semi_analytic(kmat, state, occ, wfn, basis, hfx_pot)
Calculate the exchange potential semi-analytically.
Definition: atom_utils.F:1308
atom_utils::ceri_contract
pure subroutine, public ceri_contract(jmat, erint, pmat, nsize, all_nu)
Contract Coulomb Electron Repulsion Integrals.
Definition: atom_utils.F:1900
atom_utils::coulomb_potential_analytic
subroutine, public coulomb_potential_analytic(cpot, pmat, basis, grid, maxl)
Analytically compute the Coulomb potential on an atomic radial grid.
Definition: atom_utils.F:1115
atom_utils::err_matrix
pure subroutine, public err_matrix(emat, demax, kmat, pmat, umat, upmat, nval, nbs)
Calculate the error matrix for each angular momentum.
Definition: atom_utils.F:2051
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::wigner_slater_functional
pure subroutine, public wigner_slater_functional(rho, vxc)
Calculate the functional derivative of the Wigner (correlation) - Slater (exchange) density functiona...
Definition: atom_utils.F:2159
atom_utils::atom_read_zmp_restart
subroutine, public atom_read_zmp_restart(atom, doguess, iw)
ZMP subroutine to read external restart file.
Definition: atom_utils.F:458
atom_utils::exchange_numeric
subroutine, public exchange_numeric(kmat, state, occ, wfn, basis, hfx_pot)
Calculate the exchange potential numerically.
Definition: atom_utils.F:1220
atom_utils::atom_trace
pure real(kind=dp) function, public atom_trace(opmat, pmat)
Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
Definition: atom_utils.F:1823
atom_utils::atom_read_external_density
subroutine, public atom_read_external_density(density, atom, iw)
ZMP subroutine to read external density from linear grid of density matrix.
Definition: atom_utils.F:517
atom_utils::numpot_matrix
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
Definition: atom_utils.F:1840
atom_utils::atom_solve
subroutine, public atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
Solve the generalised eigenproblem for every angular momentum.
Definition: atom_utils.F:944
atom_utils::coulomb_potential_numeric
subroutine, public coulomb_potential_numeric(cpot, density, grid)
Numerically compute the Coulomb potential on an atomic radial grid.
Definition: atom_utils.F:1077
atom_xc
routines that build the integrals of the Vxc potential calculated for the atomic code
Definition: atom_xc.F:12
atom_xc::calculate_atom_vxc_lda
subroutine, public calculate_atom_vxc_lda(xcmat, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-restricted case.
Definition: atom_xc.F:199
atom_xc::calculate_atom_ext_vxc
subroutine, public calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
ZMP subroutine for the scf external density from the external v_xc.
Definition: atom_xc.F:148
atom_xc::calculate_atom_zmp
subroutine, public calculate_atom_zmp(ext_density, atom, lprint, xcmat)
ZMP subroutine for the calculation of the effective constraint potential in the atomic code.
Definition: atom_xc.F:78
atom_xc::calculate_atom_vxc_lsd
subroutine, public calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-polarized case.
Definition: atom_xc.F:404
atom
Definition: atom.F:9
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_dkh3_atom
integer, parameter, public do_dkh3_atom
Definition: input_constants.F:165
input_constants::do_nonrel_atom
integer, parameter, public do_nonrel_atom
Definition: input_constants.F:165
input_constants::do_dkh0_atom
integer, parameter, public do_dkh0_atom
Definition: input_constants.F:165
input_constants::do_uhf_atom
integer, parameter, public do_uhf_atom
Definition: input_constants.F:159
input_constants::do_dkh2_atom
integer, parameter, public do_dkh2_atom
Definition: input_constants.F:165
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition: input_constants.F:159
input_constants::do_numeric
integer, parameter, public do_numeric
Definition: input_constants.F:149
input_constants::do_zoramp_atom
integer, parameter, public do_zoramp_atom
Definition: input_constants.F:165
input_constants::do_dkh1_atom
integer, parameter, public do_dkh1_atom
Definition: input_constants.F:165
input_constants::do_rohf_atom
integer, parameter, public do_rohf_atom
Definition: input_constants.F:159
input_constants::do_semi_analytic
integer, parameter, public do_semi_analytic
Definition: input_constants.F:149
input_constants::do_sczoramp_atom
integer, parameter, public do_sczoramp_atom
Definition: input_constants.F:165
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_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34