(git:e5fdd81)
atom_utils.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Some basic routines for atomic calculations
10 !> \author jgh
11 !> \date 01.04.2008
12 !> \version 1.0
13 !>
14 ! **************************************************************************************************
15 MODULE atom_utils
16  USE ai_onecenter, ONLY: sg_overlap,&
17  sto_overlap
18  USE ai_overlap, ONLY: overlap_ab_s,&
19  overlap_ab_sp
20  USE ao_util, ONLY: exp_radius
21  USE atom_types, ONLY: &
22  cgto_basis, gto_basis, num_basis, sto_basis, atom_basis_type, atom_gthpot_type, &
23  atom_hfx_type, atom_potential_type, atom_state, atom_type, ecp_pseudo, eri, gth_pseudo, &
24  lmat, no_pseudo, sgp_pseudo, upf_pseudo
25  USE basis_set_types, ONLY: srules
26  USE cp_files, ONLY: close_file,&
27  get_unit_number,&
28  open_file
29  USE input_constants, ONLY: do_rhf_atom,&
30  do_rks_atom,&
31  do_rohf_atom,&
32  do_uhf_atom,&
33  do_uks_atom
34  USE kahan_sum, ONLY: accurate_dot_product
35  USE kinds, ONLY: default_string_length,&
36  dp
37  USE lapack, ONLY: lapack_ssyev
38  USE mathconstants, ONLY: dfac,&
39  fac,&
40  fourpi,&
41  maxfac,&
42  pi,&
43  rootpi
44  USE mathlib, ONLY: binomial_gen,&
45  invmat_symm
46  USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
47  init_orbital_pointers
48  USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
49  init_spherical_harmonics
50  USE periodic_table, ONLY: nelem,&
51  ptable
52  USE physcon, ONLY: bohr
53  USE qs_grid_atom, ONLY: grid_atom_type
54  USE splines, ONLY: spline3ders
55  USE string_utilities, ONLY: uppercase
56 #include "./base/base_uses.f90"
57 
58  IMPLICIT NONE
59 
60  PRIVATE
61 
62  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_utils'
63 
64  PUBLIC :: atom_condnumber, atom_completeness, atom_basis_condnum
65  PUBLIC :: atom_set_occupation, get_maxl_occ, get_maxn_occ
66  PUBLIC :: atom_denmat, atom_density, atom_core_density
67  PUBLIC :: integrate_grid, atom_trace, atom_solve
68  PUBLIC :: coulomb_potential_numeric, coulomb_potential_analytic
69  PUBLIC :: exchange_numeric, exchange_semi_analytic
70  PUBLIC :: numpot_matrix, ceri_contract, eeri_contract, err_matrix
71  PUBLIC :: slater_density, wigner_slater_functional, atom_local_potential
72  PUBLIC :: atom_orbital_charge, atom_orbital_nodes, atom_consistent_method
73  PUBLIC :: atom_orbital_max, atom_wfnr0, get_rho0
74  PUBLIC :: contract2, contract4, contract2add
75 ! ZMP added public subroutines
76  PUBLIC :: atom_read_external_density
77  PUBLIC :: atom_read_external_vxc
78  PUBLIC :: atom_read_zmp_restart
79  PUBLIC :: atom_write_zmp_restart
80 
81 !-----------------------------------------------------------------------------!
82 
83  INTERFACE integrate_grid
84  MODULE PROCEDURE integrate_grid_function1, &
85  integrate_grid_function2, &
86  integrate_grid_function3
87  END INTERFACE
88 
89 ! **************************************************************************************************
90 
91 CONTAINS
92 
93 ! **************************************************************************************************
94 !> \brief Set occupation of atomic orbitals.
95 !> \param ostring list of electronic states
96 !> \param occupation ...
97 !> \param wfnocc ...
98 !> \param multiplicity ...
99 !> \par History
100 !> * 11.2009 unrestricted KS and HF methods [Juerg Hutter]
101 !> * 08.2008 created [Juerg Hutter]
102 ! **************************************************************************************************
103  SUBROUTINE atom_set_occupation(ostring, occupation, wfnocc, multiplicity)
104  CHARACTER(LEN=default_string_length), &
105  DIMENSION(:), POINTER :: ostring
106  REAL(kind=dp), DIMENSION(0:lmat, 10) :: occupation, wfnocc
107  INTEGER, INTENT(OUT), OPTIONAL :: multiplicity
108 
109  CHARACTER(len=2) :: elem
110  CHARACTER(LEN=default_string_length) :: pstring
111  INTEGER :: i, i1, i2, ielem, is, jd, jf, jp, js, k, &
112  l, mult, n, no
113  REAL(kind=dp) :: e0, el, oo
114 
115  occupation = 0._dp
116 
117  cpassert(ASSOCIATED(ostring))
118  cpassert(SIZE(ostring) > 0)
119 
120  no = SIZE(ostring)
121 
122  is = 1
123  ! look for multiplicity
124  mult = -1 !not specified
125  IF (is <= no) THEN
126  IF (index(ostring(is), "(") /= 0) THEN
127  i1 = index(ostring(is), "(")
128  i2 = index(ostring(is), ")")
129  cpassert((i2 - i1 - 1 > 0) .AND. (i2 - i1 - 1 < 3))
130  elem = ostring(is) (i1 + 1:i2 - 1)
131  IF (index(elem, "HS") /= 0) THEN
132  mult = -2 !High spin
133  ELSE IF (index(elem, "LS") /= 0) THEN
134  mult = -3 !Low spin
135  ELSE
136  READ (elem, *) mult
137  END IF
138  is = is + 1
139  END IF
140  END IF
141 
142  IF (is <= no) THEN
143  IF (index(ostring(is), "CORE") /= 0) is = is + 1 !Pseudopotential detected
144  END IF
145  IF (is <= no) THEN
146  IF (index(ostring(is), "none") /= 0) is = is + 1 !no electrons, used with CORE
147  END IF
148  IF (is <= no) THEN
149  IF (index(ostring(is), "[") /= 0) THEN
150  ! core occupation from element [XX]
151  i1 = index(ostring(is), "[")
152  i2 = index(ostring(is), "]")
153  cpassert((i2 - i1 - 1 > 0) .AND. (i2 - i1 - 1 < 3))
154  elem = ostring(is) (i1 + 1:i2 - 1)
155  ielem = 0
156  DO k = 1, nelem
157  IF (elem == ptable(k)%symbol) THEN
158  ielem = k
159  EXIT
160  END IF
161  END DO
162  cpassert(ielem /= 0)
163  DO l = 0, min(lmat, ubound(ptable(ielem)%e_conv, 1))
164  el = 2._dp*(2._dp*real(l, dp) + 1._dp)
165  e0 = ptable(ielem)%e_conv(l)
166  DO k = 1, 10
167  occupation(l, k) = min(el, e0)
168  e0 = e0 - el
169  IF (e0 <= 0._dp) EXIT
170  END DO
171  END DO
172  is = is + 1
173  END IF
174 
175  END IF
176 
177  DO i = is, no
178  pstring = ostring(i)
179  CALL uppercase(pstring)
180  js = index(pstring, "S")
181  jp = index(pstring, "P")
182  jd = index(pstring, "D")
183  jf = index(pstring, "F")
184  cpassert(js + jp + jd + jf > 0)
185  IF (js > 0) THEN
186  cpassert(jp + jd + jf == 0)
187  READ (pstring(1:js - 1), *) n
188  READ (pstring(js + 1:), *) oo
189  cpassert(n > 0)
190  cpassert(oo >= 0._dp)
191  cpassert(occupation(0, n) == 0)
192  occupation(0, n) = oo
193  END IF
194  IF (jp > 0) THEN
195  cpassert(js + jd + jf == 0)
196  READ (pstring(1:jp - 1), *) n
197  READ (pstring(jp + 1:), *) oo
198  cpassert(n > 1)
199  cpassert(oo >= 0._dp)
200  cpassert(occupation(1, n - 1) == 0)
201  occupation(1, n - 1) = oo
202  END IF
203  IF (jd > 0) THEN
204  cpassert(js + jp + jf == 0)
205  READ (pstring(1:jd - 1), *) n
206  READ (pstring(jd + 1:), *) oo
207  cpassert(n > 2)
208  cpassert(oo >= 0._dp)
209  cpassert(occupation(2, n - 2) == 0)
210  occupation(2, n - 2) = oo
211  END IF
212  IF (jf > 0) THEN
213  cpassert(js + jp + jd == 0)
214  READ (pstring(1:jf - 1), *) n
215  READ (pstring(jf + 1:), *) oo
216  cpassert(n > 3)
217  cpassert(oo >= 0._dp)
218  cpassert(occupation(3, n - 3) == 0)
219  occupation(3, n - 3) = oo
220  END IF
221 
222  END DO
223 
224  wfnocc = 0._dp
225  DO l = 0, lmat
226  k = 0
227  DO i = 1, 10
228  IF (occupation(l, i) /= 0._dp) THEN
229  k = k + 1
230  wfnocc(l, k) = occupation(l, i)
231  END IF
232  END DO
233  END DO
234 
235  !Check for consistency with multiplicity
236  IF (mult /= -1) THEN
237  ! count open shells
238  js = 0
239  DO l = 0, lmat
240  k = 2*(2*l + 1)
241  DO i = 1, 10
242  IF (wfnocc(l, i) /= 0._dp .AND. wfnocc(l, i) /= real(k, dp)) THEN
243  js = js + 1
244  i1 = l
245  i2 = i
246  END IF
247  END DO
248  END DO
249 
250  IF (js == 0 .AND. mult == -2) mult = 1
251  IF (js == 0 .AND. mult == -3) mult = 1
252  IF (js == 0) THEN
253  cpassert(mult == 1)
254  END IF
255  IF (js == 1) THEN
256  l = i1
257  i = i2
258  k = nint(wfnocc(l, i))
259  IF (k > (2*l + 1)) k = 2*(2*l + 1) - k
260  IF (mult == -2) mult = k + 1
261  IF (mult == -3) mult = mod(k, 2) + 1
262  cpassert(mod(k + 1 - mult, 2) == 0)
263  END IF
264  IF (js > 1 .AND. mult /= -2) THEN
265  cpassert(mult == -2)
266  END IF
267 
268  END IF
269 
270  IF (PRESENT(multiplicity)) multiplicity = mult
271 
272  END SUBROUTINE atom_set_occupation
273 
274 ! **************************************************************************************************
275 !> \brief Return the maximum orbital quantum number of occupied orbitals.
276 !> \param occupation ...
277 !> \return ...
278 !> \par History
279 !> * 08.2008 created [Juerg Hutter]
280 ! **************************************************************************************************
281  PURE FUNCTION get_maxl_occ(occupation) RESULT(maxl)
282  REAL(kind=dp), DIMENSION(0:lmat, 10), INTENT(IN) :: occupation
283  INTEGER :: maxl
284 
285  INTEGER :: l
286 
287  maxl = 0
288  DO l = 0, lmat
289  IF (sum(occupation(l, :)) /= 0._dp) maxl = l
290  END DO
291 
292  END FUNCTION get_maxl_occ
293 
294 ! **************************************************************************************************
295 !> \brief Return the maximum principal quantum number of occupied orbitals.
296 !> \param occupation ...
297 !> \return ...
298 !> \par History
299 !> * 08.2008 created [Juerg Hutter]
300 ! **************************************************************************************************
301  PURE FUNCTION get_maxn_occ(occupation) RESULT(maxn)
302  REAL(kind=dp), DIMENSION(0:lmat, 10), INTENT(IN) :: occupation
303  INTEGER, DIMENSION(0:lmat) :: maxn
304 
305  INTEGER :: k, l
306 
307  maxn = 0
308  DO l = 0, lmat
309  DO k = 1, 10
310  IF (occupation(l, k) /= 0._dp) maxn(l) = maxn(l) + 1
311  END DO
312  END DO
313 
314  END FUNCTION get_maxn_occ
315 
316 ! **************************************************************************************************
317 !> \brief Calculate a density matrix using atomic orbitals.
318 !> \param pmat electron density matrix
319 !> \param wfn atomic wavefunctions
320 !> \param nbas number of basis functions
321 !> \param occ occupation numbers
322 !> \param maxl maximum angular momentum to consider
323 !> \param maxn maximum principal quantum number for each angular momentum
324 !> \par History
325 !> * 08.2008 created [Juerg Hutter]
326 ! **************************************************************************************************
327  PURE SUBROUTINE atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
328  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: pmat
329  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: wfn
330  INTEGER, DIMENSION(0:lmat), INTENT(IN) :: nbas
331  REAL(kind=dp), DIMENSION(0:, :), INTENT(IN) :: occ
332  INTEGER, INTENT(IN) :: maxl
333  INTEGER, DIMENSION(0:lmat), INTENT(IN) :: maxn
334 
335  INTEGER :: i, j, k, l, n
336 
337  pmat = 0._dp
338  n = SIZE(wfn, 2)
339  DO l = 0, maxl
340  DO i = 1, min(n, maxn(l))
341  DO k = 1, nbas(l)
342  DO j = 1, nbas(l)
343  pmat(j, k, l) = pmat(j, k, l) + occ(l, i)*wfn(j, i, l)*wfn(k, i, l)
344  END DO
345  END DO
346  END DO
347  END DO
348 
349  END SUBROUTINE atom_denmat
350 
351 ! **************************************************************************************************
352 !> \brief Map the electron density on an atomic radial grid.
353 !> \param density computed electron density
354 !> \param pmat electron density matrix
355 !> \param basis atomic basis set
356 !> \param maxl maximum angular momentum to consider
357 !> \param typ type of the matrix to map:
358 !> RHO -- density matrix;
359 !> DER -- first derivatives of the electron density;
360 !> KIN -- kinetic energy density;
361 !> LAP -- second derivatives of the electron density.
362 !> \param rr abscissa on the radial grid (required for typ == 'KIN')
363 !> \par History
364 !> * 08.2008 created [Juerg Hutter]
365 ! **************************************************************************************************
366  SUBROUTINE atom_density(density, pmat, basis, maxl, typ, rr)
367  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: density
368  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
369  TYPE(atom_basis_type), INTENT(IN) :: basis
370  INTEGER, INTENT(IN) :: maxl
371  CHARACTER(LEN=*), OPTIONAL :: typ
372  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: rr
373 
374  CHARACTER(LEN=3) :: my_typ
375  INTEGER :: i, j, l, n
376  REAL(kind=dp) :: ff
377 
378  my_typ = "RHO"
379  IF (PRESENT(typ)) my_typ = typ(1:3)
380  IF (my_typ == "KIN") THEN
381  cpassert(PRESENT(rr))
382  END IF
383 
384  density = 0._dp
385  DO l = 0, maxl
386  n = basis%nbas(l)
387  DO i = 1, n
388  DO j = i, n
389  ff = pmat(i, j, l)
390  IF (i /= j) ff = 2._dp*pmat(i, j, l)
391  IF (my_typ == "RHO") THEN
392  density(:) = density(:) + ff*basis%bf(:, i, l)*basis%bf(:, j, l)
393  ELSE IF (my_typ == "DER") THEN
394  density(:) = density(:) + ff*basis%dbf(:, i, l)*basis%bf(:, j, l) &
395  + ff*basis%bf(:, i, l)*basis%dbf(:, j, l)
396  ELSE IF (my_typ == "KIN") THEN
397  density(:) = density(:) + 0.5_dp*ff*( &
398  basis%dbf(:, i, l)*basis%dbf(:, j, l) + &
399  REAL(l*(l + 1), dp)*basis%bf(:, i, l)*basis%bf(:, j, l)/rr(:))
400  ELSE IF (my_typ == "LAP") THEN
401  density(:) = density(:) + ff*basis%ddbf(:, i, l)*basis%bf(:, j, l) &
402  + ff*basis%bf(:, i, l)*basis%ddbf(:, j, l) &
403  + 2._dp*ff*basis%dbf(:, i, l)*basis%bf(:, j, l)/rr(:) &
404  + 2._dp*ff*basis%bf(:, i, l)*basis%dbf(:, j, l)/rr(:)
405  ELSE
406  cpabort("")
407  END IF
408  END DO
409  END DO
410  END DO
411  ! this factor from the product of two spherical harmonics
412  density = density/fourpi
413 
414  END SUBROUTINE atom_density
415 
416 ! **************************************************************************************************
417 !> \brief ZMP subroutine to write external restart file.
418 !> \param atom information about the atomic kind
419 !> \date 07.10.2013
420 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
421 !> \version 1.0
422 ! **************************************************************************************************
423  SUBROUTINE atom_write_zmp_restart(atom)
424  TYPE(atom_type), INTENT(IN) :: atom
425 
426  INTEGER :: extunit, i, k, l, n
427 
428  extunit = get_unit_number()
429  CALL open_file(file_name=atom%zmp_restart_file, file_status="UNKNOWN", &
430  file_form="FORMATTED", file_action="WRITE", &
431  unit_number=extunit)
432 
433  n = SIZE(atom%orbitals%wfn, 2)
434  WRITE (extunit, *) atom%basis%nbas
435  DO l = 0, atom%state%maxl_occ
436  DO i = 1, min(n, atom%state%maxn_occ(l))
437  DO k = 1, atom%basis%nbas(l)
438  WRITE (extunit, *) atom%orbitals%wfn(k, i, l)
439  END DO
440  END DO
441  END DO
442 
443  CALL close_file(unit_number=extunit)
444 
445  END SUBROUTINE atom_write_zmp_restart
446 
447 ! **************************************************************************************************
448 !> \brief ZMP subroutine to read external restart file.
449 !> \param atom information about the atomic kind
450 !> \param doguess flag that indicates that the restart file has not been read,
451 !> so the initial guess is required
452 !> \param iw output file unit
453 !> \date 07.10.2013
454 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
455 !> \version 1.0
456 ! **************************************************************************************************
457  SUBROUTINE atom_read_zmp_restart(atom, doguess, iw)
458  TYPE(atom_type), INTENT(INOUT) :: atom
459  LOGICAL, INTENT(INOUT) :: doguess
460  INTEGER, INTENT(IN) :: iw
461 
462  INTEGER :: er, extunit, i, k, l, maxl, n
463  INTEGER, DIMENSION(0:lmat) :: maxn, nbas
464 
465  INQUIRE (file=atom%zmp_restart_file, exist=atom%doread)
466 
467  IF (atom%doread) THEN
468  WRITE (iw, fmt="(' ZMP | Restart file found')")
469  extunit = get_unit_number()
470  CALL open_file(file_name=atom%zmp_restart_file, file_status="OLD", &
471  file_form="FORMATTED", file_action="READ", &
472  unit_number=extunit)
473 
474  READ (extunit, *, iostat=er) nbas
475 
476  IF (er .NE. 0) THEN
477  WRITE (iw, fmt="(' ZMP | ERROR! Restart file unreadable')")
478  WRITE (iw, fmt="(' ZMP | ERROR! Starting ZMP calculation form initial atomic guess')")
479  doguess = .true.
480  atom%doread = .false.
481  ELSE IF (nbas(1) .NE. atom%basis%nbas(1)) THEN
482  WRITE (iw, fmt="(' ZMP | ERROR! Restart file contains a different basis set')")
483  WRITE (iw, fmt="(' ZMP | ERROR! Starting ZMP calculation form initial atomic guess')")
484  doguess = .true.
485  atom%doread = .false.
486  ELSE
487  nbas = atom%basis%nbas
488  maxl = atom%state%maxl_occ
489  maxn = atom%state%maxn_occ
490  n = SIZE(atom%orbitals%wfn, 2)
491  DO l = 0, atom%state%maxl_occ
492  DO i = 1, min(n, atom%state%maxn_occ(l))
493  DO k = 1, atom%basis%nbas(l)
494  READ (extunit, *) atom%orbitals%wfn(k, i, l)
495  END DO
496  END DO
497  END DO
498  doguess = .false.
499  END IF
500  CALL close_file(unit_number=extunit)
501  ELSE
502  WRITE (iw, fmt="(' ZMP | WARNING! Restart file not found')")
503  WRITE (iw, fmt="(' ZMP | WARNING! Starting ZMP calculation form initial atomic guess')")
504  END IF
505  END SUBROUTINE atom_read_zmp_restart
506 
507 ! **************************************************************************************************
508 !> \brief ZMP subroutine to read external density from linear grid of density matrix.
509 !> \param density external density
510 !> \param atom information about the atomic kind
511 !> \param iw output file unit
512 !> \date 07.10.2013
513 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
514 !> \version 1.0
515 ! **************************************************************************************************
516  SUBROUTINE atom_read_external_density(density, atom, iw)
517  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: density
518  TYPE(atom_type), INTENT(INOUT) :: atom
519  INTEGER, INTENT(IN) :: iw
520 
521  CHARACTER(LEN=default_string_length) :: filename
522  INTEGER :: extunit, ir, j, k, l, maxl_occ, maxnbas, &
523  nbas, nr
524  LOGICAL :: ldm
525  REAL(kind=dp) :: rr
526  REAL(kind=dp), ALLOCATABLE :: pmatread(:, :, :)
527 
528  filename = atom%ext_file
529  ldm = atom%dm
530  extunit = get_unit_number()
531 
532  CALL open_file(file_name=filename, file_status="OLD", &
533  file_form="FORMATTED", file_action="READ", &
534  unit_number=extunit)
535 
536  IF (.NOT. ldm) THEN
537  READ (extunit, *) nr
538 
539  IF (nr .NE. atom%basis%grid%nr) THEN
540  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! External grid dimension ',I4,' differs from internal grid ',I4)") &
541  nr, atom%basis%grid%nr
542  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Stopping ZMP calculation')")
543  cpabort("")
544  END IF
545 
546  DO ir = 1, nr
547  READ (extunit, *) rr, density(ir)
548  IF (abs(rr - atom%basis%grid%rad(ir)) .GT. atom%zmpgrid_tol) THEN
549  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Grid points do not coincide: ')")
550  IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T20,"R_out[bohr]",T36,"R_in[bohr]",T61,"R_diff[bohr]")')
551  IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T14,E24.15,T39,E24.15,T64,E24.15)') &
552  rr, atom%basis%grid%rad(ir), abs(rr - atom%basis%grid%rad(ir))
553  cpabort("")
554  END IF
555  END DO
556  CALL close_file(unit_number=extunit)
557  ELSE
558  READ (extunit, *) maxl_occ
559  maxnbas = maxval(atom%basis%nbas)
560  ALLOCATE (pmatread(maxnbas, maxnbas, 0:maxl_occ))
561  pmatread = 0.0
562  DO l = 0, maxl_occ
563  nbas = atom%basis%nbas(l)
564  READ (extunit, *) ! Read empty line
565  DO k = 1, nbas
566  READ (extunit, *) (pmatread(j, k, l), j=1, k)
567  DO j = 1, k
568  pmatread(k, j, l) = pmatread(j, k, l)
569  END DO
570  END DO
571  END DO
572 
573  CALL close_file(unit_number=extunit)
574 
575  CALL atom_density(density, pmatread, atom%basis, maxl_occ, typ="RHO")
576 
577  extunit = get_unit_number()
578  CALL open_file(file_name="rho_target.dat", file_status="UNKNOWN", &
579  file_form="FORMATTED", file_action="WRITE", unit_number=extunit)
580 
581  IF (iw > 0) WRITE (iw, fmt="(' ZMP | Writing target density from density matrix')")
582 
583  WRITE (extunit, fmt='("# Target density built from density matrix : ",A20)') filename
584  WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]")')
585 
586  nr = atom%basis%grid%nr
587 
588  DO ir = 1, nr
589  WRITE (extunit, fmt='(T1,E24.15,T26,E24.15)') &
590  atom%basis%grid%rad(ir), density(ir)
591  END DO
592  DEALLOCATE (pmatread)
593 
594  CALL close_file(unit_number=extunit)
595 
596  END IF
597 
598  END SUBROUTINE atom_read_external_density
599 
600 ! **************************************************************************************************
601 !> \brief ZMP subroutine to read external v_xc in the atomic code.
602 !> \param vxc external exchange-correlation potential
603 !> \param atom information about the atomic kind
604 !> \param iw output file unit
605 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
606 ! **************************************************************************************************
607  SUBROUTINE atom_read_external_vxc(vxc, atom, iw)
608  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: vxc
609  TYPE(atom_type), INTENT(INOUT) :: atom
610  INTEGER, INTENT(IN) :: iw
611 
612  CHARACTER(LEN=default_string_length) :: adum, filename
613  INTEGER :: extunit, ir, nr
614  REAL(kind=dp) :: rr
615 
616  filename = atom%ext_vxc_file
617  extunit = get_unit_number()
618 
619  CALL open_file(file_name=filename, file_status="OLD", &
620  file_form="FORMATTED", file_action="READ", &
621  unit_number=extunit)
622 
623  READ (extunit, *)
624  READ (extunit, *) adum, nr
625 
626  IF (nr .NE. atom%basis%grid%nr) THEN
627  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! External grid dimension ',I4,' differs from internal grid ',I4)") &
628  nr, atom%basis%grid%nr
629  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Stopping ZMP calculation')")
630  cpabort("")
631  END IF
632  DO ir = 1, nr
633  READ (extunit, *) rr, vxc(ir)
634  IF (abs(rr - atom%basis%grid%rad(ir)) .GT. atom%zmpvxcgrid_tol) THEN
635  IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Grid points do not coincide: ')")
636  IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T20,"R_out[bohr]",T36,"R_in[bohr]",T61,"R_diff[bohr]")')
637  IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T14,E24.15,T39,E24.15,T64,E24.15)') &
638  rr, atom%basis%grid%rad(ir), abs(rr - atom%basis%grid%rad(ir))
639  cpabort("")
640  END IF
641  END DO
642 
643  END SUBROUTINE atom_read_external_vxc
644 
645 ! **************************************************************************************************
646 !> \brief ...
647 !> \param charge ...
648 !> \param wfn ...
649 !> \param rcov ...
650 !> \param l ...
651 !> \param basis ...
652 ! **************************************************************************************************
653  PURE SUBROUTINE atom_orbital_charge(charge, wfn, rcov, l, basis)
654  REAL(kind=dp), INTENT(OUT) :: charge
655  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wfn
656  REAL(kind=dp), INTENT(IN) :: rcov
657  INTEGER, INTENT(IN) :: l
658  TYPE(atom_basis_type), INTENT(IN) :: basis
659 
660  INTEGER :: i, j, m, n
661  REAL(kind=dp) :: ff
662  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: den
663 
664  charge = 0._dp
665  m = SIZE(basis%bf, 1)
666  ALLOCATE (den(m))
667  n = basis%nbas(l)
668  den = 0._dp
669  DO i = 1, n
670  DO j = 1, n
671  ff = wfn(i)*wfn(j)
672  den(1:m) = den(1:m) + ff*basis%bf(1:m, i, l)*basis%bf(1:m, j, l)
673  END DO
674  END DO
675  DO i = 1, m
676  IF (basis%grid%rad(i) > rcov) den(i) = 0._dp
677  END DO
678  charge = sum(den(1:m)*basis%grid%wr(1:m))
679  DEALLOCATE (den)
680 
681  END SUBROUTINE atom_orbital_charge
682 
683 ! **************************************************************************************************
684 !> \brief ...
685 !> \param corden ...
686 !> \param potential ...
687 !> \param typ ...
688 !> \param rr ...
689 !> \par History
690 !> * 01.2017 rewritten [Juerg Hutter]
691 !> * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
692 ! **************************************************************************************************
693  SUBROUTINE atom_core_density(corden, potential, typ, rr)
694  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: corden
695  TYPE(atom_potential_type), INTENT(IN) :: potential
696  CHARACTER(LEN=*), OPTIONAL :: typ
697  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rr
698 
699  CHARACTER(LEN=3) :: my_typ
700  INTEGER :: i, j, m, n
701  LOGICAL :: reverse
702  REAL(kind=dp) :: a, a2, cval, fb
703  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc, rhoc, rval
704 
705  my_typ = "RHO"
706  IF (PRESENT(typ)) my_typ = typ(1:3)
707 
708  SELECT CASE (potential%ppot_type)
709  CASE (no_pseudo, ecp_pseudo)
710  ! we do nothing
711  CASE (gth_pseudo)
712  IF (potential%gth_pot%nlcc) THEN
713  m = SIZE(corden)
714  ALLOCATE (fe(m), rc(m))
715  n = potential%gth_pot%nexp_nlcc
716  DO i = 1, n
717  a = potential%gth_pot%alpha_nlcc(i)
718  a2 = a*a
719  ! note that all terms are computed with rc, not rr
720  rc(:) = rr(:)/a
721  fe(:) = exp(-0.5_dp*rc(:)*rc(:))
722  DO j = 1, potential%gth_pot%nct_nlcc(i)
723  cval = potential%gth_pot%cval_nlcc(j, i)
724  IF (my_typ == "RHO") THEN
725  corden(:) = corden(:) + fe(:)*rc**(2*j - 2)*cval
726  ELSE IF (my_typ == "DER") THEN
727  corden(:) = corden(:) - fe(:)*rc**(2*j - 1)*cval/a
728  IF (j > 1) THEN
729  corden(:) = corden(:) + real(2*j - 2, dp)*fe(:)*rc**(2*j - 3)*cval/a
730  END IF
731  ELSE IF (my_typ == "LAP") THEN
732  fb = 2._dp*cval/a
733  corden(:) = corden(:) - fb*fe(:)/rr(:)*rc**(2*j - 1)
734  corden(:) = corden(:) + fe(:)*rc**(2*j)*cval/a2
735  IF (j > 1) THEN
736  corden(:) = corden(:) + fb*real(2*j - 2, dp)*fe(:)/rr(:)*rc**(2*j - 3)
737  corden(:) = corden(:) + real((2*j - 2)*(2*j - 3), dp)*fe(:)*rc**(2*j - 4)*cval/a2
738  corden(:) = corden(:) - real(2*j - 2, dp)*fe(:)*rc**(2*j - 2)*cval/a2
739  END IF
740  ELSE
741  cpabort("")
742  END IF
743  END DO
744  END DO
745  DEALLOCATE (fe, rc)
746  END IF
747  CASE (upf_pseudo)
748  IF (potential%upf_pot%core_correction) THEN
749  m = SIZE(corden)
750  n = potential%upf_pot%mesh_size
751  reverse = .false.
752  IF (rr(1) > rr(m)) reverse = .true.
753  ALLOCATE (rhoc(m), rval(m))
754  IF (reverse) THEN
755  DO i = 1, m
756  rval(i) = rr(m - i + 1)
757  END DO
758  ELSE
759  rval(1:m) = rr(1:m)
760  END IF
761  IF (my_typ == "RHO") THEN
762  CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
763  ynew=rhoc(1:m))
764  ELSE IF (my_typ == "DER") THEN
765  CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
766  dynew=rhoc(1:m))
767  ELSE IF (my_typ == "LAP") THEN
768  CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
769  d2ynew=rhoc(1:m))
770  ELSE
771  cpabort("")
772  END IF
773  IF (reverse) THEN
774  DO i = 1, m
775  rval(i) = rr(m - i + 1)
776  corden(i) = corden(i) + rhoc(m - i + 1)
777  END DO
778  ELSE
779  corden(1:m) = corden(1:m) + rhoc(1:m)
780  END IF
781  DEALLOCATE (rhoc, rval)
782  END IF
783  CASE (sgp_pseudo)
784  IF (potential%sgp_pot%has_nlcc) THEN
785  cpabort("not implemented")
786  END IF
787  CASE DEFAULT
788  cpabort("Unknown PP type")
789  END SELECT
790 
791  END SUBROUTINE atom_core_density
792 
793 ! **************************************************************************************************
794 !> \brief ...
795 !> \param locpot ...
796 !> \param gthpot ...
797 !> \param rr ...
798 ! **************************************************************************************************
799  PURE SUBROUTINE atom_local_potential(locpot, gthpot, rr)
800  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: locpot
801  TYPE(atom_gthpot_type), INTENT(IN) :: gthpot
802  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rr
803 
804  INTEGER :: i, j, m, n
805  REAL(kind=dp) :: a
806  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc
807 
808  m = SIZE(locpot)
809  ALLOCATE (fe(m), rc(m))
810  rc(:) = rr(:)/gthpot%rc
811  DO i = 1, m
812  locpot(i) = -gthpot%zion*erf(rc(i)/sqrt(2._dp))/rr(i)
813  END DO
814  n = gthpot%ncl
815  fe(:) = exp(-0.5_dp*rc(:)*rc(:))
816  DO i = 1, n
817  locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cl(i)
818  END DO
819  IF (gthpot%lpotextended) THEN
820  DO j = 1, gthpot%nexp_lpot
821  a = gthpot%alpha_lpot(j)
822  rc(:) = rr(:)/a
823  fe(:) = exp(-0.5_dp*rc(:)*rc(:))
824  n = gthpot%nct_lpot(j)
825  DO i = 1, n
826  locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cval_lpot(i, j)
827  END DO
828  END DO
829  END IF
830  DEALLOCATE (fe, rc)
831 
832  END SUBROUTINE atom_local_potential
833 
834 ! **************************************************************************************************
835 !> \brief ...
836 !> \param rmax ...
837 !> \param wfn ...
838 !> \param rcov ...
839 !> \param l ...
840 !> \param basis ...
841 ! **************************************************************************************************
842  PURE SUBROUTINE atom_orbital_max(rmax, wfn, rcov, l, basis)
843  REAL(kind=dp), INTENT(OUT) :: rmax
844  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wfn
845  REAL(kind=dp), INTENT(IN) :: rcov
846  INTEGER, INTENT(IN) :: l
847  TYPE(atom_basis_type), INTENT(IN) :: basis
848 
849  INTEGER :: i, m, n
850  REAL(kind=dp) :: ff
851  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dorb
852 
853  m = SIZE(basis%bf, 1)
854  ALLOCATE (dorb(m))
855  n = basis%nbas(l)
856  dorb = 0._dp
857  DO i = 1, n
858  ff = wfn(i)
859  dorb(1:m) = dorb(1:m) + ff*basis%dbf(1:m, i, l)
860  END DO
861  rmax = -1._dp
862  DO i = 1, m - 1
863  IF (basis%grid%rad(i) < 2*rcov) THEN
864  IF (dorb(i)*dorb(i + 1) < 0._dp) THEN
865  rmax = max(rmax, basis%grid%rad(i))
866  END IF
867  END IF
868  END DO
869  DEALLOCATE (dorb)
870 
871  END SUBROUTINE atom_orbital_max
872 
873 ! **************************************************************************************************
874 !> \brief ...
875 !> \param node ...
876 !> \param wfn ...
877 !> \param rcov ...
878 !> \param l ...
879 !> \param basis ...
880 ! **************************************************************************************************
881  PURE SUBROUTINE atom_orbital_nodes(node, wfn, rcov, l, basis)
882  INTEGER, INTENT(OUT) :: node
883  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wfn
884  REAL(kind=dp), INTENT(IN) :: rcov
885  INTEGER, INTENT(IN) :: l
886  TYPE(atom_basis_type), INTENT(IN) :: basis
887 
888  INTEGER :: i, m, n
889  REAL(kind=dp) :: ff
890  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: orb
891 
892  node = 0
893  m = SIZE(basis%bf, 1)
894  ALLOCATE (orb(m))
895  n = basis%nbas(l)
896  orb = 0._dp
897  DO i = 1, n
898  ff = wfn(i)
899  orb(1:m) = orb(1:m) + ff*basis%bf(1:m, i, l)
900  END DO
901  DO i = 1, m - 1
902  IF (basis%grid%rad(i) < rcov) THEN
903  IF (orb(i)*orb(i + 1) < 0._dp) node = node + 1
904  END IF
905  END DO
906  DEALLOCATE (orb)
907 
908  END SUBROUTINE atom_orbital_nodes
909 
910 ! **************************************************************************************************
911 !> \brief ...
912 !> \param value ...
913 !> \param wfn ...
914 !> \param basis ...
915 ! **************************************************************************************************
916  PURE SUBROUTINE atom_wfnr0(value, wfn, basis)
917  REAL(kind=dp), INTENT(OUT) :: value
918  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: wfn
919  TYPE(atom_basis_type), INTENT(IN) :: basis
920 
921  INTEGER :: i, m, n
922 
923  value = 0._dp
924  m = maxval(minloc(basis%grid%rad))
925  n = basis%nbas(0)
926  DO i = 1, n
927  value = value + wfn(i)*basis%bf(m, i, 0)
928  END DO
929  END SUBROUTINE atom_wfnr0
930 
931 ! **************************************************************************************************
932 !> \brief Solve the generalised eigenproblem for every angular momentum.
933 !> \param hmat Hamiltonian matrix
934 !> \param umat transformation matrix which reduces the overlap matrix to its unitary form
935 !> \param orb atomic wavefunctions
936 !> \param ener atomic orbital energies
937 !> \param nb number of contracted basis functions
938 !> \param nv number of linear-independent contracted basis functions
939 !> \param maxl maximum angular momentum to consider
940 !> \par History
941 !> * 08.2008 created [Juerg Hutter]
942 ! **************************************************************************************************
943  SUBROUTINE atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
944  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: hmat, umat
945  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: orb
946  REAL(kind=dp), DIMENSION(:, 0:), INTENT(INOUT) :: ener
947  INTEGER, DIMENSION(0:), INTENT(IN) :: nb, nv
948  INTEGER, INTENT(IN) :: maxl
949 
950  INTEGER :: info, l, lwork, m, n
951  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
952  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a
953 
954  cpassert(all(nb >= nv))
955 
956  orb = 0._dp
957  DO l = 0, maxl
958  n = nb(l)
959  m = nv(l)
960  IF (n > 0 .AND. m > 0) THEN
961  lwork = 10*m
962  ALLOCATE (a(n, n), w(n), work(lwork))
963  a(1:m, 1:m) = matmul(transpose(umat(1:n, 1:m, l)), matmul(hmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
964  CALL lapack_ssyev("V", "U", m, a(1:m, 1:m), m, w(1:m), work, lwork, info)
965  a(1:n, 1:m) = matmul(umat(1:n, 1:m, l), a(1:m, 1:m))
966 
967  m = min(m, SIZE(orb, 2))
968  orb(1:n, 1:m, l) = a(1:n, 1:m)
969  ener(1:m, l) = w(1:m)
970 
971  DEALLOCATE (a, w, work)
972  END IF
973  END DO
974 
975  END SUBROUTINE atom_solve
976 
977 ! **************************************************************************************************
978 !> \brief THIS FUNCTION IS NO LONGER IN USE.
979 !> \param fun ...
980 !> \param deps ...
981 !> \return ...
982 ! **************************************************************************************************
983  FUNCTION prune_grid(fun, deps) RESULT(nc)
984  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: fun
985  REAL(kind=dp), INTENT(IN), OPTIONAL :: deps
986  INTEGER :: nc
987 
988  INTEGER :: i, nr
989  REAL(kind=dp) :: meps
990 
991  meps = 1.e-14_dp
992  IF (PRESENT(deps)) meps = deps
993 
994  nr = SIZE(fun)
995  nc = 0
996  DO i = nr, 1, -1
997  IF (abs(fun(i)) > meps) THEN
998  nc = i
999  EXIT
1000  END IF
1001  END DO
1002 
1003  END FUNCTION prune_grid
1004 
1005 ! **************************************************************************************************
1006 !> \brief Integrate a function on an atomic radial grid.
1007 !> \param fun function to integrate
1008 !> \param grid atomic radial grid
1009 !> \return value of the integral
1010 !> \par History
1011 !> * 08.2008 created [Juerg Hutter]
1012 ! **************************************************************************************************
1013  PURE FUNCTION integrate_grid_function1(fun, grid) RESULT(integral)
1014  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: fun
1015  TYPE(grid_atom_type), INTENT(IN) :: grid
1016  REAL(kind=dp) :: integral
1017 
1018  INTEGER :: nc
1019 
1020  nc = SIZE(fun)
1021  integral = sum(fun(1:nc)*grid%wr(1:nc))
1022 
1023  END FUNCTION integrate_grid_function1
1024 
1025 ! **************************************************************************************************
1026 !> \brief Integrate the product of two functions on an atomic radial grid.
1027 !> \param fun1 first function
1028 !> \param fun2 second function
1029 !> \param grid atomic radial grid
1030 !> \return value of the integral
1031 !> \par History
1032 !> * 08.2008 created [Juerg Hutter]
1033 ! **************************************************************************************************
1034  PURE FUNCTION integrate_grid_function2(fun1, fun2, grid) RESULT(integral)
1035  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2
1036  TYPE(grid_atom_type), INTENT(IN) :: grid
1037  REAL(kind=dp) :: integral
1038 
1039  INTEGER :: nc
1040 
1041  nc = min(SIZE(fun1), SIZE(fun2))
1042  integral = sum(fun1(1:nc)*fun2(1:nc)*grid%wr(1:nc))
1043 
1044  END FUNCTION integrate_grid_function2
1045 
1046 ! **************************************************************************************************
1047 !> \brief Integrate the product of three functions on an atomic radial grid.
1048 !> \param fun1 first function
1049 !> \param fun2 second function
1050 !> \param fun3 third function
1051 !> \param grid atomic radial grid
1052 !> \return value of the integral
1053 !> \par History
1054 !> * 08.2008 created [Juerg Hutter]
1055 ! **************************************************************************************************
1056  PURE FUNCTION integrate_grid_function3(fun1, fun2, fun3, grid) RESULT(integral)
1057  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2, fun3
1058  TYPE(grid_atom_type), INTENT(IN) :: grid
1059  REAL(kind=dp) :: integral
1060 
1061  INTEGER :: nc
1062 
1063  nc = min(SIZE(fun1), SIZE(fun2), SIZE(fun3))
1064  integral = sum(fun1(1:nc)*fun2(1:nc)*fun3(1:nc)*grid%wr(1:nc))
1065 
1066  END FUNCTION integrate_grid_function3
1067 
1068 ! **************************************************************************************************
1069 !> \brief Numerically compute the Coulomb potential on an atomic radial grid.
1070 !> \param cpot Coulomb potential on the radial grid
1071 !> \param density electron density on the radial grid
1072 !> \param grid atomic radial grid
1073 !> \par History
1074 !> * 08.2008 created [Juerg Hutter]
1075 ! **************************************************************************************************
1076  SUBROUTINE coulomb_potential_numeric(cpot, density, grid)
1077  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1078  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density
1079  TYPE(grid_atom_type), INTENT(IN) :: grid
1080 
1081  INTEGER :: i, nc
1082  REAL(dp) :: int1, int2
1083  REAL(dp), DIMENSION(:), POINTER :: r, wr
1084 
1085  nc = min(SIZE(cpot), SIZE(density))
1086  r => grid%rad
1087  wr => grid%wr
1088 
1089  int1 = fourpi*integrate_grid(density, grid)
1090  int2 = 0._dp
1091  cpot(nc:) = int1/r(nc:)
1092  ! IF (log_unit>0) WRITE(log_unit,"(A,2F10.8)") "r", int1, cpot(nc:)
1093 
1094  ! test that grid is decreasing
1095  cpassert(r(1) > r(nc))
1096  DO i = 1, nc
1097  cpot(i) = int1/r(i) + int2
1098  int1 = int1 - fourpi*density(i)*wr(i)
1099  int2 = int2 + fourpi*density(i)*wr(i)/r(i)
1100  END DO
1101 
1102  END SUBROUTINE coulomb_potential_numeric
1103 
1104 ! **************************************************************************************************
1105 !> \brief Analytically compute the Coulomb potential on an atomic radial grid.
1106 !> \param cpot Coulomb potential on the radial grid
1107 !> \param pmat density matrix
1108 !> \param basis atomic basis set
1109 !> \param grid atomic radial grid
1110 !> \param maxl maximum angular momentum to consider
1111 !> \par History
1112 !> * 08.2008 created [Juerg Hutter]
1113 ! **************************************************************************************************
1114  SUBROUTINE coulomb_potential_analytic(cpot, pmat, basis, grid, maxl)
1115  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1116  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1117  TYPE(atom_basis_type), INTENT(IN) :: basis
1118  TYPE(grid_atom_type) :: grid
1119  INTEGER, INTENT(IN) :: maxl
1120 
1121  INTEGER :: i, j, k, l, m, n
1122  REAL(kind=dp) :: a, b, ff, oab, sab
1123  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: erfa, expa, z
1124  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: unp
1125 
1126  m = SIZE(cpot)
1127  ALLOCATE (erfa(1:m), expa(1:m), z(1:m))
1128 
1129  cpot = 0._dp
1130 
1131  DO l = 0, maxl
1132  IF (maxval(abs(pmat(:, :, l))) < 1.e-14_dp) cycle
1133  SELECT CASE (basis%basis_type)
1134  CASE DEFAULT
1135  cpabort("")
1136  CASE (gto_basis)
1137  DO i = 1, basis%nbas(l)
1138  DO j = i, basis%nbas(l)
1139  IF (abs(pmat(i, j, l)) < 1.e-14_dp) cycle
1140  ff = pmat(i, j, l)
1141  IF (i /= j) ff = 2._dp*ff
1142  a = basis%am(i, l)
1143  b = basis%am(j, l)
1144  sab = sqrt(a + b)
1145  oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1146  z(:) = sab*grid%rad(:)
1147  DO k = 1, SIZE(erfa)
1148  erfa(k) = oab*erf(z(k))/grid%rad(k)
1149  END DO
1150  expa(:) = exp(-z(:)**2)*ff/(a + b)**(l + 1)
1151  SELECT CASE (l)
1152  CASE DEFAULT
1153  cpabort("")
1154  CASE (0)
1155  cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1156  CASE (1)
1157  cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1158  CASE (2)
1159  cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1160  CASE (3)
1161  cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1162  END SELECT
1163  END DO
1164  END DO
1165  CASE (cgto_basis)
1166  n = basis%nprim(l)
1167  m = basis%nbas(l)
1168  ALLOCATE (unp(n, n))
1169 
1170  unp(1:n, 1:n) = matmul(matmul(basis%cm(1:n, 1:m, l), pmat(1:m, 1:m, l)), &
1171  transpose(basis%cm(1:n, 1:m, l)))
1172  DO i = 1, basis%nprim(l)
1173  DO j = i, basis%nprim(l)
1174  IF (abs(unp(i, j)) < 1.e-14_dp) cycle
1175  ff = unp(i, j)
1176  IF (i /= j) ff = 2._dp*ff
1177  a = basis%am(i, l)
1178  b = basis%am(j, l)
1179  sab = sqrt(a + b)
1180  oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1181  z(:) = sab*grid%rad(:)
1182  DO k = 1, SIZE(erfa)
1183  erfa(k) = oab*erf(z(k))/grid%rad(k)
1184  END DO
1185  expa(:) = exp(-z(:)**2)*ff/(a + b)**(l + 1)
1186  SELECT CASE (l)
1187  CASE DEFAULT
1188  cpabort("")
1189  CASE (0)
1190  cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1191  CASE (1)
1192  cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1193  CASE (2)
1194  cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1195  CASE (3)
1196  cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1197  END SELECT
1198  END DO
1199  END DO
1200 
1201  DEALLOCATE (unp)
1202  END SELECT
1203  END DO
1204  DEALLOCATE (erfa, expa, z)
1205 
1206  END SUBROUTINE coulomb_potential_analytic
1207 
1208 ! **************************************************************************************************
1209 !> \brief Calculate the exchange potential numerically.
1210 !> \param kmat Kohn-Sham matrix
1211 !> \param state electronic state
1212 !> \param occ occupation numbers
1213 !> \param wfn wavefunctions
1214 !> \param basis atomic basis set
1215 !> \param hfx_pot potential parameters from Hartree-Fock section
1216 !> \par History
1217 !> * 01.2009 created [Juerg Hutter]
1218 ! **************************************************************************************************
1219  SUBROUTINE exchange_numeric(kmat, state, occ, wfn, basis, hfx_pot)
1220  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1221  TYPE(atom_state), INTENT(IN) :: state
1222  REAL(kind=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1223  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wfn
1224  TYPE(atom_basis_type), INTENT(IN) :: basis
1225  TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1226 
1227  INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1228  norb, nr, nu
1229  REAL(kind=dp) :: almn
1230  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi, pot
1231  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1232  REAL(kind=dp), DIMENSION(0:maxfac) :: arho
1233 
1234  arho = 0._dp
1235  DO ll = 0, maxfac, 2
1236  lh = ll/2
1237  arho(ll) = fac(ll)/fac(lh)**2
1238  END DO
1239 
1240  kmat = 0._dp
1241 
1242  nr = basis%grid%nr
1243  ALLOCATE (nai(nr), nbi(nr), cpot(nr), pot(nr))
1244 
1245  DO lad = 0, state%maxl_calc
1246  DO lbc = 0, state%maxl_occ
1247  norb = state%maxn_occ(lbc)
1248  nbas = basis%nbas(lbc)
1249  ! calculate orbitals for angmom lbc
1250  ALLOCATE (orb(nr, norb))
1251  orb = 0._dp
1252  DO i = 1, norb
1253  DO k = 1, nbas
1254  orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1255  END DO
1256  END DO
1257  DO nu = abs(lad - lbc), lad + lbc, 2
1258  almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu)*arho(lad + lbc - nu)/(real(lad + lbc + nu + 1, dp) &
1259  *arho(lad + lbc + nu))
1260  almn = -0.5_dp*almn
1261 
1262  DO ia = 1, basis%nbas(lad)
1263  DO i = 1, norb
1264  nai(:) = orb(:, i)*basis%bf(:, ia, lad)
1265  cpot = 0.0_dp
1266  IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1267  CALL potential_coulomb_numeric(pot, nai, nu, basis%grid)
1268  cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_coulomb
1269  END IF
1270  IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1271  IF (hfx_pot%do_gh) THEN
1272  CALL potential_longrange_numeric_gh(pot, nai, nu, basis%grid, hfx_pot%omega, &
1273  hfx_pot%kernel(:, :, nu))
1274  ELSE
1275  CALL potential_longrange_numeric(pot, nai, nu, basis%grid, hfx_pot%omega, &
1276  hfx_pot%kernel(:, :, nu))
1277  END IF
1278  cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_longrange
1279  END IF
1280  DO ib = 1, basis%nbas(lad)
1281  kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1282  integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1283  END DO
1284  END DO
1285  END DO
1286 
1287  END DO
1288  DEALLOCATE (orb)
1289  END DO
1290  END DO
1291 
1292  DEALLOCATE (nai, nbi, cpot)
1293 
1294  END SUBROUTINE exchange_numeric
1295 
1296 ! **************************************************************************************************
1297 !> \brief Calculate the exchange potential semi-analytically.
1298 !> \param kmat Kohn-Sham matrix
1299 !> \param state electronic state
1300 !> \param occ occupation numbers
1301 !> \param wfn wavefunctions
1302 !> \param basis atomic basis set
1303 !> \param hfx_pot properties of the Hartree-Fock potential
1304 !> \par History
1305 !> * 01.2009 created [Juerg Hutter]
1306 ! **************************************************************************************************
1307  SUBROUTINE exchange_semi_analytic(kmat, state, occ, wfn, basis, hfx_pot)
1308  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1309  TYPE(atom_state), INTENT(IN) :: state
1310  REAL(kind=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1311  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: wfn
1312  TYPE(atom_basis_type), INTENT(IN) :: basis
1313  TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1314 
1315  INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1316  norb, nr, nu
1317  REAL(kind=dp) :: almn
1318  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi
1319  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1320  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pot
1321  REAL(kind=dp), DIMENSION(0:maxfac) :: arho
1322 
1323  arho = 0._dp
1324  DO ll = 0, maxfac, 2
1325  lh = ll/2
1326  arho(ll) = fac(ll)/fac(lh)**2
1327  END DO
1328 
1329  kmat = 0._dp
1330 
1331  nr = basis%grid%nr
1332  nbas = maxval(basis%nbas)
1333  ALLOCATE (pot(nr, nbas, nbas))
1334  ALLOCATE (nai(nr), nbi(nr), cpot(nr))
1335 
1336  DO lad = 0, state%maxl_calc
1337  DO lbc = 0, state%maxl_occ
1338  norb = state%maxn_occ(lbc)
1339  nbas = basis%nbas(lbc)
1340  ! calculate orbitals for angmom lbc
1341  ALLOCATE (orb(nr, norb))
1342  orb = 0._dp
1343  DO i = 1, norb
1344  DO k = 1, nbas
1345  orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1346  END DO
1347  END DO
1348  DO nu = abs(lad - lbc), lad + lbc, 2
1349  almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu) &
1350  *arho(lad + lbc - nu)/(real(lad + lbc + nu + 1, dp)*arho(lad + lbc + nu))
1351  almn = -0.5_dp*almn
1352  ! calculate potential for basis function pair (lad,lbc)
1353  pot = 0._dp
1354  CALL potential_analytic(pot, lad, lbc, nu, basis, hfx_pot)
1355  DO ia = 1, basis%nbas(lad)
1356  DO i = 1, norb
1357  cpot = 0._dp
1358  DO k = 1, nbas
1359  cpot(:) = cpot(:) + pot(:, ia, k)*wfn(k, i, lbc)
1360  END DO
1361  DO ib = 1, basis%nbas(lad)
1362  kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1363  integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1364  END DO
1365  END DO
1366  END DO
1367  END DO
1368  DEALLOCATE (orb)
1369  END DO
1370  END DO
1371 
1372  DEALLOCATE (pot)
1373  DEALLOCATE (nai, nbi, cpot)
1374 
1375  END SUBROUTINE exchange_semi_analytic
1376 
1377 ! **************************************************************************************************
1378 !> \brief Calculate the electrostatic potential using numerical differentiation.
1379 !> \param cpot computed potential
1380 !> \param density electron density on the atomic radial grid
1381 !> \param nu integer parameter [ABS(la-lb) .. la+lb];
1382 !> see potential_analytic() and exchange_numeric()
1383 !> \param grid atomic radial grid
1384 !> \par History
1385 !> * 01.2009 created [Juerg Hutter]
1386 ! **************************************************************************************************
1387  SUBROUTINE potential_coulomb_numeric(cpot, density, nu, grid)
1388  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cpot
1389  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density
1390  INTEGER, INTENT(IN) :: nu
1391  TYPE(grid_atom_type), INTENT(IN) :: grid
1392 
1393  INTEGER :: i, nc
1394  REAL(dp) :: int1, int2
1395  REAL(dp), DIMENSION(:), POINTER :: r, wr
1396 
1397  nc = min(SIZE(cpot), SIZE(density))
1398  r => grid%rad
1399  wr => grid%wr
1400 
1401  int1 = integrate_grid(density, r**nu, grid)
1402  int2 = 0._dp
1403  cpot(nc:) = int1/r(nc:)**(nu + 1)
1404 
1405  ! test that grid is decreasing
1406  cpassert(r(1) > r(nc))
1407  DO i = 1, nc
1408  cpot(i) = int1/r(i)**(nu + 1) + int2*r(i)**nu
1409  int1 = int1 - r(i)**(nu)*density(i)*wr(i)
1410  int2 = int2 + density(i)*wr(i)/r(i)**(nu + 1)
1411  END DO
1412 
1413  END SUBROUTINE potential_coulomb_numeric
1414 
1415 ! **************************************************************************************************
1416 !> \brief ...
1417 !> \param cpot ...
1418 !> \param density ...
1419 !> \param nu ...
1420 !> \param grid ...
1421 !> \param omega ...
1422 !> \param kernel ...
1423 ! **************************************************************************************************
1424  SUBROUTINE potential_longrange_numeric(cpot, density, nu, grid, omega, kernel)
1425  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cpot
1426  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density
1427  INTEGER, INTENT(IN) :: nu
1428  TYPE(grid_atom_type), INTENT(IN) :: grid
1429  REAL(kind=dp), INTENT(IN) :: omega
1430  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1431 
1432  INTEGER :: nc
1433  REAL(dp), DIMENSION(:), POINTER :: r, wr
1434  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1435 
1436  nc = min(SIZE(cpot), SIZE(density))
1437  r => grid%rad
1438  wr => grid%wr
1439 
1440  ALLOCATE (work_inp(nc), work_out(nc))
1441 
1442  cpot = 0.0_dp
1443 
1444  ! First Bessel transform
1445  work_inp(:nc) = density(:nc)*wr(:nc)
1446  CALL dsymv('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1447 
1448  ! Second Bessel transform
1449  work_inp(:nc) = work_out(:nc)*exp(-r(:nc)**2)/r(:nc)**2*wr(:nc)
1450  CALL dsymv('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1451 
1452  cpot(:nc) = work_out(:nc)*(2.0_dp*real(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1453 
1454  END SUBROUTINE potential_longrange_numeric
1455 
1456 ! **************************************************************************************************
1457 !> \brief ...
1458 !> \param cpot ...
1459 !> \param density ...
1460 !> \param nu ...
1461 !> \param grid ...
1462 !> \param omega ...
1463 !> \param kernel ...
1464 ! **************************************************************************************************
1465  SUBROUTINE potential_longrange_numeric_gh(cpot, density, nu, grid, omega, kernel)
1466  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cpot
1467  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density
1468  INTEGER, INTENT(IN) :: nu
1469  TYPE(grid_atom_type), INTENT(IN) :: grid
1470  REAL(kind=dp), INTENT(IN) :: omega
1471  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1472 
1473  INTEGER :: n_max, nc, nc_kernel, nr_kernel
1474  REAL(dp), DIMENSION(:), POINTER :: wr
1475  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1476 
1477  nc = min(SIZE(cpot), SIZE(density))
1478  wr => grid%wr
1479 
1480  nc_kernel = SIZE(kernel, 1)
1481  nr_kernel = SIZE(kernel, 2)
1482  n_max = max(nc, nc_kernel, nr_kernel)
1483 
1484  ALLOCATE (work_inp(n_max), work_out(n_max))
1485 
1486  cpot = 0.0_dp
1487 
1488  ! First Bessel transform
1489  work_inp(:nc) = density(:nc)*wr(:nc)
1490  CALL dgemv('T', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1491 
1492  ! Second Bessel transform
1493  work_inp(:nr_kernel) = work_out(:nr_kernel)
1494  CALL dgemv('N', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1495 
1496  cpot(:nc) = work_out(:nc)*(2.0_dp*real(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1497 
1498  END SUBROUTINE potential_longrange_numeric_gh
1499 
1500 ! **************************************************************************************************
1501 !> \brief Calculate the electrostatic potential using analytical expressions.
1502 !> The interaction potential has the form
1503 !> V(r)=scale_coulomb*1/r+scale_lr*erf(omega*r)/r
1504 !> \param cpot computed potential
1505 !> \param la angular momentum of the calculated state
1506 !> \param lb angular momentum of the occupied state
1507 !> \param nu integer parameter [ABS(la-lb) .. la+lb] with the parity of 'la+lb'
1508 !> \param basis atomic basis set
1509 !> \param hfx_pot properties of the Hartree-Fock potential
1510 !> \par History
1511 !> * 01.2009 created [Juerg Hutter]
1512 ! **************************************************************************************************
1513  SUBROUTINE potential_analytic(cpot, la, lb, nu, basis, hfx_pot)
1514  REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT) :: cpot
1515  INTEGER, INTENT(IN) :: la, lb, nu
1516  TYPE(atom_basis_type), INTENT(IN) :: basis
1517  TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1518 
1519  INTEGER :: i, j, k, l, ll, m
1520  REAL(kind=dp) :: a, b
1521  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: erfa, pot
1522 
1523  m = SIZE(cpot, 1)
1524  ALLOCATE (erfa(1:m))
1525 
1526  ll = la + lb
1527 
1528  cpot = 0._dp
1529 
1530  SELECT CASE (basis%basis_type)
1531  CASE DEFAULT
1532  cpabort("")
1533  CASE (gto_basis)
1534  DO i = 1, basis%nbas(la)
1535  DO j = 1, basis%nbas(lb)
1536  a = basis%am(i, la)
1537  b = basis%am(j, lb)
1538 
1539  IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1540  CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1541 
1542  cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_coulomb
1543  END IF
1544 
1545  IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1546  CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1547 
1548  cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_longrange
1549  END IF
1550  END DO
1551  END DO
1552  CASE (cgto_basis)
1553  ALLOCATE (pot(1:m))
1554 
1555  DO i = 1, basis%nprim(la)
1556  DO j = 1, basis%nprim(lb)
1557  a = basis%am(i, la)
1558  b = basis%am(j, lb)
1559 
1560  pot = 0.0_dp
1561 
1562  IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1563  CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1564 
1565  pot(:) = pot(:) + erfa(:)*hfx_pot%scale_coulomb
1566  END IF
1567 
1568  IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1569  CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1570 
1571  pot(:) = pot(:) + erfa(:)*hfx_pot%scale_longrange
1572  END IF
1573 
1574  DO k = 1, basis%nbas(la)
1575  DO l = 1, basis%nbas(lb)
1576  cpot(:, k, l) = cpot(:, k, l) + pot(:)*basis%cm(i, k, la)*basis%cm(j, l, lb)
1577  END DO
1578  END DO
1579  END DO
1580  END DO
1581 
1582  END SELECT
1583 
1584  DEALLOCATE (erfa)
1585 
1586  END SUBROUTINE potential_analytic
1587 
1588 ! **************************************************************************************************
1589 !> \brief ...
1590 !> \param erfa Array will contain the potential
1591 !> \param a Exponent of first Gaussian charge distribution
1592 !> \param b Exponent of second Gaussian charge distribution
1593 !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1594 !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1595 !> \param rad Radial grid
1596 ! **************************************************************************************************
1597  SUBROUTINE potential_coulomb_analytic(erfa, a, b, ll, nu, rad)
1598  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: erfa
1599  REAL(kind=dp), INTENT(IN) :: a, b
1600  INTEGER, INTENT(IN) :: ll, nu
1601  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rad
1602 
1603  INTEGER :: nr
1604  REAL(kind=dp) :: oab, sab
1605  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1606 
1607  nr = SIZE(rad)
1608  ALLOCATE (expa(nr), z(nr))
1609 
1610  sab = sqrt(a + b)
1611  oab = dfac(ll + nu + 1)*rootpi/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1612  z(:) = sab*rad(:)
1613  erfa(:) = oab*erf(z(:))/z(:)**(nu + 1)
1614  expa(:) = exp(-z(:)**2)/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1615  SELECT CASE (ll)
1616  CASE DEFAULT
1617  cpabort("")
1618  CASE (0)
1619  cpassert(nu == 0)
1620  CASE (1)
1621  cpassert(nu == 1)
1622  erfa(:) = erfa(:) - 6._dp*expa(:)/z(:)
1623  CASE (2)
1624  SELECT CASE (nu)
1625  CASE DEFAULT
1626  cpabort("")
1627  CASE (0)
1628  erfa(:) = erfa(:) - 2._dp*expa(:)
1629  CASE (2)
1630  erfa(:) = erfa(:) - expa(:)*(20._dp + 30._dp/z(:)**2)
1631  END SELECT
1632  CASE (3)
1633  SELECT CASE (nu)
1634  CASE DEFAULT
1635  cpabort("")
1636  CASE (1)
1637  erfa(:) = erfa(:) - expa(:)*(12._dp*z(:) + 30._dp/z(:))
1638  CASE (3)
1639  erfa(:) = erfa(:) - expa(:)*(56._dp*z(:) + 140._dp/z(:) + 210._dp/z(:)**3)
1640  END SELECT
1641  CASE (4)
1642  SELECT CASE (nu)
1643  CASE DEFAULT
1644  cpabort("")
1645  CASE (0)
1646  erfa(:) = erfa(:) - expa(:)*(4._dp*z(:)**2 + 14._dp)
1647  CASE (2)
1648  erfa(:) = erfa(:) - expa(:)*(40._dp*z(:)**2 + 140._dp + 210._dp/z(:)**2)
1649  CASE (4)
1650  erfa(:) = erfa(:) - expa(:)*(144._dp*z(:)**2 + 504._dp + 1260._dp/z(:)**2 + 1890._dp/z(:)**4)
1651  END SELECT
1652  CASE (5)
1653  SELECT CASE (nu)
1654  CASE DEFAULT
1655  cpabort("")
1656  CASE (1)
1657  erfa(:) = erfa(:) - expa(:)*(24._dp*z(:)**3 + 108._dp*z(:) + 210._dp/z(:))
1658  CASE (3)
1659  erfa(:) = erfa(:) - expa(:)*(112._dp*z(:)**3 + 504._dp*z(:) + 1260._dp/z(:) + 1890._dp/z(:)**3)
1660  CASE (5)
1661  erfa(:) = erfa(:) - expa(:)*(352._dp*z(:)**3 + 1584._dp*z(:) + 5544._dp/z(:) + &
1662  13860._dp/z(:)**3 + 20790._dp/z(:)**5)
1663  END SELECT
1664  CASE (6)
1665  SELECT CASE (nu)
1666  CASE DEFAULT
1667  cpabort("")
1668  CASE (0)
1669  erfa(:) = erfa(:) - expa(:)*(8._dp*z(:)**4 + 44._dp*z(:)**2 + 114._dp)
1670  CASE (2)
1671  erfa(:) = erfa(:) - expa(:)*(80._dp*z(:)**4 + 440._dp*z(:)**2 + 1260._dp + 1896._dp/z(:)**2)
1672  CASE (4)
1673  erfa(:) = erfa(:) - expa(:)*(288._dp*z(:)**4 + 1584._dp*z(:)**2 + 5544._dp + &
1674  13860._dp/z(:)**2 + 20790._dp/z(:)**4)
1675  CASE (6)
1676  erfa(:) = erfa(:) - expa(:)*(832._dp*z(:)**4 + 4576._dp*z(:)**2 + 20592._dp + &
1677  72072._dp/z(:)**2 + 180180._dp/z(:)**4 + 270270._dp/z(:)**6)
1678  END SELECT
1679  END SELECT
1680 
1681  DEALLOCATE (expa, z)
1682 
1683  END SUBROUTINE potential_coulomb_analytic
1684 
1685 ! **************************************************************************************************
1686 !> \brief This routine calculates the longrange Coulomb potential of a product of two Gaussian with given angular momentum
1687 !> \param erfa Array will contain the potential
1688 !> \param a Exponent of first Gaussian charge distribution
1689 !> \param b Exponent of second Gaussian charge distribution
1690 !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1691 !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1692 !> \param rad Radial grid
1693 !> \param omega Range-separation parameter
1694 ! **************************************************************************************************
1695  PURE SUBROUTINE potential_longrange_analytic(erfa, a, b, ll, nu, rad, omega)
1696  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: erfa
1697  REAL(kind=dp), INTENT(IN) :: a, b
1698  INTEGER, INTENT(IN) :: ll, nu
1699  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rad
1700  REAL(kind=dp), INTENT(IN) :: omega
1701 
1702  INTEGER :: i, lambda, nr
1703  REAL(kind=dp) :: ab, oab, pab, prel, sab
1704  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1705 
1706  IF (mod(ll - nu, 2) == 0 .AND. ll >= nu .AND. nu >= 0) THEN
1707  nr = SIZE(rad)
1708  ALLOCATE (expa(nr), z(nr))
1709 
1710  lambda = (ll - nu)/2
1711  ab = a + b
1712  sab = sqrt(ab)
1713  pab = omega**2*ab/(omega**2 + ab)
1714  prel = pab/ab
1715  z(:) = sqrt(pab)*rad(:)
1716  oab = fac(lambda)/sqrt(ab)**(ll + 2)/4.0_dp*sqrt(prel)**(nu + 1)*(2.0_dp*real(nu, kind=dp) + 1.0_dp)
1717  expa(:) = exp(-z(:)**2)
1718  lambda = (ll - nu)/2
1719 
1720  IF (lambda > 0) THEN
1721  erfa = 0.0_dp
1722  DO i = 1, lambda
1723  erfa = erfa + (-prel)**i/real(i, kind=dp)*binomial_gen(lambda + nu + 0.5_dp, lambda - i)* &
1724  assoc_laguerre(z, real(nu, kind=dp) + 0.5_dp, i - 1)
1725  END DO
1726  erfa = erfa*expa*z**nu
1727 
1728  erfa = erfa + 2.0_dp*binomial_gen(lambda + nu + 0.5_dp, lambda)*znfn(z, expa, nu)
1729  ELSE
1730  erfa = 2.0_dp*znfn(z, expa, nu)
1731  END IF
1732 
1733  erfa = erfa*oab
1734 
1735  DEALLOCATE (expa, z)
1736  ELSE
1737  ! Contribution to potential is zero (properties of spherical harmonics)
1738  erfa = 0.0_dp
1739  END IF
1740 
1741  END SUBROUTINE potential_longrange_analytic
1742 
1743 ! **************************************************************************************************
1744 !> \brief Boys' function times z**n
1745 !> \param z ...
1746 !> \param expa ...
1747 !> \param n ...
1748 !> \return ...
1749 ! **************************************************************************************************
1750  ELEMENTAL FUNCTION znfn(z, expa, n) RESULT(res)
1751 
1752  REAL(kind=dp), INTENT(IN) :: z, expa
1753  INTEGER, INTENT(IN) :: n
1754  REAL(kind=dp) :: res
1755 
1756  INTEGER :: i
1757  REAL(kind=dp) :: z_exp
1758 
1759  IF (n >= 0) THEN
1760  IF (abs(z) < 1.0e-20) THEN
1761  res = z**n/(2.0_dp*real(n, kind=dp) + 1.0_dp)
1762  ELSE IF (n == 0) THEN
1763  res = rootpi/2.0_dp*erf(z)/z
1764  ELSE
1765  res = rootpi/4.0_dp*erf(z)/z**2 - expa/z/2.0_dp
1766  z_exp = -expa*0.5_dp
1767 
1768  DO i = 2, n
1769  res = (real(i, kind=dp) - 0.5_dp)*res/z + z_exp
1770  z_exp = z_exp*z
1771  END DO
1772  END IF
1773  ELSE ! Set it to zero (no Boys' function, to keep the ELEMENTAL keyword)
1774  res = 0.0_dp
1775  END IF
1776 
1777  END FUNCTION znfn
1778 
1779 ! **************************************************************************************************
1780 !> \brief ...
1781 !> \param z ...
1782 !> \param a ...
1783 !> \param n ...
1784 !> \return ...
1785 ! **************************************************************************************************
1786  ELEMENTAL FUNCTION assoc_laguerre(z, a, n) RESULT(res)
1787 
1788  REAL(kind=dp), INTENT(IN) :: z, a
1789  INTEGER, INTENT(IN) :: n
1790  REAL(kind=dp) :: res
1791 
1792  INTEGER :: i
1793  REAL(kind=dp) :: f0, f1
1794 
1795  IF (n == 0) THEN
1796  res = 1.0_dp
1797  ELSE IF (n == 1) THEN
1798  res = a + 1.0_dp - z
1799  ELSE IF (n > 0) THEN
1800  f0 = 1.0_dp
1801  f1 = a + 1.0_dp - z
1802 
1803  DO i = 3, n
1804  res = (2.0_dp + (a - 1.0_dp - z)/real(i, dp))*f1 - (1.0_dp + (a - 1.0_dp)/real(i, dp))*f0
1805  f0 = f1
1806  f1 = res
1807  END DO
1808  ELSE ! n is negative, set it zero (no polynomials, to keep the ELEMENTAL keyword)
1809  res = 0.0_dp
1810  END IF
1811 
1812  END FUNCTION assoc_laguerre
1813 
1814 ! **************************************************************************************************
1815 !> \brief Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
1816 !> \param opmat operator matrix (e.g. Kohn-Sham matrix or overlap matrix)
1817 !> \param pmat density matrix
1818 !> \return value of trace
1819 !> \par History
1820 !> * 08.2008 created [Juerg Hutter]
1821 ! **************************************************************************************************
1822  PURE FUNCTION atom_trace(opmat, pmat) RESULT(trace)
1823  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: opmat, pmat
1824  REAL(kind=dp) :: trace
1825 
1826  trace = accurate_dot_product(opmat, pmat)
1827 
1828  END FUNCTION atom_trace
1829 
1830 ! **************************************************************************************************
1831 !> \brief Calculate a potential matrix by integrating the potential on an atomic radial grid.
1832 !> \param imat potential matrix
1833 !> \param cpot potential on the atomic radial grid
1834 !> \param basis atomic basis set
1835 !> \param derivatives order of derivatives
1836 !> \par History
1837 !> * 08.2008 created [Juerg Hutter]
1838 ! **************************************************************************************************
1839  SUBROUTINE numpot_matrix(imat, cpot, basis, derivatives)
1840  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: imat
1841  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: cpot
1842  TYPE(atom_basis_type), INTENT(INOUT) :: basis
1843  INTEGER, INTENT(IN) :: derivatives
1844 
1845  INTEGER :: i, j, l, n
1846 
1847  SELECT CASE (derivatives)
1848  CASE (0)
1849  DO l = 0, lmat
1850  n = basis%nbas(l)
1851  DO i = 1, n
1852  DO j = i, n
1853  imat(i, j, l) = imat(i, j, l) + &
1854  integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
1855  imat(j, i, l) = imat(i, j, l)
1856  END DO
1857  END DO
1858  END DO
1859  CASE (1)
1860  DO l = 0, lmat
1861  n = basis%nbas(l)
1862  DO i = 1, n
1863  DO j = i, n
1864  imat(i, j, l) = imat(i, j, l) + &
1865  integrate_grid(cpot, basis%dbf(:, i, l), basis%bf(:, j, l), basis%grid)
1866  imat(i, j, l) = imat(i, j, l) + &
1867  integrate_grid(cpot, basis%bf(:, i, l), basis%dbf(:, j, l), basis%grid)
1868  imat(j, i, l) = imat(i, j, l)
1869  END DO
1870  END DO
1871  END DO
1872  CASE (2)
1873  DO l = 0, lmat
1874  n = basis%nbas(l)
1875  DO i = 1, n
1876  DO j = i, n
1877  imat(i, j, l) = imat(i, j, l) + &
1878  integrate_grid(cpot, basis%dbf(:, i, l), basis%dbf(:, j, l), basis%grid)
1879  imat(j, i, l) = imat(i, j, l)
1880  END DO
1881  END DO
1882  END DO
1883  CASE DEFAULT
1884  cpabort("")
1885  END SELECT
1886 
1887  END SUBROUTINE numpot_matrix
1888 
1889 ! **************************************************************************************************
1890 !> \brief Contract Coulomb Electron Repulsion Integrals.
1891 !> \param jmat ...
1892 !> \param erint ...
1893 !> \param pmat ...
1894 !> \param nsize ...
1895 !> \param all_nu ...
1896 !> \par History
1897 !> * 08.2008 created [Juerg Hutter]
1898 ! **************************************************************************************************
1899  PURE SUBROUTINE ceri_contract(jmat, erint, pmat, nsize, all_nu)
1900  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: jmat
1901  TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1902  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1903  INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1904  LOGICAL, INTENT(IN), OPTIONAL :: all_nu
1905 
1906  INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, ll, &
1907  n1, n2
1908  LOGICAL :: have_all_nu
1909  REAL(kind=dp) :: eint, f1, f2
1910 
1911  IF (PRESENT(all_nu)) THEN
1912  have_all_nu = all_nu
1913  ELSE
1914  have_all_nu = .false.
1915  END IF
1916 
1917  jmat(:, :, :) = 0._dp
1918  ll = 0
1919  DO l1 = 0, lmat
1920  n1 = nsize(l1)
1921  DO l2 = 0, l1
1922  n2 = nsize(l2)
1923  ll = ll + 1
1924  ij1 = 0
1925  DO i1 = 1, n1
1926  DO j1 = i1, n1
1927  ij1 = ij1 + 1
1928  f1 = 1._dp
1929  IF (i1 /= j1) f1 = 2._dp
1930  ij2 = 0
1931  DO i2 = 1, n2
1932  DO j2 = i2, n2
1933  ij2 = ij2 + 1
1934  f2 = 1._dp
1935  IF (i2 /= j2) f2 = 2._dp
1936  eint = erint(ll)%int(ij1, ij2)
1937  IF (l1 == l2) THEN
1938  jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1939  ELSE
1940  jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1941  jmat(i2, j2, l2) = jmat(i2, j2, l2) + f1*pmat(i1, j1, l1)*eint
1942  END IF
1943  END DO
1944  END DO
1945  END DO
1946  END DO
1947  IF (have_all_nu) THEN
1948  ! skip integral blocks with nu/=0
1949  ll = ll + l2
1950  END IF
1951  END DO
1952  END DO
1953  DO l1 = 0, lmat
1954  n1 = nsize(l1)
1955  DO i1 = 1, n1
1956  DO j1 = i1, n1
1957  jmat(j1, i1, l1) = jmat(i1, j1, l1)
1958  END DO
1959  END DO
1960  END DO
1961 
1962  END SUBROUTINE ceri_contract
1963 
1964 ! **************************************************************************************************
1965 !> \brief Contract exchange Electron Repulsion Integrals.
1966 !> \param kmat ...
1967 !> \param erint ...
1968 !> \param pmat ...
1969 !> \param nsize ...
1970 !> \par History
1971 !> * 08.2008 created [Juerg Hutter]
1972 ! **************************************************************************************************
1973  PURE SUBROUTINE eeri_contract(kmat, erint, pmat, nsize)
1974  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1975  TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1976  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1977  INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1978 
1979  INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, lh, &
1980  ll, n1, n2, nu
1981  REAL(kind=dp) :: almn, eint, f1, f2
1982  REAL(kind=dp), DIMENSION(0:maxfac) :: arho
1983 
1984  arho = 0._dp
1985  DO ll = 0, maxfac, 2
1986  lh = ll/2
1987  arho(ll) = fac(ll)/fac(lh)**2
1988  END DO
1989 
1990  kmat(:, :, :) = 0._dp
1991  ll = 0
1992  DO l1 = 0, lmat
1993  n1 = nsize(l1)
1994  DO l2 = 0, l1
1995  n2 = nsize(l2)
1996  DO nu = abs(l1 - l2), l1 + l2, 2
1997  almn = arho(-l1 + l2 + nu)*arho(l1 - l2 + nu)*arho(l1 + l2 - nu)/(real(l1 + l2 + nu + 1, dp)*arho(l1 + l2 + nu))
1998  almn = -0.5_dp*almn
1999  ll = ll + 1
2000  ij1 = 0
2001  DO i1 = 1, n1
2002  DO j1 = i1, n1
2003  ij1 = ij1 + 1
2004  f1 = 1._dp
2005  IF (i1 /= j1) f1 = 2._dp
2006  ij2 = 0
2007  DO i2 = 1, n2
2008  DO j2 = i2, n2
2009  ij2 = ij2 + 1
2010  f2 = 1._dp
2011  IF (i2 /= j2) f2 = 2._dp
2012  eint = erint(ll)%int(ij1, ij2)
2013  IF (l1 == l2) THEN
2014  kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2015  ELSE
2016  kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2017  kmat(i2, j2, l2) = kmat(i2, j2, l2) + f1*almn*pmat(i1, j1, l1)*eint
2018  END IF
2019  END DO
2020  END DO
2021  END DO
2022  END DO
2023  END DO
2024  END DO
2025  END DO
2026  DO l1 = 0, lmat
2027  n1 = nsize(l1)
2028  DO i1 = 1, n1
2029  DO j1 = i1, n1
2030  kmat(j1, i1, l1) = kmat(i1, j1, l1)
2031  END DO
2032  END DO
2033  END DO
2034 
2035  END SUBROUTINE eeri_contract
2036 
2037 ! **************************************************************************************************
2038 !> \brief Calculate the error matrix for each angular momentum.
2039 !> \param emat error matrix
2040 !> \param demax the largest absolute value of error matrix elements
2041 !> \param kmat Kohn-Sham matrix
2042 !> \param pmat electron density matrix
2043 !> \param umat transformation matrix which reduce overlap matrix to its unitary form
2044 !> \param upmat transformation matrix which reduce density matrix to its unitary form
2045 !> \param nval number of linear-independent contracted basis functions
2046 !> \param nbs number of contracted basis functions
2047 !> \par History
2048 !> * 08.2008 created [Juerg Hutter]
2049 ! **************************************************************************************************
2050  PURE SUBROUTINE err_matrix(emat, demax, kmat, pmat, umat, upmat, nval, nbs)
2051  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: emat
2052  REAL(kind=dp), INTENT(OUT) :: demax
2053  REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN) :: kmat, pmat, umat, upmat
2054  INTEGER, DIMENSION(0:), INTENT(IN) :: nval, nbs
2055 
2056  INTEGER :: l, m, n
2057  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tkmat, tpmat
2058 
2059  emat = 0._dp
2060  DO l = 0, lmat
2061  n = nval(l)
2062  m = nbs(l)
2063  IF (m > 0) THEN
2064  ALLOCATE (tkmat(1:m, 1:m), tpmat(1:m, 1:m))
2065  tkmat = 0._dp
2066  tpmat = 0._dp
2067  tkmat(1:m, 1:m) = matmul(transpose(umat(1:n, 1:m, l)), matmul(kmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
2068  tpmat(1:m, 1:m) = matmul(transpose(umat(1:n, 1:m, l)), matmul(pmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
2069  tpmat(1:m, 1:m) = matmul(upmat(1:m, 1:m, l), matmul(tpmat(1:m, 1:m), upmat(1:m, 1:m, l)))
2070 
2071  emat(1:m, 1:m, l) = matmul(tkmat(1:m, 1:m), tpmat(1:m, 1:m)) - matmul(tpmat(1:m, 1:m), tkmat(1:m, 1:m))
2072 
2073  DEALLOCATE (tkmat, tpmat)
2074  END IF
2075  END DO
2076  demax = maxval(abs(emat))
2077 
2078  END SUBROUTINE err_matrix
2079 
2080 ! **************************************************************************************************
2081 !> \brief Calculate Slater density on a radial grid.
2082 !> \param density1 alpha-spin electron density
2083 !> \param density2 beta-spin electron density
2084 !> \param zcore nuclear charge
2085 !> \param state electronic state
2086 !> \param grid atomic radial grid
2087 !> \par History
2088 !> * 06.2018 bugfix [Rustam Khaliullin]
2089 !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
2090 !> * 12.2008 created [Juerg Hutter]
2091 !> \note An initial electron density to guess atomic orbitals.
2092 ! **************************************************************************************************
2093  SUBROUTINE slater_density(density1, density2, zcore, state, grid)
2094  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: density1, density2
2095  INTEGER, INTENT(IN) :: zcore
2096  TYPE(atom_state), INTENT(IN) :: state
2097  TYPE(grid_atom_type), INTENT(IN) :: grid
2098 
2099  INTEGER :: counter, i, l, mc, mm(0:lmat), mo, n, ns
2100  INTEGER, DIMENSION(lmat+1, 20) :: ne
2101  REAL(kind=dp) :: a, pf
2102 
2103  ! fill out a table with the total number of electrons
2104  ! core + valence. format ne(l,n)
2105  ns = SIZE(ne, 2)
2106  ne = 0
2107  mm = get_maxn_occ(state%core)
2108  DO l = 0, lmat
2109  mo = state%maxn_occ(l)
2110  IF (sum(state%core(l, :)) == 0) THEN
2111  cpassert(ns >= l + mo)
2112  DO counter = 1, mo
2113  ne(l + 1, l + counter) = nint(state%occ(l, counter))
2114  END DO
2115  ELSE
2116  mc = mm(l) ! number of levels in the core
2117  cpassert(sum(state%occ(l, 1:mc)) == 0)
2118  cpassert(ns >= l + mc)
2119  DO counter = 1, mc
2120  ne(l + 1, l + counter) = nint(state%core(l, counter))
2121  END DO
2122  cpassert(ns >= l + mc + mo)
2123  DO counter = mc + 1, mc + mo
2124  ne(l + 1, l + counter) = nint(state%occ(l, counter))
2125  END DO
2126  END IF
2127  END DO
2128 
2129  density1 = 0._dp
2130  density2 = 0._dp
2131  DO l = 0, state%maxl_occ
2132  DO i = 1, SIZE(state%occ, 2)
2133  IF (state%occ(l, i) > 0._dp) THEN
2134  n = i + l
2135  a = srules(zcore, ne, n, l)
2136  pf = 1._dp/sqrt(fac(2*n))*(2._dp*a)**(n + 0.5_dp)
2137  IF (state%multiplicity == -1) THEN
2138  density1(:) = density1(:) + state%occ(l, i)/fourpi*(grid%rad(:)**(n - 1)*exp(-a*grid%rad(:))*pf)**2
2139  ELSE
2140  density1(:) = density1(:) + state%occa(l, i)/fourpi*(grid%rad(:)**(n - 1)*exp(-a*grid%rad(:))*pf)**2
2141  density2(:) = density2(:) + state%occb(l, i)/fourpi*(grid%rad(:)**(n - 1)*exp(-a*grid%rad(:))*pf)**2
2142  END IF
2143  END IF
2144  END DO
2145  END DO
2146 
2147  END SUBROUTINE slater_density
2148 
2149 ! **************************************************************************************************
2150 !> \brief Calculate the functional derivative of the Wigner (correlation) - Slater (exchange)
2151 !> density functional.
2152 !> \param rho electron density on a radial grid
2153 !> \param vxc (output) exchange-correlation potential
2154 !> \par History
2155 !> * 12.2008 created [Juerg Hutter]
2156 !> \note A model XC-potential to guess atomic orbitals.
2157 ! **************************************************************************************************
2158  PURE SUBROUTINE wigner_slater_functional(rho, vxc)
2159  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rho
2160  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: vxc
2161 
2162  INTEGER :: i
2163  REAL(kind=dp) :: ec, ex, rs, vc, vx
2164 
2165  vxc = 0._dp
2166  DO i = 1, SIZE(rho)
2167  IF (rho(i) > 1.e-20_dp) THEN
2168  ! 3/4 * (3/pi)^{1/3} == 0.7385588
2169  ex = -0.7385588_dp*rho(i)**0.333333333_dp
2170  vx = 1.333333333_dp*ex
2171  rs = (3._dp/fourpi/rho(i))**0.333333333_dp
2172  ec = -0.88_dp/(rs + 7.8_dp)
2173  vc = ec*(1._dp + rs/(3._dp*(rs + 7.8_dp)))
2174  vxc(i) = vx + vc
2175  END IF
2176  END DO
2177 
2178  END SUBROUTINE wigner_slater_functional
2179 
2180 ! **************************************************************************************************
2181 !> \brief Check that the atomic multiplicity is consistent with the electronic structure method.
2182 !> \param method electronic structure method
2183 !> \param multiplicity multiplicity
2184 !> \return consistency status
2185 !> \par History
2186 !> * 11.2009 unrestricted KS and HF methods [Juerg Hutter]
2187 ! **************************************************************************************************
2188  PURE FUNCTION atom_consistent_method(method, multiplicity) RESULT(consistent)
2189  INTEGER, INTENT(IN) :: method, multiplicity
2190  LOGICAL :: consistent
2191 
2192  ! multiplicity == -1 means it has not been specified explicitly;
2193  ! see the source code of the subroutine atom_set_occupation() for further details.
2194  SELECT CASE (method)
2195  CASE DEFAULT
2196  consistent = .false.
2197  CASE (do_rks_atom)
2198  consistent = (multiplicity == -1)
2199  CASE (do_uks_atom)
2200  consistent = (multiplicity /= -1)
2201  CASE (do_rhf_atom)
2202  consistent = (multiplicity == -1)
2203  CASE (do_uhf_atom)
2204  consistent = (multiplicity /= -1)
2205  CASE (do_rohf_atom)
2206  consistent = .false.
2207  END SELECT
2208 
2209  END FUNCTION atom_consistent_method
2210 
2211 ! **************************************************************************************************
2212 !> \brief Calculate the total electron density at R=0.
2213 !> \param atom information about the atomic kind
2214 !> \param rho0 (output) computed density
2215 !> \par History
2216 !> * 05.2016 created [Juerg Hutter]
2217 ! **************************************************************************************************
2218  SUBROUTINE get_rho0(atom, rho0)
2219  TYPE(atom_type), INTENT(IN) :: atom
2220  REAL(kind=dp), INTENT(OUT) :: rho0
2221 
2222  INTEGER :: m0, m1, m2, method, nr
2223  LOGICAL :: nlcc, spinpol
2224  REAL(kind=dp) :: d0, d1, d2, r0, r1, r2, w0, w1
2225  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: xfun
2226  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rho
2227 
2228  method = atom%method_type
2229  SELECT CASE (method)
2230  CASE (do_rks_atom, do_rhf_atom)
2231  spinpol = .false.
2232  CASE (do_uks_atom, do_uhf_atom)
2233  spinpol = .true.
2234  CASE (do_rohf_atom)
2235  cpabort("")
2236  CASE DEFAULT
2237  cpabort("")
2238  END SELECT
2239 
2240  nr = atom%basis%grid%nr
2241  nlcc = .false.
2242  IF (atom%potential%ppot_type == gth_pseudo) THEN
2243  nlcc = atom%potential%gth_pot%nlcc
2244  ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
2245  nlcc = atom%potential%upf_pot%core_correction
2246  ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
2247  nlcc = atom%potential%sgp_pot%has_nlcc
2248  END IF
2249  IF (nlcc) THEN
2250  ALLOCATE (xfun(nr))
2251  END IF
2252 
2253  m0 = maxval(minloc(atom%basis%grid%rad))
2254  IF (m0 == nr) THEN
2255  m1 = m0 - 1
2256  m2 = m0 - 2
2257  ELSE IF (m0 == 1) THEN
2258  m1 = 2
2259  m2 = 3
2260  ELSE
2261  cpabort("GRID Definition incompatible")
2262  END IF
2263  r0 = atom%basis%grid%rad(m0)
2264  r1 = atom%basis%grid%rad(m1)
2265  r2 = atom%basis%grid%rad(m2)
2266  w0 = r1/(r1 - r0)
2267  w1 = 1 - w0
2268 
2269  IF (spinpol) THEN
2270  ALLOCATE (rho(nr, 2))
2271  CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
2272  CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
2273  IF (nlcc) THEN
2274  xfun = 0.0_dp
2275  CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2276  rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
2277  rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
2278  END IF
2279  rho(:, 1) = rho(:, 1) + rho(:, 2)
2280  ELSE
2281  ALLOCATE (rho(nr, 1))
2282  CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
2283  IF (nlcc) THEN
2284  CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2285  END IF
2286  END IF
2287  d0 = rho(m0, 1)
2288  d1 = rho(m1, 1)
2289  d2 = rho(m2, 1)
2290 
2291  rho0 = w0*d0 + w1*d1
2292  rho0 = max(rho0, 0.0_dp)
2293 
2294  DEALLOCATE (rho)
2295  IF (nlcc) THEN
2296  DEALLOCATE (xfun)
2297  END IF
2298 
2299  END SUBROUTINE get_rho0
2300 
2301 ! **************************************************************************************************
2302 !> \brief Print condition numbers of the given atomic basis set.
2303 !> \param basis atomic basis set
2304 !> \param crad atomic covalent radius
2305 !> \param iw output file unit
2306 !> \par History
2307 !> * 05.2016 created [Juerg Hutter]
2308 ! **************************************************************************************************
2309  SUBROUTINE atom_condnumber(basis, crad, iw)
2310  TYPE(atom_basis_type), POINTER :: basis
2311  REAL(kind=dp) :: crad
2312  INTEGER, INTENT(IN) :: iw
2313 
2314  INTEGER :: i
2315  REAL(kind=dp) :: ci
2316  REAL(kind=dp), DIMENSION(10) :: cnum, rad
2317 
2318  WRITE (iw, '(/,A,F8.4)') " Basis Set Condition Numbers: 2*covalent radius=", 2*crad
2319  CALL init_orbital_pointers(lmat)
2320  CALL init_spherical_harmonics(lmat, 0)
2321  cnum = 0.0_dp
2322  DO i = 1, 9
2323  ci = 2.0_dp*(0.85_dp + i*0.05_dp)
2324  rad(i) = crad*ci
2325  CALL atom_basis_condnum(basis, rad(i), cnum(i))
2326  WRITE (iw, '(A,F15.3,T50,A,F14.4)') " Lattice constant:", &
2327  rad(i), "Condition number:", cnum(i)
2328  END DO
2329  rad(10) = 0.01_dp
2330  CALL atom_basis_condnum(basis, rad(10), cnum(10))
2331  WRITE (iw, '(A,A,T50,A,F14.4)') " Lattice constant:", &
2332  " Inf", "Condition number:", cnum(i)
2333  CALL deallocate_orbital_pointers
2334  CALL deallocate_spherical_harmonics
2335 
2336  END SUBROUTINE atom_condnumber
2337 
2338 ! **************************************************************************************************
2339 !> \brief Estimate completeness of the given atomic basis set.
2340 !> \param basis atomic basis set
2341 !> \param zv atomic number
2342 !> \param iw output file unit
2343 ! **************************************************************************************************
2344  SUBROUTINE atom_completeness(basis, zv, iw)
2345  TYPE(atom_basis_type), POINTER :: basis
2346  INTEGER, INTENT(IN) :: zv, iw
2347 
2348  INTEGER :: i, j, l, ll, m, n, nbas, nl, nr
2349  INTEGER, DIMENSION(0:lmat) :: nelem, nlmax, nlmin
2350  INTEGER, DIMENSION(lmat+1, 7) :: ne
2351  REAL(kind=dp) :: al, c1, c2, pf
2352  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sfun
2353  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: bmat
2354  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: omat
2355  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sint
2356  REAL(kind=dp), DIMENSION(0:lmat, 10) :: snl
2357  REAL(kind=dp), DIMENSION(2) :: sse
2358  REAL(kind=dp), DIMENSION(lmat+1, 7) :: sexp
2359 
2360  ne = 0
2361  nelem = 0
2362  nelem(0:3) = ptable(zv)%e_conv(0:3)
2363  DO l = 0, lmat
2364  ll = 2*(2*l + 1)
2365  DO i = 1, 7
2366  IF (nelem(l) >= ll) THEN
2367  ne(l + 1, i) = ll
2368  nelem(l) = nelem(l) - ll
2369  ELSE IF (nelem(l) > 0) THEN
2370  ne(l + 1, i) = nelem(l)
2371  nelem(l) = 0
2372  ELSE
2373  EXIT
2374  END IF
2375  END DO
2376  END DO
2377 
2378  nlmin = 1
2379  nlmax = 1
2380  DO l = 0, lmat
2381  nlmin(l) = l + 1
2382  DO i = 1, 7
2383  IF (ne(l + 1, i) > 0) THEN
2384  nlmax(l) = i + l
2385  END IF
2386  END DO
2387  nlmax(l) = max(nlmax(l), nlmin(l) + 1)
2388  END DO
2389 
2390  ! Slater exponents
2391  sexp = 0.0_dp
2392  DO l = 0, lmat
2393  sse(1) = 0.05_dp
2394  sse(2) = 10.0_dp
2395  DO i = l + 1, 7
2396  sexp(l + 1, i) = srules(zv, ne, i, l)
2397  IF (ne(l + 1, i - l) > 0) THEN
2398  sse(1) = max(sse(1), sexp(l + 1, i))
2399  sse(2) = min(sse(2), sexp(l + 1, i))
2400  END IF
2401  END DO
2402  DO i = 1, 10
2403  snl(l, i) = abs(2._dp*sse(1) - 0.5_dp*sse(2))/9._dp*real(i - 1, kind=dp) + 0.5_dp*min(sse(1), sse(2))
2404  END DO
2405  END DO
2406 
2407  nbas = maxval(basis%nbas)
2408  ALLOCATE (omat(nbas, nbas, 0:lmat))
2409  nr = SIZE(basis%bf, 1)
2410  ALLOCATE (sfun(nr), sint(10, 2, nbas, 0:lmat))
2411  sint = 0._dp
2412 
2413  ! calculate overlaps between test functions and basis
2414  DO l = 0, lmat
2415  DO i = 1, 10
2416  al = snl(l, i)
2417  nl = nlmin(l)
2418  pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
2419  sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*exp(-al*basis%grid%rad(1:nr))
2420  DO j = 1, basis%nbas(l)
2421  sint(i, 1, j, l) = sum(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2422  END DO
2423  nl = nlmax(l)
2424  pf = (2._dp*al)**nl*sqrt(2._dp*al/fac(2*nl))
2425  sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*exp(-al*basis%grid%rad(1:nr))
2426  DO j = 1, basis%nbas(l)
2427  sint(i, 2, j, l) = sum(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2428  END DO
2429  END DO
2430  END DO
2431 
2432  DO l = 0, lmat
2433  n = basis%nbas(l)
2434  IF (n <= 0) cycle
2435  m = basis%nprim(l)
2436  SELECT CASE (basis%basis_type)
2437  CASE DEFAULT
2438  cpabort("")
2439  CASE (gto_basis)
2440  CALL sg_overlap(omat(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
2441  CASE (cgto_basis)
2442  ALLOCATE (bmat(m, m))
2443  CALL sg_overlap(bmat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
2444  CALL contract2(omat(1:n, 1:n, l), bmat(1:m, 1:m), basis%cm(1:m, 1:n, l))
2445  DEALLOCATE (bmat)
2446  CASE (sto_basis)
2447  CALL sto_overlap(omat(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
2448  basis%ns(1:n, l), basis%as(1:n, l))
2449  CASE (num_basis)
2450  cpabort("")
2451  END SELECT
2452  CALL invmat_symm(omat(1:n, 1:n, l))
2453  END DO
2454 
2455  WRITE (iw, '(/,A)') " Basis Set Completeness Estimates"
2456  DO l = 0, lmat
2457  n = basis%nbas(l)
2458  IF (n <= 0) cycle
2459  WRITE (iw, '(A,I3)') " L-quantum number: ", l
2460  WRITE (iw, '(A,T31,A,I2,T61,A,I2)') " Slater Exponent", "Completeness n-qm=", nlmin(l), &
2461  "Completeness n-qm=", nlmax(l)
2462  DO i = 10, 1, -1
2463  c1 = dot_product(sint(i, 1, 1:n, l), matmul(omat(1:n, 1:n, l), sint(i, 1, 1:n, l)))
2464  c2 = dot_product(sint(i, 2, 1:n, l), matmul(omat(1:n, 1:n, l), sint(i, 2, 1:n, l)))
2465  WRITE (iw, "(T6,F14.6,T41,F10.6,T71,F10.6)") snl(l, i), c1, c2
2466  END DO
2467  END DO
2468 
2469  DEALLOCATE (omat, sfun, sint)
2470 
2471  END SUBROUTINE atom_completeness
2472 
2473 ! **************************************************************************************************
2474 !> \brief Calculate the condition number of the given atomic basis set.
2475 !> \param basis atomic basis set
2476 !> \param rad lattice constant (e.g. doubled atomic covalent radius)
2477 !> \param cnum (output) condition number
2478 ! **************************************************************************************************
2479  SUBROUTINE atom_basis_condnum(basis, rad, cnum)
2480  TYPE(atom_basis_type) :: basis
2481  REAL(kind=dp), INTENT(IN) :: rad
2482  REAL(kind=dp), INTENT(OUT) :: cnum
2483 
2484  INTEGER :: ia, ib, imax, info, ix, iy, iz, ja, jb, &
2485  ka, kb, l, la, lb, lwork, na, nb, &
2486  nbas, nna, nnb
2487  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ibptr
2488  REAL(kind=dp) :: r1, r2, reps, rmax
2489  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weig, work
2490  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: smat
2491  REAL(kind=dp), DIMENSION(2*lmat+1, 2*lmat+1) :: sab
2492  REAL(kind=dp), DIMENSION(3) :: rab
2493  REAL(kind=dp), DIMENSION(:), POINTER :: zeta, zetb
2494 
2495  ! Number of spherical Gaussian orbitals with angular momentum lmat: nso(lmat) = 2*lmat+1
2496 
2497  ! total number of basis functions
2498  nbas = 0
2499  DO l = 0, lmat
2500  nbas = nbas + basis%nbas(l)*(2*l + 1)
2501  END DO
2502 
2503  ALLOCATE (smat(nbas, nbas), ibptr(nbas, 0:lmat))
2504  smat = 0.0_dp
2505  ibptr = 0
2506  na = 0
2507  DO l = 0, lmat
2508  DO ia = 1, basis%nbas(l)
2509  ibptr(ia, l) = na + 1
2510  na = na + (2*l + 1)
2511  END DO
2512  END DO
2513 
2514  reps = 1.e-14_dp
2515  IF (basis%basis_type == gto_basis .OR. &
2516  basis%basis_type == cgto_basis) THEN
2517  DO la = 0, lmat
2518  na = basis%nprim(la)
2519  nna = 2*la + 1
2520  IF (na == 0) cycle
2521  zeta => basis%am(:, la)
2522  DO lb = 0, lmat
2523  nb = basis%nprim(lb)
2524  nnb = 2*lb + 1
2525  IF (nb == 0) cycle
2526  zetb => basis%am(:, lb)
2527  DO ia = 1, na
2528  DO ib = 1, nb
2529  IF (rad < 0.1_dp) THEN
2530  imax = 0
2531  ELSE
2532  r1 = exp_radius(la, zeta(ia), reps, 1.0_dp)
2533  r2 = exp_radius(lb, zetb(ib), reps, 1.0_dp)
2534  rmax = max(2._dp*r1, 2._dp*r2)
2535  imax = int(rmax/rad) + 1
2536  END IF
2537  IF (imax > 1) THEN
2538  CALL overlap_ab_sp(la, zeta(ia), lb, zetb(ib), rad, sab)
2539  IF (basis%basis_type == gto_basis) THEN
2540  ja = ibptr(ia, la)
2541  jb = ibptr(ib, lb)
2542  smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + sab(1:nna, 1:nnb)
2543  ELSEIF (basis%basis_type == cgto_basis) THEN
2544  DO ka = 1, basis%nbas(la)
2545  DO kb = 1, basis%nbas(lb)
2546  ja = ibptr(ka, la)
2547  jb = ibptr(kb, lb)
2548  smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2549  sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2550  END DO
2551  END DO
2552  END IF
2553  ELSE
2554  DO ix = -imax, imax
2555  rab(1) = rad*ix
2556  DO iy = -imax, imax
2557  rab(2) = rad*iy
2558  DO iz = -imax, imax
2559  rab(3) = rad*iz
2560  CALL overlap_ab_s(la, zeta(ia), lb, zetb(ib), rab, sab)
2561  IF (basis%basis_type == gto_basis) THEN
2562  ja = ibptr(ia, la)
2563  jb = ibptr(ib, lb)
2564  smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) &
2565  + sab(1:nna, 1:nnb)
2566  ELSEIF (basis%basis_type == cgto_basis) THEN
2567  DO ka = 1, basis%nbas(la)
2568  DO kb = 1, basis%nbas(lb)
2569  ja = ibptr(ka, la)
2570  jb = ibptr(kb, lb)
2571  smat(ja:ja + nna - 1, jb:jb + nnb - 1) = &
2572  smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2573  sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2574  END DO
2575  END DO
2576  END IF
2577  END DO
2578  END DO
2579  END DO
2580  END IF
2581  END DO
2582  END DO
2583  END DO
2584  END DO
2585  ELSE
2586  cpabort("Condition number not available for this basis type")
2587  END IF
2588 
2589  info = 0
2590  lwork = max(nbas*nbas, nbas + 100)
2591  ALLOCATE (weig(nbas), work(lwork))
2592 
2593  CALL lapack_ssyev("N", "U", nbas, smat, nbas, weig, work, lwork, info)
2594  cpassert(info == 0)
2595  IF (weig(1) < 0.0_dp) THEN
2596  cnum = 100._dp
2597  ELSE
2598  cnum = abs(weig(nbas)/weig(1))
2599  cnum = log10(cnum)
2600  END IF
2601 
2602  DEALLOCATE (smat, ibptr, weig, work)
2603 
2604  END SUBROUTINE atom_basis_condnum
2605 
2606 ! **************************************************************************************************
2607 !> \brief Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
2608 !> \param int (output) contracted matrix
2609 !> \param omat uncontracted matrix
2610 !> \param cm matrix of contraction coefficients
2611 ! **************************************************************************************************
2612  SUBROUTINE contract2(int, omat, cm)
2613  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2614  REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2615 
2616  CHARACTER(len=*), PARAMETER :: routinen = 'contract2'
2617 
2618  INTEGER :: handle, m, n
2619 
2620  CALL timeset(routinen, handle)
2621 
2622  n = SIZE(int, 1)
2623  m = SIZE(omat, 1)
2624 
2625  int(1:n, 1:n) = matmul(transpose(cm(1:m, 1:n)), matmul(omat(1:m, 1:m), cm(1:m, 1:n)))
2626 
2627  CALL timestop(handle)
2628 
2629  END SUBROUTINE contract2
2630 
2631 ! **************************************************************************************************
2632 !> \brief Same as contract2(), but add the new contracted matrix to the old one
2633 !> instead of overwriting it.
2634 !> \param int (input/output) contracted matrix to be updated
2635 !> \param omat uncontracted matrix
2636 !> \param cm matrix of contraction coefficients
2637 ! **************************************************************************************************
2638  SUBROUTINE contract2add(int, omat, cm)
2639  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2640  REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2641 
2642  CHARACTER(len=*), PARAMETER :: routinen = 'contract2add'
2643 
2644  INTEGER :: handle, m, n
2645 
2646  CALL timeset(routinen, handle)
2647 
2648  n = SIZE(int, 1)
2649  m = SIZE(omat, 1)
2650 
2651  int(1:n, 1:n) = int(1:n, 1:n) + matmul(transpose(cm(1:m, 1:n)), matmul(omat(1:m, 1:m), cm(1:m, 1:n)))
2652 
2653  CALL timestop(handle)
2654 
2655  END SUBROUTINE contract2add
2656 
2657 ! **************************************************************************************************
2658 !> \brief Contract a matrix of Electron Repulsion Integrals (ERI-s).
2659 !> \param eri (output) contracted ERI-s
2660 !> \param omat uncontracted ERI-s
2661 !> \param cm1 first matrix of contraction coefficients
2662 !> \param cm2 second matrix of contraction coefficients
2663 ! **************************************************************************************************
2664  SUBROUTINE contract4(eri, omat, cm1, cm2)
2665  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: eri
2666  REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm1, cm2
2667 
2668  CHARACTER(len=*), PARAMETER :: routinen = 'contract4'
2669 
2670  INTEGER :: handle, i1, i2, m1, m2, mm1, mm2, n1, &
2671  n2, nn1, nn2
2672  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: amat, atran, bmat, btran, hint
2673 
2674  CALL timeset(routinen, handle)
2675 
2676  m1 = SIZE(cm1, 1)
2677  n1 = SIZE(cm1, 2)
2678  m2 = SIZE(cm2, 1)
2679  n2 = SIZE(cm2, 2)
2680  nn1 = SIZE(eri, 1)
2681  nn2 = SIZE(eri, 2)
2682  mm1 = SIZE(omat, 1)
2683  mm2 = SIZE(omat, 2)
2684 
2685  ALLOCATE (amat(m1, m1), atran(n1, n1), bmat(m2, m2), btran(n2, n2))
2686  ALLOCATE (hint(mm1, nn2))
2687 
2688  DO i1 = 1, mm1
2689  CALL iunpack(bmat(1:m2, 1:m2), omat(i1, 1:mm2), m2)
2690  CALL contract2(btran(1:n2, 1:n2), bmat(1:m2, 1:m2), cm2(1:m2, 1:n2))
2691  CALL ipack(btran(1:n2, 1:n2), hint(i1, 1:nn2), n2)
2692  END DO
2693 
2694  DO i2 = 1, nn2
2695  CALL iunpack(amat(1:m1, 1:m1), hint(1:mm1, i2), m1)
2696  CALL contract2(atran(1:n1, 1:n1), amat(1:m1, 1:m1), cm1(1:m1, 1:n1))
2697  CALL ipack(atran(1:n1, 1:n1), eri(1:nn1, i2), n1)
2698  END DO
2699 
2700  DEALLOCATE (amat, atran, bmat, btran)
2701  DEALLOCATE (hint)
2702 
2703  CALL timestop(handle)
2704 
2705  END SUBROUTINE contract4
2706 
2707 ! **************************************************************************************************
2708 !> \brief Pack an n x n square real matrix into a 1-D vector.
2709 !> \param mat matrix to pack
2710 !> \param vec resulting vector
2711 !> \param n size of the matrix
2712 ! **************************************************************************************************
2713  PURE SUBROUTINE ipack(mat, vec, n)
2714  REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat
2715  REAL(dp), DIMENSION(:), INTENT(INOUT) :: vec
2716  INTEGER, INTENT(IN) :: n
2717 
2718  INTEGER :: i, ij, j
2719 
2720  ij = 0
2721  DO i = 1, n
2722  DO j = i, n
2723  ij = ij + 1
2724  vec(ij) = mat(i, j)
2725  END DO
2726  END DO
2727 
2728  END SUBROUTINE ipack
2729 
2730 ! **************************************************************************************************
2731 !> \brief Unpack a 1-D real vector as a n x n square matrix.
2732 !> \param mat resulting matrix
2733 !> \param vec vector to unpack
2734 !> \param n size of the matrix
2735 ! **************************************************************************************************
2736  PURE SUBROUTINE iunpack(mat, vec, n)
2737  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: mat
2738  REAL(dp), DIMENSION(:), INTENT(IN) :: vec
2739  INTEGER, INTENT(IN) :: n
2740 
2741  INTEGER :: i, ij, j
2742 
2743  ij = 0
2744  DO i = 1, n
2745  DO j = i, n
2746  ij = ij + 1
2747  mat(i, j) = vec(ij)
2748  mat(j, i) = vec(ij)
2749  END DO
2750  END DO
2751 
2752  END SUBROUTINE iunpack
2753 
2754 END MODULE atom_utils
imax
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
Definition: dbm_distribution.c:73
ai_onecenter
Definition: ai_onecenter.F:28
ai_onecenter::sto_overlap
subroutine, public sto_overlap(smat, na, pa, nb, pb)
...
Definition: ai_onecenter.F:667
ai_onecenter::sg_overlap
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition: ai_onecenter.F:64
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
ai_overlap::overlap_ab_s
subroutine, public overlap_ab_s(la, zeta, lb, zetb, rab, sab)
Calculation of the two-center overlap integrals [a|b] over Spherical Gaussian-type functions.
Definition: ai_overlap.F:1965
ai_overlap::overlap_ab_sp
subroutine, public overlap_ab_sp(la, zeta, lb, zetb, alat, sab)
Calculation of the overlap integrals [a|b] over cubic periodic Spherical Gaussian-type functions.
Definition: ai_overlap.F:2013
ao_util
All kind of helpful little routines.
Definition: ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition: ao_util.F:96
atom_types
Define the atom type and its sub types.
Definition: atom_types.F:15
atom_types::num_basis
integer, parameter, public num_basis
Definition: atom_types.F:69
atom_types::cgto_basis
integer, parameter, public cgto_basis
Definition: atom_types.F:69
atom_types::gto_basis
integer, parameter, public gto_basis
Definition: atom_types.F:69
atom_types::sto_basis
integer, parameter, public sto_basis
Definition: atom_types.F:69
atom_types::lmat
integer, parameter, public lmat
Definition: atom_types.F:67
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_condnumber
subroutine, public atom_condnumber(basis, crad, iw)
Print condition numbers of the given atomic basis set.
Definition: atom_utils.F:2310
atom_utils::atom_denmat
pure subroutine, public atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
Calculate a density matrix using atomic orbitals.
Definition: atom_utils.F:328
atom_utils::atom_local_potential
pure subroutine, public atom_local_potential(locpot, gthpot, rr)
...
Definition: atom_utils.F:800
atom_utils::contract2
subroutine, public contract2(int, omat, cm)
Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
Definition: atom_utils.F:2613
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::atom_orbital_charge
pure subroutine, public atom_orbital_charge(charge, wfn, rcov, l, basis)
...
Definition: atom_utils.F:654
atom_utils::atom_consistent_method
pure logical function, public atom_consistent_method(method, multiplicity)
Check that the atomic multiplicity is consistent with the electronic structure method.
Definition: atom_utils.F:2189
atom_utils::atom_core_density
subroutine, public atom_core_density(corden, potential, typ, rr)
...
Definition: atom_utils.F:694
atom_utils::get_maxn_occ
pure integer function, dimension(0:lmat), public get_maxn_occ(occupation)
Return the maximum principal quantum number of occupied orbitals.
Definition: atom_utils.F:302
atom_utils::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::atom_basis_condnum
subroutine, public atom_basis_condnum(basis, rad, cnum)
Calculate the condition number of the given atomic basis set.
Definition: atom_utils.F:2480
atom_utils::contract4
subroutine, public contract4(eri, omat, cm1, cm2)
Contract a matrix of Electron Repulsion Integrals (ERI-s).
Definition: atom_utils.F:2665
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_orbital_max
pure subroutine, public atom_orbital_max(rmax, wfn, rcov, l, basis)
...
Definition: atom_utils.F:843
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::get_rho0
subroutine, public get_rho0(atom, rho0)
Calculate the total electron density at R=0.
Definition: atom_utils.F:2219
atom_utils::atom_completeness
subroutine, public atom_completeness(basis, zv, iw)
Estimate completeness of the given atomic basis set.
Definition: atom_utils.F:2345
atom_utils::atom_write_zmp_restart
subroutine, public atom_write_zmp_restart(atom)
ZMP subroutine to write external restart file.
Definition: atom_utils.F:424
atom_utils::atom_set_occupation
subroutine, public atom_set_occupation(ostring, occupation, wfnocc, multiplicity)
Set occupation of atomic orbitals.
Definition: atom_utils.F:104
atom_utils::contract2add
subroutine, public contract2add(int, omat, cm)
Same as contract2(), but add the new contracted matrix to the old one instead of overwriting it.
Definition: atom_utils.F:2639
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_utils::atom_orbital_nodes
pure subroutine, public atom_orbital_nodes(node, wfn, rcov, l, basis)
...
Definition: atom_utils.F:882
atom_utils::get_maxl_occ
pure integer function, public get_maxl_occ(occupation)
Return the maximum orbital quantum number of occupied orbitals.
Definition: atom_utils.F:282
atom_utils::atom_wfnr0
pure subroutine, public atom_wfnr0(value, wfn, basis)
...
Definition: atom_utils.F:917
atom
Definition: atom.F:9
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::srules
pure real(dp) function, public srules(z, ne, n, l)
...
Definition: basis_set_types.F:2509
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::get_unit_number
integer function, public get_unit_number(file_name)
Returns the first logical unit that is not preconnected.
Definition: cp_files.F:237
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_rhf_atom
integer, parameter, public do_rhf_atom
Definition: input_constants.F:159
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition: input_constants.F:159
input_constants::do_uhf_atom
integer, parameter, public do_uhf_atom
Definition: input_constants.F:159
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition: input_constants.F:159
input_constants::do_rohf_atom
integer, parameter, public do_rohf_atom
Definition: input_constants.F:159
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
lapack
Interface to the LAPACK F77 library.
Definition: lapack.F:17
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition: mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition: mathconstants.F:114
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
mathconstants::maxfac
integer, parameter, public maxfac
Definition: mathconstants.F:36
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::invmat_symm
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:579
mathlib::binomial_gen
elemental real(kind=dp) function, public binomial_gen(z, k)
The generalized binomial coefficient z over k for 0 <= k <= n is calculated. (z) z*(z-1)*....
Definition: mathlib.F:231
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
orbital_pointers::deallocate_orbital_pointers
subroutine, public deallocate_orbital_pointers()
Deallocate the orbital pointers.
Definition: orbital_pointers.F:211
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition: orbital_transformation_matrices.F:17
orbital_transformation_matrices::init_spherical_harmonics
subroutine, public init_spherical_harmonics(maxl, output_unit)
Initialize or update the orbital transformation matrices.
Definition: orbital_transformation_matrices.F:256
orbital_transformation_matrices::deallocate_spherical_harmonics
subroutine, public deallocate_spherical_harmonics()
Deallocate the orbital transformation matrices.
Definition: orbital_transformation_matrices.F:226
periodic_table
Periodic Table related data definitions.
Definition: periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition: periodic_table.F:65
periodic_table::nelem
integer, parameter, public nelem
Definition: periodic_table.F:63
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
qs_grid_atom
Definition: qs_grid_atom.F:8
splines
Simple splines Splines are fully specified by the interpolation points, except that at the ends,...
Definition: splines.F:28
splines::spline3ders
subroutine, public spline3ders(x, y, xnew, ynew, dynew, d2ynew)
...
Definition: splines.F:82
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361