(git:374b731)
Loading...
Searching...
No Matches
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
8MODULE 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
47CONTAINS
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
794END 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::sgp_pseudo
integer, parameter, public sgp_pseudo
Definition input_constants.F:174
input_constants::do_dkh3_atom
integer, parameter, public do_dkh3_atom
Definition input_constants.F:165
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
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::upf_pseudo
integer, parameter, public upf_pseudo
Definition input_constants.F:174
input_constants::do_dkh2_atom
integer, parameter, public do_dkh2_atom
Definition input_constants.F:165
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
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
atom_optimization::atom_history_type
Definition atom_optimization.F:31
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:290
atom_types::opgrid_type
Operator grids.
Definition atom_types.F:255
atom_types::opmat_type
Operator matrices.
Definition atom_types.F:248
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126