(git:f2099e5)
Loading...
Searching...
No Matches
atom_grb.F
Go to the documentation of this file.
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
8MODULE atom_grb
9 USE ai_onecenter, ONLY: sg_conf,&
10 sg_kinetic,&
11 sg_nuclear,&
12 sg_overlap
13 USE atom_electronic_structure, ONLY: calculate_atom
14 USE atom_operators, ONLY: atom_int_release,&
15 atom_int_setup,&
16 atom_ppint_release,&
17 atom_ppint_setup,&
18 atom_relint_release,&
19 atom_relint_setup
20 USE atom_types, ONLY: &
21 cgto_basis, gto_basis, atom_basis_type, atom_integrals, atom_orbitals, atom_p_type, &
22 atom_potential_type, atom_state, atom_type, create_atom_orbs, create_atom_type, lmat, &
23 release_atom_basis, release_atom_type, set_atom
24 USE atom_utils, ONLY: atom_basis_condnum,&
25 atom_density
26 USE cp_files, ONLY: close_file,&
27 open_file
28 USE input_constants, ONLY: barrier_conf,&
29 do_analytic,&
30 do_rhf_atom,&
31 do_rks_atom,&
32 do_rohf_atom,&
33 do_uhf_atom,&
34 do_uks_atom
35 USE input_section_types, ONLY: section_vals_get_subs_vals,&
36 section_vals_type,&
37 section_vals_val_get
38 USE kinds, ONLY: default_string_length,&
39 dp
40 USE mathconstants, ONLY: dfac,&
41 rootpi
42 USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
43 init_orbital_pointers
44 USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
45 init_spherical_harmonics
46 USE periodic_table, ONLY: ptable
47 USE physcon, ONLY: bohr
48 USE powell, ONLY: opt_state_type,&
49 powell_optimize
50 USE qs_grid_atom, ONLY: allocate_grid_atom,&
51 create_grid_atom
52#include "./base/base_uses.f90"
53
54 IMPLICIT NONE
55
56 TYPE basis_p_type
57 TYPE(atom_basis_type), POINTER :: basis => null()
58 END TYPE basis_p_type
59
60 PRIVATE
61 PUBLIC :: atom_grb_construction
62
63 CHARACTER(len=*), PARAMETER, PRIVATE :: modulen = 'atom_grb'
64
65CONTAINS
66
67! **************************************************************************************************
68!> \brief Construct geometrical response basis set.
69!> \param atom_info information about the atomic kind. Two-dimensional array of size
70!> (electronic-configuration, electronic-structure-method)
71!> \param atom_section ATOM input section
72!> \param iw output file unit
73!> \par History
74!> * 11.2016 created [Juerg Hutter]
75! **************************************************************************************************
76 SUBROUTINE atom_grb_construction(atom_info, atom_section, iw)
77
78 TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
79 TYPE(section_vals_type), POINTER :: atom_section
80 INTEGER, INTENT(IN) :: iw
81
82 REAL(kind=dp), PARAMETER :: error_threshold = 1.0e-12_dp
83
84 CHARACTER(len=default_string_length) :: abas, basname
85 CHARACTER(len=default_string_length), DIMENSION(1) :: basline
86 CHARACTER(len=default_string_length), DIMENSION(3) :: headline
87 INTEGER :: i, ider, is, iunit, j, k, l, lhomo, ll, &
88 lval, m, maxl, mb, method, mo, n, &
89 nder, ngp, nhomo, nr, num_gto, &
90 num_pol, quadtype, s1, s2
91 INTEGER, DIMENSION(0:7) :: nbas
92 INTEGER, DIMENSION(0:lmat) :: next_bas, next_prim
93 INTEGER, DIMENSION(:), POINTER :: num_bas
94 REAL(kind=dp) :: al, amin, aval, cnum, crad, cradx, cval, delta, dene, ear, emax, &
95 energy_ex(0:lmat), energy_ref, energy_vb(0:lmat), expzet, fhomo, o, prefac, rconf, rk, &
96 rmax, scon, zeta, zval
97 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ale, alp, rho
98 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
99 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ebasis, pbasis, qbasis, rbasis
100 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
101 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ovlp
102 TYPE(atom_basis_type), POINTER :: basis, basis_grb, basis_ref, basis_vrb
103 TYPE(atom_integrals), POINTER :: atint
104 TYPE(atom_orbitals), POINTER :: orbitals
105 TYPE(atom_state), POINTER :: state
106 TYPE(atom_type), POINTER :: atom, atom_ref, atom_test
107 TYPE(basis_p_type), DIMENSION(0:10) :: vbasis
108 TYPE(section_vals_type), POINTER :: grb_section, powell_section
109
110 IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T28,A,/," ",79("*"))') "GEOMETRICAL RESPONSE BASIS"
111
112 DO i = 0, 10
113 NULLIFY (vbasis(i)%basis)
114 END DO
115 ! make some basic checks
116 is = SIZE(atom_info)
117 IF (iw > 0 .AND. is > 1) THEN
118 WRITE (iw, '(/,A,/)') " WARNING: Only use first electronic structure/method for basis set generation"
119 END IF
120 atom_ref => atom_info(1, 1)%atom
121
122 ! check method
123 method = atom_ref%method_type
124 SELECT CASE (method)
125 CASE (do_rks_atom, do_rhf_atom)
126 ! restricted methods are okay
127 CASE (do_uks_atom, do_uhf_atom, do_rohf_atom)
128 cpabort("Unrestricted methods not allowed for GRB generation")
129 CASE DEFAULT
130 cpabort("Unknown method for GRB generation")
131 END SELECT
132
133 ! input for basis optimization
134 grb_section => section_vals_get_subs_vals(atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
135
136 ! generate an atom type
137 NULLIFY (atom)
138 CALL create_atom_type(atom)
139 CALL copy_atom_basics(atom_ref, atom, state=.true., potential=.true., optimization=.true., xc=.true.)
140 ! set confinement potential
141 atom%potential%confinement = .true.
142 atom%potential%conf_type = barrier_conf
143 atom%potential%acon = 200._dp
144 atom%potential%rcon = 4._dp
145 CALL section_vals_val_get(grb_section, "CONFINEMENT", r_val=scon)
146 atom%potential%scon = scon
147 ! generate main block geometrical exponents
148 basis_ref => atom_ref%basis
149 ALLOCATE (basis)
150 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
151 ! get information on quadrature type and number of grid points
152 ! allocate and initialize the atomic grid
153 NULLIFY (basis%grid)
154 CALL allocate_grid_atom(basis%grid)
155 CALL section_vals_val_get(grb_section, "QUADRATURE", i_val=quadtype)
156 CALL section_vals_val_get(grb_section, "GRID_POINTS", i_val=ngp)
157 IF (ngp <= 0) THEN
158 cpabort("# point radial grid < 0")
159 END IF
160 CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
161 basis%grid%nr = ngp
162 !
163 maxl = atom%state%maxl_occ
164 basis%basis_type = gto_basis
165 CALL section_vals_val_get(grb_section, "NUM_GTO_CORE", i_val=num_gto)
166 basis%nbas = 0
167 basis%nbas(0:maxl) = num_gto
168 basis%nprim = basis%nbas
169 CALL section_vals_val_get(grb_section, "GEOMETRICAL_FACTOR", r_val=cval)
170 CALL section_vals_val_get(grb_section, "GEO_START_VALUE", r_val=aval)
171 m = maxval(basis%nbas)
172 ALLOCATE (basis%am(m, 0:lmat))
173 basis%am = 0._dp
174 DO l = 0, lmat
175 DO i = 1, basis%nbas(l)
176 ll = i - 1
177 basis%am(i, l) = aval*cval**(ll)
178 END DO
179 END DO
180
181 basis%eps_eig = basis_ref%eps_eig
182 basis%geometrical = .true.
183 basis%aval = aval
184 basis%cval = cval
185 basis%start = 0
186
187 ! initialize basis function on a radial grid
188 nr = basis%grid%nr
189 m = maxval(basis%nbas)
190 ALLOCATE (basis%bf(nr, m, 0:lmat))
191 ALLOCATE (basis%dbf(nr, m, 0:lmat))
192 ALLOCATE (basis%ddbf(nr, m, 0:lmat))
193 basis%bf = 0._dp
194 basis%dbf = 0._dp
195 basis%ddbf = 0._dp
196 DO l = 0, lmat
197 DO i = 1, basis%nbas(l)
198 al = basis%am(i, l)
199 DO k = 1, nr
200 rk = basis%grid%rad(k)
201 ear = exp(-al*basis%grid%rad(k)**2)
202 basis%bf(k, i, l) = rk**l*ear
203 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
204 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
205 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
206 END DO
207 END DO
208 END DO
209
210 NULLIFY (orbitals)
211 mo = maxval(atom%state%maxn_calc)
212 mb = maxval(basis%nbas)
213 CALL create_atom_orbs(orbitals, mb, mo)
214 CALL set_atom(atom, orbitals=orbitals)
215
216 powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
217 CALL atom_fit_grb(atom, basis, iw, powell_section)
218 CALL set_atom(atom, basis=basis)
219
220 ! generate response contractions
221 CALL section_vals_val_get(grb_section, "DELTA_CHARGE", r_val=delta)
222 CALL section_vals_val_get(grb_section, "DERIVATIVES", i_val=nder)
223 IF (iw > 0) THEN
224 WRITE (iw, '(/,A,T76,I5)') " Generate Response Basis Sets with Order ", nder
225 END IF
226
227 state => atom%state
228 ! find HOMO
229 lhomo = -1
230 nhomo = -1
231 emax = -huge(1._dp)
232 DO l = 0, state%maxl_occ
233 DO i = 1, state%maxn_occ(l)
234 IF (atom%orbitals%ener(i, l) > emax) THEN
235 lhomo = l
236 nhomo = i
237 emax = atom%orbitals%ener(i, l)
238 fhomo = state%occupation(l, i)
239 END IF
240 END DO
241 END DO
242
243 s1 = SIZE(atom%orbitals%wfn, 1)
244 s2 = SIZE(atom%orbitals%wfn, 2)
245 ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
246 s2 = maxval(state%maxn_occ) + nder
247 ALLOCATE (rbasis(s1, s2, 0:lmat), qbasis(s1, s2, 0:lmat))
248 rbasis = 0._dp
249 qbasis = 0._dp
250
251 ! calculate integrals
252 ALLOCATE (atint)
253 CALL atom_int_setup(atint, basis, potential=atom%potential, eri_coulomb=.false., eri_exchange=.false.)
254 CALL atom_ppint_setup(atint, basis, potential=atom%potential)
255 IF (atom%pp_calc) THEN
256 NULLIFY (atint%tzora, atint%hdkh)
257 ELSE
258 ! relativistic correction terms
259 CALL atom_relint_setup(atint, basis, atom%relativistic, zcore=real(atom%z, dp))
260 END IF
261 CALL set_atom(atom, integrals=atint)
262
263 CALL calculate_atom(atom, iw=0)
264 DO ider = -nder, nder
265 dene = real(ider, kind=dp)*delta
266 cpassert(fhomo > abs(dene))
267 state%occupation(lhomo, nhomo) = fhomo + dene
268 CALL calculate_atom(atom, iw=0, noguess=.true.)
269 wfn(:, :, :, ider) = atom%orbitals%wfn
270 state%occupation(lhomo, nhomo) = fhomo
271 END DO
272 IF (iw > 0) THEN
273 WRITE (iw, '(A,T76,I5)') " Total number of electronic structure calculations ", 2*nder + 1
274 END IF
275
276 ovlp => atom%integrals%ovlp
277
278 DO l = 0, state%maxl_occ
279 IF (iw > 0) THEN
280 WRITE (iw, '(A,T76,I5)') " Response derivatives for l quantum number ", l
281 END IF
282 ! occupied states
283 DO i = 1, max(state%maxn_occ(l), 1)
284 rbasis(:, i, l) = wfn(:, i, l, 0)
285 END DO
286 ! differentiation
287 DO ider = 1, nder
288 i = max(state%maxn_occ(l), 1)
289 SELECT CASE (ider)
290 CASE (1)
291 rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
292 CASE (2)
293 rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
294 CASE (3)
295 rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
296 + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
297 CASE DEFAULT
298 cpabort("Only 1, 2, 3 are supported as the number of response derivatives")
299 END SELECT
300 END DO
301
302 ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
303 n = state%maxn_occ(l) + nder
304 m = atom%basis%nbas(l)
305 DO i = 1, n
306 DO j = 1, i - 1
307 o = dot_product(rbasis(1:m, j, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
308 rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
309 END DO
310 o = dot_product(rbasis(1:m, i, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
311 rbasis(1:m, i, l) = rbasis(1:m, i, l)/sqrt(o)
312 END DO
313
314 ! check
315 ALLOCATE (amat(n, n))
316 amat(1:n, 1:n) = matmul(transpose(rbasis(1:m, 1:n, l)), matmul(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
317 DO i = 1, n
318 amat(i, i) = amat(i, i) - 1._dp
319 END DO
320 IF (maxval(abs(amat)) > error_threshold) THEN
321 IF (iw > 0) WRITE (iw, '(A,G20.10)') " Orthogonality error ", maxval(abs(amat))
322 END IF
323 DEALLOCATE (amat)
324
325 ! Quickstep normalization
326 expzet = 0.25_dp*real(2*l + 3, dp)
327 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
328 DO i = 1, m
329 zeta = (2._dp*atom%basis%am(i, l))**expzet
330 qbasis(i, 1:n, l) = rbasis(i, 1:n, l)*prefac/zeta
331 END DO
332
333 END DO
334
335 ! check for condition numbers
336 IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Valence Response Basis Sets"
337 CALL init_orbital_pointers(lmat)
338 CALL init_spherical_harmonics(lmat, 0)
339 DO ider = 0, nder
340 NULLIFY (basis_vrb)
341 ALLOCATE (basis_vrb)
342 NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
343 basis_vrb%dbf, basis_vrb%ddbf)
344 ! allocate and initialize the atomic grid
345 NULLIFY (basis_vrb%grid)
346 CALL allocate_grid_atom(basis_vrb%grid)
347 CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
348 basis_vrb%grid%nr = ngp
349 !
350 basis_vrb%eps_eig = basis_ref%eps_eig
351 basis_vrb%geometrical = .false.
352 basis_vrb%basis_type = cgto_basis
353 basis_vrb%nprim = basis%nprim
354 basis_vrb%nbas = 0
355 DO l = 0, state%maxl_occ
356 basis_vrb%nbas(l) = state%maxn_occ(l) + ider
357 END DO
358 m = maxval(basis_vrb%nprim)
359 n = maxval(basis_vrb%nbas)
360 ALLOCATE (basis_vrb%am(m, 0:lmat))
361 basis_vrb%am = basis%am
362 ! contractions
363 ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
364 DO l = 0, state%maxl_occ
365 m = basis_vrb%nprim(l)
366 n = basis_vrb%nbas(l)
367 basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
368 END DO
369
370 ! initialize basis function on a radial grid
371 nr = basis_vrb%grid%nr
372 m = maxval(basis_vrb%nbas)
373 ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
374 ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
375 ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
376 basis_vrb%bf = 0._dp
377 basis_vrb%dbf = 0._dp
378 basis_vrb%ddbf = 0._dp
379 DO l = 0, lmat
380 DO i = 1, basis_vrb%nprim(l)
381 al = basis_vrb%am(i, l)
382 DO k = 1, nr
383 rk = basis_vrb%grid%rad(k)
384 ear = exp(-al*basis_vrb%grid%rad(k)**2)
385 DO j = 1, basis_vrb%nbas(l)
386 basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
387 basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
388 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
389 basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
390 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + &
391 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
392 END DO
393 END DO
394 END DO
395 END DO
396
397 IF (iw > 0) THEN
398 CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
399 WRITE (iw, '(A,A)') " Basis set ", trim(abas)
400 END IF
401 crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
402 cradx = crad*1.00_dp
403 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
404 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
405 cradx = crad*1.10_dp
406 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
407 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
408 cradx = crad*1.20_dp
409 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
410 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
411 vbasis(ider)%basis => basis_vrb
412 END DO
413 CALL deallocate_orbital_pointers
414 CALL deallocate_spherical_harmonics
415
416 ! get density maximum
417 ALLOCATE (rho(basis%grid%nr))
418 CALL calculate_atom(atom, iw=0, noguess=.true.)
419 CALL atom_density(rho(:), atom%orbitals%pmat, atom%basis, maxl, typ="RHO")
420 n = sum(maxloc(rho(:)))
421 rmax = basis%grid%rad(n)
422 IF (rmax < 0.1_dp) rmax = 1.0_dp
423 DEALLOCATE (rho)
424
425 ! generate polarization sets
426 maxl = atom%state%maxl_occ
427 CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto)
428 num_pol = num_gto
429 IF (num_gto > 0) THEN
430 IF (iw > 0) THEN
431 WRITE (iw, '(/,A)') " Polarization basis set "
432 END IF
433 ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto))
434 pbasis = 0.0_dp
435 ! optimize exponents
436 lval = maxl + 1
437 zval = sqrt(real(2*lval + 2, dp))*real(lval + 1, dp)/(2._dp*rmax)
438 aval = atom%basis%am(1, 0)
439 cval = 2.5_dp
440 rconf = atom%potential%scon
441 CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section)
442 ! calculate contractions
443 DO i = 1, num_gto
444 alp(i) = aval*cval**(i - 1)
445 END DO
446 ALLOCATE (rho(num_gto))
447 DO l = maxl + 1, min(maxl + num_gto, 7)
448 zval = sqrt(real(2*l + 2, dp))*real(l + 1, dp)/(2._dp*rmax)
449 CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l))
450 IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
451 " Polarization basis set contraction for lval=", l, "zval=", zval
452 END DO
453 DEALLOCATE (rho)
454 END IF
455
456 ! generate valence expansion sets
457 maxl = atom%state%maxl_occ
458 CALL section_vals_val_get(grb_section, "NUM_GTO_EXTENDED", i_val=num_gto)
459 CALL section_vals_val_get(grb_section, "EXTENSION_BASIS", i_vals=num_bas)
460 next_bas(0:lmat) = 0
461 IF (num_bas(1) == -1) THEN
462 DO l = 0, maxl
463 next_bas(l) = maxl - l + 1
464 END DO
465 ELSE
466 n = min(SIZE(num_bas, 1), 4)
467 next_bas(0:n - 1) = num_bas(1:n)
468 END IF
469 next_prim = 0
470 DO l = 0, lmat
471 IF (next_bas(l) > 0) next_prim(l) = num_gto
472 END DO
473 IF (iw > 0) THEN
474 CALL basis_label(abas, next_prim, next_bas)
475 WRITE (iw, '(/,A,A)') " Extension basis set ", trim(abas)
476 END IF
477 n = maxval(next_prim)
478 m = maxval(next_bas)
479 ALLOCATE (ebasis(n, n, 0:lmat), ale(n))
480 basis_vrb => vbasis(0)%basis
481 amin = atom%basis%aval/atom%basis%cval**1.5_dp
482 DO i = 1, n
483 ale(i) = amin*atom%basis%cval**(i - 1)
484 END DO
485 ebasis = 0._dp
486 ALLOCATE (rho(n))
487 rconf = 2.0_dp*atom%potential%scon
488 DO l = 0, lmat
489 IF (next_bas(l) < 1) cycle
490 zval = sqrt(real(2*l + 2, dp))*real(l + 1, dp)/(2._dp*rmax)
491 CALL hydrogenic(zval, rconf, l, ale, n, rho, ebasis(:, :, l))
492 IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
493 " Extension basis set contraction for lval=", l, "zval=", zval
494 END DO
495 DEALLOCATE (rho)
496 ! check for condition numbers
497 IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Extended Basis Sets"
498 CALL init_orbital_pointers(lmat)
499 CALL init_spherical_harmonics(lmat, 0)
500 DO ider = 0, nder
501 NULLIFY (basis_vrb)
502 ALLOCATE (basis_vrb)
503 NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
504 basis_vrb%dbf, basis_vrb%ddbf)
505 ! allocate and initialize the atomic grid
506 NULLIFY (basis_vrb%grid)
507 CALL allocate_grid_atom(basis_vrb%grid)
508 CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
509 basis_vrb%grid%nr = ngp
510 !
511 basis_vrb%eps_eig = basis_ref%eps_eig
512 basis_vrb%geometrical = .false.
513 basis_vrb%basis_type = cgto_basis
514 basis_vrb%nprim = basis%nprim + next_prim
515 basis_vrb%nbas = 0
516 DO l = 0, state%maxl_occ
517 basis_vrb%nbas(l) = state%maxn_occ(l) + ider + next_bas(l)
518 END DO
519 m = maxval(basis_vrb%nprim)
520 ALLOCATE (basis_vrb%am(m, 0:lmat))
521 ! exponents
522 m = SIZE(basis%am, 1)
523 basis_vrb%am(1:m, :) = basis%am(1:m, :)
524 n = SIZE(ale, 1)
525 DO l = 0, state%maxl_occ
526 basis_vrb%am(m + 1:m + n, l) = ale(1:n)
527 END DO
528 ! contractions
529 m = maxval(basis_vrb%nprim)
530 n = maxval(basis_vrb%nbas)
531 ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
532 basis_vrb%cm = 0.0_dp
533 DO l = 0, state%maxl_occ
534 m = basis%nprim(l)
535 n = state%maxn_occ(l) + ider
536 basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
537 basis_vrb%cm(m + 1:m + next_prim(l), n + 1:n + next_bas(l), l) = ebasis(1:next_prim(l), 1:next_bas(l), l)
538 END DO
539
540 ! initialize basis function on a radial grid
541 nr = basis_vrb%grid%nr
542 m = maxval(basis_vrb%nbas)
543 ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
544 ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
545 ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
546 basis_vrb%bf = 0._dp
547 basis_vrb%dbf = 0._dp
548 basis_vrb%ddbf = 0._dp
549 DO l = 0, lmat
550 DO i = 1, basis_vrb%nprim(l)
551 al = basis_vrb%am(i, l)
552 DO k = 1, nr
553 rk = basis_vrb%grid%rad(k)
554 ear = exp(-al*basis_vrb%grid%rad(k)**2)
555 DO j = 1, basis_vrb%nbas(l)
556 basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
557 basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
558 + (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
559 basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
560 (real(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1, dp)*rk**(l) + &
561 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
562 END DO
563 END DO
564 END DO
565 END DO
566
567 IF (iw > 0) THEN
568 CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
569 WRITE (iw, '(A,A)') " Basis set ", trim(abas)
570 END IF
571 crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
572 cradx = crad*1.00_dp
573 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
574 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
575 cradx = crad*1.10_dp
576 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
577 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
578 cradx = crad*1.20_dp
579 CALL atom_basis_condnum(basis_vrb, cradx, cnum)
580 IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
581 vbasis(nder + 1 + ider)%basis => basis_vrb
582 END DO
583 CALL deallocate_orbital_pointers
584 CALL deallocate_spherical_harmonics
585
586 ! Tests for energy
587 energy_ref = atom_ref%energy%etot
588 IF (iw > 0) WRITE (iw, '(/,A,A)') " Basis set tests "
589 IF (iw > 0) WRITE (iw, '(T10,A,T59,F22.9)') " Reference Energy [a.u.] ", energy_ref
590 DO ider = 0, 2*nder + 1
591 ! generate an atom type
592 NULLIFY (atom_test)
593 CALL create_atom_type(atom_test)
594 CALL copy_atom_basics(atom_ref, atom_test, state=.true., potential=.true., &
595 optimization=.true., xc=.true.)
596 basis_grb => vbasis(ider)%basis
597 NULLIFY (orbitals)
598 mo = maxval(atom_test%state%maxn_calc)
599 mb = maxval(basis_grb%nbas)
600 CALL create_atom_orbs(orbitals, mb, mo)
601 CALL set_atom(atom_test, orbitals=orbitals, basis=basis_grb)
602 ! calculate integrals
603 ALLOCATE (atint)
604 CALL atom_int_setup(atint, basis_grb, potential=atom_test%potential, eri_coulomb=.false., eri_exchange=.false.)
605 CALL atom_ppint_setup(atint, basis_grb, potential=atom_test%potential)
606 IF (atom_test%pp_calc) THEN
607 NULLIFY (atint%tzora, atint%hdkh)
608 ELSE
609 ! relativistic correction terms
610 CALL atom_relint_setup(atint, basis_grb, atom_test%relativistic, zcore=real(atom_test%z, dp))
611 END IF
612 CALL set_atom(atom_test, integrals=atint)
613 !
614 CALL calculate_atom(atom_test, iw=0)
615 IF (ider <= nder) THEN
616 energy_vb(ider) = atom_test%energy%etot
617 IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (VB)", ider, " Energy [a.u.] ", &
618 energy_ref - energy_vb(ider), energy_vb(ider)
619 ELSE
620 i = ider - nder - 1
621 energy_ex(i) = atom_test%energy%etot
622 IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (EX)", i, " Energy [a.u.] ", &
623 energy_ref - energy_ex(i), energy_ex(i)
624 END IF
625 CALL atom_int_release(atint)
626 CALL atom_ppint_release(atint)
627 CALL atom_relint_release(atint)
628 DEALLOCATE (atom_test%state, atom_test%potential, atint)
629 CALL release_atom_type(atom_test)
630 END DO
631
632 ! Quickstep normalization polarization basis
633 DO l = 0, 7
634 expzet = 0.25_dp*real(2*l + 3, dp)
635 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
636 DO i = 1, num_pol
637 zeta = (2._dp*alp(i))**expzet
638 pbasis(i, 1:num_pol, l) = pbasis(i, 1:num_pol, l)*prefac/zeta
639 END DO
640 END DO
641 ! Quickstep normalization extended basis
642 DO l = 0, lmat
643 expzet = 0.25_dp*real(2*l + 3, dp)
644 prefac = sqrt(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
645 DO i = 1, next_prim(l)
646 zeta = (2._dp*ale(i))**expzet
647 ebasis(i, 1:next_bas(l), l) = ebasis(i, 1:next_bas(l), l)*prefac/zeta
648 END DO
649 END DO
650
651 ! Print basis sets
652 CALL section_vals_val_get(grb_section, "NAME_BODY", c_val=basname)
653 CALL open_file(file_name="GRB_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iunit)
654 ! header info
655 headline = ""
656 headline(1) = "#"
657 headline(2) = "# Generated with CP2K Atom Code"
658 headline(3) = "#"
659 CALL grb_print_basis(header=headline, iunit=iunit)
660 ! valence basis
661 basline(1) = ""
662 WRITE (basline(1), "(T2,A)") adjustl(ptable(atom_ref%z)%symbol)
663 DO ider = 0, nder
664 basline(1) = ""
665 WRITE (basline(1), "(T2,A,T5,A,I1)") adjustl(ptable(atom_ref%z)%symbol), trim(adjustl(basname))//"-VAL-", ider
666 CALL grb_print_basis(header=basline, nprim=vbasis(ider)%basis%nprim(0), nbas=vbasis(ider)%basis%nbas, &
667 al=vbasis(ider)%basis%am(:, 0), gcc=qbasis, iunit=iunit)
668 END DO
669 ! polarization basis
670 maxl = atom_ref%state%maxl_occ
671 DO l = maxl + 1, min(maxl + num_pol, 7)
672 nbas = 0
673 DO i = maxl + 1, l
674 nbas(i) = l - i + 1
675 END DO
676 i = l - maxl
677 basline(1) = ""
678 WRITE (basline(1), "(T2,A,T5,A,I1)") adjustl(ptable(atom_ref%z)%symbol), trim(adjustl(basname))//"-POL-", i
679 CALL grb_print_basis(header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit)
680 END DO
681 ! extension set
682 IF (sum(next_bas) > 0) THEN
683 basline(1) = ""
684 WRITE (basline(1), "(T2,A,T5,A)") adjustl(ptable(atom_ref%z)%symbol), trim(adjustl(basname))//"-EXT"
685 CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit)
686 END IF
687 !
688 CALL close_file(unit_number=iunit)
689
690 ! clean up
691 IF (ALLOCATED(pbasis)) DEALLOCATE (pbasis)
692 IF (ALLOCATED(alp)) DEALLOCATE (alp)
693 IF (ALLOCATED(ebasis)) DEALLOCATE (ebasis)
694 DEALLOCATE (wfn, rbasis, qbasis, ale)
695
696 DO ider = 0, 10
697 IF (ASSOCIATED(vbasis(ider)%basis)) THEN
698 CALL release_atom_basis(vbasis(ider)%basis)
699 DEALLOCATE (vbasis(ider)%basis)
700 END IF
701 END DO
702
703 CALL atom_int_release(atom%integrals)
704 CALL atom_ppint_release(atom%integrals)
705 CALL atom_relint_release(atom%integrals)
706 CALL release_atom_basis(basis)
707 DEALLOCATE (atom%potential, atom%state, atom%integrals, basis)
708 CALL release_atom_type(atom)
709
710 IF (iw > 0) WRITE (iw, '(" ",79("*"))')
711
712 END SUBROUTINE atom_grb_construction
713
714! **************************************************************************************************
715!> \brief Print geometrical response basis set.
716!> \param header banner to print on top of the basis set
717!> \param nprim number of primitive exponents
718!> \param nbas number of basis functions for the given angular momentum
719!> \param al list of the primitive exponents
720!> \param gcc array of contraction coefficients of size
721!> (index-of-the-primitive-exponent, index-of-the-contraction-set, angular-momentum)
722!> \param iunit output file unit
723!> \par History
724!> * 11.2016 created [Juerg Hutter]
725! **************************************************************************************************
726 SUBROUTINE grb_print_basis(header, nprim, nbas, al, gcc, iunit)
727 CHARACTER(len=*), DIMENSION(:), INTENT(IN), &
728 OPTIONAL :: header
729 INTEGER, INTENT(IN), OPTIONAL :: nprim
730 INTEGER, DIMENSION(0:), INTENT(IN), OPTIONAL :: nbas
731 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: al
732 REAL(kind=dp), DIMENSION(:, :, 0:), INTENT(IN), &
733 OPTIONAL :: gcc
734 INTEGER, INTENT(IN) :: iunit
735
736 INTEGER :: i, j, l, lmax, lmin, nval
737
738 IF (PRESENT(header)) THEN
739 DO i = 1, SIZE(header, 1)
740 IF (header(i) /= "") THEN
741 WRITE (iunit, "(A)") trim(header(i))
742 END IF
743 END DO
744 END IF
745
746 IF (PRESENT(nprim)) THEN
747 IF (nprim > 0) THEN
748 cpassert(PRESENT(nbas))
749 cpassert(PRESENT(al))
750 cpassert(PRESENT(gcc))
751
752 DO i = lbound(nbas, 1), ubound(nbas, 1)
753 IF (nbas(i) > 0) THEN
754 lmin = i
755 EXIT
756 END IF
757 END DO
758 DO i = ubound(nbas, 1), lbound(nbas, 1), -1
759 IF (nbas(i) > 0) THEN
760 lmax = i
761 EXIT
762 END IF
763 END DO
764
765 nval = lmax
766 WRITE (iunit, *) " 1"
767 WRITE (iunit, "(40I3)") nval, lmin, lmax, nprim, (nbas(l), l=lmin, lmax)
768 DO i = nprim, 1, -1
769 WRITE (iunit, "(G20.12)", advance="no") al(i)
770 DO l = lmin, lmax
771 DO j = 1, nbas(l)
772 WRITE (iunit, "(F16.10)", advance="no") gcc(i, j, l)
773 END DO
774 END DO
775 WRITE (iunit, *)
776 END DO
777 WRITE (iunit, *)
778 END IF
779 END IF
780
781 END SUBROUTINE grb_print_basis
782
783! **************************************************************************************************
784!> \brief Compose the basis set label:
785!> (np(0)'s'np(1)'p'...) -> [nb(0)'s'nb(1)'p'...] .
786!> \param label basis set label
787!> \param np number of primitive basis functions per angular momentum
788!> \param nb number of contracted basis functions per angular momentum
789!> \par History
790!> * 11.2016 created [Juerg Hutter]
791! **************************************************************************************************
792 SUBROUTINE basis_label(label, np, nb)
793 CHARACTER(len=*), INTENT(out) :: label
794 INTEGER, DIMENSION(0:), INTENT(in) :: np, nb
795
796 INTEGER :: i, l, lmax
797 CHARACTER(len=1), DIMENSION(0:7), PARAMETER :: lq = ["s", "p", "d", "f", "g", "h", "i", "k"]
798
799 label = ""
800 lmax = min(ubound(np, 1), ubound(nb, 1), 7)
801 i = 1
802 label(i:i) = "("
803 DO l = 0, lmax
804 IF (np(l) > 0) THEN
805 i = i + 1
806 IF (np(l) > 9) THEN
807 WRITE (label(i:i + 1), "(I2)") np(l)
808 i = i + 2
809 ELSE
810 WRITE (label(i:i), "(I1)") np(l)
811 i = i + 1
812 END IF
813 label(i:i) = lq(l)
814 END IF
815 END DO
816 i = i + 1
817 label(i:i + 6) = ") --> ["
818 i = i + 6
819 DO l = 0, lmax
820 IF (nb(l) > 0) THEN
821 i = i + 1
822 IF (nb(l) > 9) THEN
823 WRITE (label(i:i + 1), "(I2)") nb(l)
824 i = i + 2
825 ELSE
826 WRITE (label(i:i), "(I1)") nb(l)
827 i = i + 1
828 END IF
829 label(i:i) = lq(l)
830 END IF
831 END DO
832 i = i + 1
833 label(i:i) = "]"
834
835 END SUBROUTINE basis_label
836
837! **************************************************************************************************
838!> \brief Compute the total energy for the given atomic kind and basis set.
839!> \param atom information about the atomic kind
840!> \param basis basis set to fit
841!> \param afun (output) atomic total energy
842!> \param iw output file unit
843!> \par History
844!> * 11.2016 created [Juerg Hutter]
845! **************************************************************************************************
846 SUBROUTINE grb_fit(atom, basis, afun, iw)
847 TYPE(atom_type), POINTER :: atom
848 TYPE(atom_basis_type), POINTER :: basis
849 REAL(dp), INTENT(OUT) :: afun
850 INTEGER, INTENT(IN) :: iw
851
852 INTEGER :: do_eric, do_erie, reltyp, zval
853 LOGICAL :: eri_c, eri_e
854 TYPE(atom_integrals), POINTER :: atint
855 TYPE(atom_potential_type), POINTER :: pot
856
857 ALLOCATE (atint)
858 ! calculate integrals
859 NULLIFY (pot)
860 eri_c = .false.
861 eri_e = .false.
862 pot => atom%potential
863 zval = atom%z
864 reltyp = atom%relativistic
865 do_eric = atom%coulomb_integral_type
866 do_erie = atom%exchange_integral_type
867 IF (do_eric == do_analytic) eri_c = .true.
868 IF (do_erie == do_analytic) eri_e = .true.
869 ! general integrals
870 CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
871 ! potential
872 CALL atom_ppint_setup(atint, basis, potential=pot)
873 IF (atom%pp_calc) THEN
874 NULLIFY (atint%tzora, atint%hdkh)
875 ELSE
876 ! relativistic correction terms
877 CALL atom_relint_setup(atint, basis, reltyp, zcore=real(zval, dp))
878 END IF
879 CALL set_atom(atom, basis=basis)
880 CALL set_atom(atom, integrals=atint)
881 CALL calculate_atom(atom, iw)
882 afun = atom%energy%etot
883 CALL atom_int_release(atint)
884 CALL atom_ppint_release(atint)
885 CALL atom_relint_release(atint)
886 DEALLOCATE (atint)
887 END SUBROUTINE grb_fit
888
889! **************************************************************************************************
890!> \brief Copy basic information about the atomic kind.
891!> \param atom_ref atom to copy
892!> \param atom new atom to create
893!> \param state also copy electronic state and occupation numbers
894!> \param potential also copy pseudo-potential
895!> \param optimization also copy optimization procedure
896!> \param xc also copy the XC input section
897!> \par History
898!> * 11.2016 created [Juerg Hutter]
899! **************************************************************************************************
900 SUBROUTINE copy_atom_basics(atom_ref, atom, state, potential, optimization, xc)
901 TYPE(atom_type), POINTER :: atom_ref, atom
902 LOGICAL, INTENT(IN), OPTIONAL :: state, potential, optimization, xc
903
904 atom%z = atom_ref%z
905 atom%zcore = atom_ref%zcore
906 atom%pp_calc = atom_ref%pp_calc
907 atom%method_type = atom_ref%method_type
908 atom%relativistic = atom_ref%relativistic
909 atom%coulomb_integral_type = atom_ref%coulomb_integral_type
910 atom%exchange_integral_type = atom_ref%exchange_integral_type
911
912 NULLIFY (atom%potential, atom%state, atom%xc_section)
913 NULLIFY (atom%basis, atom%integrals, atom%orbitals, atom%fmat)
914
915 IF (PRESENT(state)) THEN
916 IF (state) THEN
917 ALLOCATE (atom%state)
918 atom%state = atom_ref%state
919 END IF
920 END IF
921
922 IF (PRESENT(potential)) THEN
923 IF (potential) THEN
924 ALLOCATE (atom%potential)
925 atom%potential = atom_ref%potential
926 END IF
927 END IF
928
929 IF (PRESENT(optimization)) THEN
930 IF (optimization) THEN
931 atom%optimization = atom_ref%optimization
932 END IF
933 END IF
934
935 IF (PRESENT(xc)) THEN
936 IF (xc) THEN
937 atom%xc_section => atom_ref%xc_section
938 END IF
939 END IF
940
941 END SUBROUTINE copy_atom_basics
942
943! **************************************************************************************************
944!> \brief Optimise a geometrical response basis set.
945!> \param atom information about the atomic kind
946!> \param basis basis set to fit
947!> \param iunit output file unit
948!> \param powell_section POWELL input section
949!> \par History
950!> * 11.2016 created [Juerg Hutter]
951! **************************************************************************************************
952 SUBROUTINE atom_fit_grb(atom, basis, iunit, powell_section)
953 TYPE(atom_type), POINTER :: atom
954 TYPE(atom_basis_type), POINTER :: basis
955 INTEGER, INTENT(IN) :: iunit
956 TYPE(section_vals_type), POINTER :: powell_section
957
958 INTEGER :: i, k, l, ll, n10, nr
959 REAL(kind=dp) :: al, cnum, crad, cradx, ear, fopt, rk
960 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: x
961 TYPE(opt_state_type) :: ostate
962
963 cpassert(basis%geometrical)
964
965 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
966 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
967 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
968
969 ostate%nvar = 2
970 ALLOCATE (x(2))
971 x(1) = sqrt(basis%aval)
972 x(2) = sqrt(basis%cval)
973
974 ostate%nf = 0
975 ostate%iprint = 1
976 ostate%unit = iunit
977
978 ostate%state = 0
979 IF (iunit > 0) THEN
980 WRITE (iunit, '(/," POWELL| Start optimization procedure")')
981 WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
982 END IF
983 n10 = max(ostate%maxfun/100, 1)
984
985 fopt = huge(0._dp)
986
987 DO
988
989 IF (ostate%state == 2) THEN
990 basis%am = 0._dp
991 DO l = 0, lmat
992 DO i = 1, basis%nbas(l)
993 ll = i - 1 + basis%start(l)
994 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
995 END DO
996 END DO
997 basis%aval = x(1)*x(1)
998 basis%cval = x(2)*x(2)
999 basis%bf = 0._dp
1000 basis%dbf = 0._dp
1001 basis%ddbf = 0._dp
1002 nr = basis%grid%nr
1003 DO l = 0, lmat
1004 DO i = 1, basis%nbas(l)
1005 al = basis%am(i, l)
1006 DO k = 1, nr
1007 rk = basis%grid%rad(k)
1008 ear = exp(-al*basis%grid%rad(k)**2)
1009 basis%bf(k, i, l) = rk**l*ear
1010 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1011 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
1012 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1013 END DO
1014 END DO
1015 END DO
1016 CALL grb_fit(atom, basis, ostate%f, 0)
1017 fopt = min(fopt, ostate%f)
1018 END IF
1019
1020 IF (ostate%state == -1) EXIT
1021
1022 CALL powell_optimize(ostate%nvar, x, ostate)
1023
1024 IF (ostate%nf == 2 .AND. iunit > 0) THEN
1025 WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1026 END IF
1027 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1028 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1029 int(real(ostate%nf, dp)/real(ostate%maxfun, dp)*100._dp), fopt
1030 END IF
1031
1032 END DO
1033
1034 ostate%state = 8
1035 CALL powell_optimize(ostate%nvar, x, ostate)
1036
1037 IF (iunit > 0) THEN
1038 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1039 WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1040 END IF
1041 ! x->basis
1042 basis%am = 0._dp
1043 DO l = 0, lmat
1044 DO i = 1, basis%nbas(l)
1045 ll = i - 1 + basis%start(l)
1046 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
1047 END DO
1048 END DO
1049 basis%aval = x(1)*x(1)
1050 basis%cval = x(2)*x(2)
1051 basis%bf = 0._dp
1052 basis%dbf = 0._dp
1053 basis%ddbf = 0._dp
1054 nr = basis%grid%nr
1055 DO l = 0, lmat
1056 DO i = 1, basis%nbas(l)
1057 al = basis%am(i, l)
1058 DO k = 1, nr
1059 rk = basis%grid%rad(k)
1060 ear = exp(-al*basis%grid%rad(k)**2)
1061 basis%bf(k, i, l) = rk**l*ear
1062 basis%dbf(k, i, l) = (real(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1063 basis%ddbf(k, i, l) = (real(l*(l - 1), dp)*rk**(l - 2) - &
1064 2._dp*al*real(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1065 END DO
1066 END DO
1067 END DO
1068
1069 DEALLOCATE (x)
1070
1071 ! final result
1072 IF (iunit > 0) THEN
1073 WRITE (iunit, '(/,A)') " Optimized Geometrical GTO basis set"
1074 WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", basis%aval, &
1075 " Proportionality factor: ", basis%cval
1076 DO l = 0, lmat
1077 WRITE (iunit, '(T41,A,I2,T76,I5)') " Number of exponents for l=", l, basis%nbas(l)
1078 END DO
1079 END IF
1080
1081 IF (iunit > 0) WRITE (iunit, '(/,A)') " Condition number of uncontracted basis set"
1082 crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
1083 CALL init_orbital_pointers(lmat)
1084 CALL init_spherical_harmonics(lmat, 0)
1085 cradx = crad*1.00_dp
1086 CALL atom_basis_condnum(basis, cradx, cnum)
1087 IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1088 cradx = crad*1.10_dp
1089 CALL atom_basis_condnum(basis, cradx, cnum)
1090 IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1091 cradx = crad*1.20_dp
1092 CALL atom_basis_condnum(basis, cradx, cnum)
1093 IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1094 CALL deallocate_orbital_pointers
1095 CALL deallocate_spherical_harmonics
1096
1097 END SUBROUTINE atom_fit_grb
1098
1099! **************************************************************************************************
1100!> \brief Optimize 'aval' and 'cval' parameters which define the geometrical response basis set.
1101!> \param zval nuclear charge
1102!> \param rconf confinement radius
1103!> \param lval angular momentum
1104!> \param aval (input/output) exponent of the first Gaussian basis function in the series
1105!> \param cval (input/output) factor of geometrical series
1106!> \param nbas number of basis functions
1107!> \param iunit output file unit
1108!> \param powell_section POWELL input section
1109!> \par History
1110!> * 11.2016 created [Juerg Hutter]
1111! **************************************************************************************************
1112 SUBROUTINE atom_fit_pol(zval, rconf, lval, aval, cval, nbas, iunit, powell_section)
1113 REAL(kind=dp), INTENT(IN) :: zval, rconf
1114 INTEGER, INTENT(IN) :: lval
1115 REAL(kind=dp), INTENT(INOUT) :: aval, cval
1116 INTEGER, INTENT(IN) :: nbas, iunit
1117 TYPE(section_vals_type), POINTER :: powell_section
1118
1119 INTEGER :: i, n10
1120 REAL(kind=dp) :: fopt, x(2)
1121 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: am, ener
1122 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1123 TYPE(opt_state_type) :: ostate
1124
1125 ALLOCATE (am(nbas), ener(nbas), orb(nbas, nbas))
1126
1127 CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
1128 CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
1129 CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
1130
1131 ostate%nvar = 2
1132 x(1) = sqrt(aval)
1133 x(2) = sqrt(cval)
1134
1135 ostate%nf = 0
1136 ostate%iprint = 1
1137 ostate%unit = iunit
1138
1139 ostate%state = 0
1140 IF (iunit > 0) THEN
1141 WRITE (iunit, '(/," POWELL| Start optimization procedure")')
1142 WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
1143 END IF
1144 n10 = max(ostate%maxfun/100, 1)
1145
1146 fopt = huge(0._dp)
1147
1148 DO
1149
1150 IF (ostate%state == 2) THEN
1151 aval = x(1)*x(1)
1152 cval = x(2)*x(2)
1153 DO i = 1, nbas
1154 am(i) = aval*cval**(i - 1)
1155 END DO
1156 CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1157 ostate%f = ener(1)
1158 fopt = min(fopt, ostate%f)
1159 END IF
1160
1161 IF (ostate%state == -1) EXIT
1162
1163 CALL powell_optimize(ostate%nvar, x, ostate)
1164
1165 IF (ostate%nf == 2 .AND. iunit > 0) THEN
1166 WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1167 END IF
1168 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1169 WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1170 int(real(ostate%nf, dp)/real(ostate%maxfun, dp)*100._dp), fopt
1171 END IF
1172
1173 END DO
1174
1175 ostate%state = 8
1176 CALL powell_optimize(ostate%nvar, x, ostate)
1177
1178 IF (iunit > 0) THEN
1179 WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1180 WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1181 END IF
1182 ! x->basis
1183 aval = x(1)*x(1)
1184 cval = x(2)*x(2)
1185
1186 ! final result
1187 IF (iunit > 0) THEN
1188 WRITE (iunit, '(/,A,T51,A,T76,I5)') " Optimized Polarization basis set", &
1189 " Number of exponents:", nbas
1190 WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", aval, &
1191 " Proportionality factor: ", cval
1192 END IF
1193
1194 DEALLOCATE (am, ener, orb)
1195
1196 END SUBROUTINE atom_fit_pol
1197
1198! **************************************************************************************************
1199!> \brief Calculate orbitals of a hydrogen-like atom.
1200!> \param zval nuclear charge
1201!> \param rconf confinement radius
1202!> \param lval angular momentum
1203!> \param am list of basis functions' exponents
1204!> \param nbas number of basis functions
1205!> \param ener orbital energies
1206!> \param orb expansion coefficients of atomic wavefunctions
1207!> \par History
1208!> * 11.2016 created [Juerg Hutter]
1209! **************************************************************************************************
1210 SUBROUTINE hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1211 REAL(kind=dp), INTENT(IN) :: zval, rconf
1212 INTEGER, INTENT(IN) :: lval
1213 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: am
1214 INTEGER, INTENT(IN) :: nbas
1215 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: ener
1216 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: orb
1217
1218 INTEGER :: info, k, lwork, n
1219 REAL(kind=dp) :: cf
1220 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: w, work
1221 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: confmat, hmat, potmat, smat, tmat
1222
1223 n = nbas
1224 ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
1225 ! calclulate overlap matrix
1226 CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
1227 ! calclulate kinetic energy matrix
1228 CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
1229 ! calclulate core potential matrix
1230 CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
1231 ! calclulate confinement potential matrix
1232 cf = 0.1_dp
1233 k = 10
1234 CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
1235 ! Hamiltionian
1236 hmat(1:n, 1:n) = tmat(1:n, 1:n) - zval*potmat(1:n, 1:n) + cf*confmat(1:n, 1:n)
1237 ! solve
1238 lwork = 100*n
1239 ALLOCATE (w(n), work(lwork))
1240 CALL dsygv(1, "V", "U", n, hmat, n, smat, n, w, work, lwork, info)
1241 cpassert(info == 0)
1242 orb(1:n, 1:n) = hmat(1:n, 1:n)
1243 ener(1:n) = w(1:n)
1244 DEALLOCATE (w, work)
1245 DEALLOCATE (smat, tmat, potmat, confmat, hmat)
1246
1247 END SUBROUTINE hydrogenic
1248
1249END MODULE atom_grb
ai_onecenter
Definition ai_onecenter.F:28
ai_onecenter::sg_nuclear
subroutine, public sg_nuclear(umat, l, pa, pb)
...
Definition ai_onecenter.F:137
ai_onecenter::sg_kinetic
subroutine, public sg_kinetic(kmat, l, pa, pb)
...
Definition ai_onecenter.F:102
ai_onecenter::sg_conf
subroutine, public sg_conf(gmat, rc, k, l, pa, pb)
...
Definition ai_onecenter.F:414
ai_onecenter::sg_overlap
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition ai_onecenter.F:65
atom_electronic_structure
Definition atom_electronic_structure.F:8
atom_electronic_structure::calculate_atom
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
Definition atom_electronic_structure.F:60
atom_grb
Definition atom_grb.F:8
atom_grb::atom_grb_construction
subroutine, public atom_grb_construction(atom_info, atom_section, iw)
Construct geometrical response basis set.
Definition atom_grb.F:77
atom_operators
Calculate the atomic operator matrices.
Definition atom_operators.F:15
atom_operators::atom_ppint_release
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
Definition atom_operators.F:1160
atom_operators::atom_int_setup
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
Definition atom_operators.F:63
atom_operators::atom_relint_setup
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
Definition atom_operators.F:817
atom_operators::atom_relint_release
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
Definition atom_operators.F:1181
atom_operators::atom_ppint_setup
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
Definition atom_operators.F:336
atom_operators::atom_int_release
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Definition atom_operators.F:1123
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::create_atom_type
subroutine, public create_atom_type(atom)
...
Definition atom_types.F:965
atom_types::cgto_basis
integer, parameter, public cgto_basis
Definition atom_types.F:69
atom_types::gto_basis
integer, parameter, public gto_basis
Definition atom_types.F:69
atom_types::release_atom_type
subroutine, public release_atom_type(atom)
...
Definition atom_types.F:989
atom_types::lmat
integer, parameter, public lmat
Definition atom_types.F:67
atom_types::set_atom
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
Definition atom_types.F:1030
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:931
atom_types::create_atom_orbs
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Definition atom_types.F:1076
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::atom_density
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition atom_utils.F:366
atom_utils::atom_basis_condnum
subroutine, public atom_basis_condnum(basis, rad, cnum)
Calculate the condition number of the given atomic basis set.
Definition atom_utils.F:2491
atom
Definition atom.F:9
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:323
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:123
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_rhf_atom
integer, parameter, public do_rhf_atom
Definition input_constants.F:170
input_constants::do_rks_atom
integer, parameter, public do_rks_atom
Definition input_constants.F:170
input_constants::do_analytic
integer, parameter, public do_analytic
Definition input_constants.F:160
input_constants::do_uhf_atom
integer, parameter, public do_uhf_atom
Definition input_constants.F:170
input_constants::do_uks_atom
integer, parameter, public do_uks_atom
Definition input_constants.F:170
input_constants::barrier_conf
integer, parameter, public barrier_conf
Definition input_constants.F:191
input_constants::do_rohf_atom
integer, parameter, public do_rohf_atom
Definition input_constants.F:170
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::deallocate_orbital_pointers
subroutine, public deallocate_orbital_pointers()
Deallocate the orbital pointers.
Definition orbital_pointers.F:211
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::init_spherical_harmonics
subroutine, public init_spherical_harmonics(maxl, output_unit)
Initialize or update the orbital transformation matrices.
Definition orbital_transformation_matrices.F:255
orbital_transformation_matrices::deallocate_spherical_harmonics
subroutine, public deallocate_spherical_harmonics()
Deallocate the orbital transformation matrices.
Definition orbital_transformation_matrices.F:225
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
powell
Definition powell.F:9
powell::powell_optimize
subroutine, public powell_optimize(n, x, optstate)
...
Definition powell.F:52
qs_grid_atom
Definition qs_grid_atom.F:8
qs_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition qs_grid_atom.F:74
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition qs_grid_atom.F:189
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
atom_grb::basis_p_type
Definition atom_grb.F:56
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_integrals
Definition atom_types.F:216
atom_types::atom_orbitals
Holds atomic orbitals and energies.
Definition atom_types.F:237
atom_types::atom_p_type
Definition atom_types.F:323
atom_types::atom_potential_type
Definition atom_types.F:171
atom_types::atom_state
Provides all information on states and occupation.
Definition atom_types.F:198
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:293
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
powell::opt_state_type
Definition powell.F:18