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