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