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 : MODULE 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 :
65 : CONTAINS
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 2 : 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 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ale, alp, rho
98 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
99 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ebasis, pbasis, qbasis, rbasis
100 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
101 2 : 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 24 : TYPE(basis_p_type), DIMENSION(0:10) :: vbasis
108 : TYPE(section_vals_type), POINTER :: grb_section, powell_section
109 :
110 2 : IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T28,A,/," ",79("*"))') "GEOMETRICAL RESPONSE BASIS"
111 :
112 24 : DO i = 0, 10
113 24 : NULLIFY (vbasis(i)%basis)
114 : END DO
115 : ! make some basic checks
116 6 : is = SIZE(atom_info)
117 2 : IF (iw > 0 .AND. is > 1) THEN
118 0 : WRITE (iw, '(/,A,/)') " WARNING: Only use first electronic structure/method for basis set generation"
119 : END IF
120 2 : atom_ref => atom_info(1, 1)%atom
121 :
122 : ! check method
123 2 : method = atom_ref%method_type
124 0 : 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 0 : CPABORT("Unrestricted methods not allowed for GRB generation")
129 : CASE DEFAULT
130 2 : CPABORT("Unknown method for GRB generation")
131 : END SELECT
132 :
133 : ! input for basis optimization
134 2 : grb_section => section_vals_get_subs_vals(atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
135 :
136 : ! generate an atom type
137 2 : NULLIFY (atom)
138 2 : CALL create_atom_type(atom)
139 2 : CALL copy_atom_basics(atom_ref, atom, state=.TRUE., potential=.TRUE., optimization=.TRUE., xc=.TRUE.)
140 : ! set confinement potential
141 2 : atom%potential%confinement = .TRUE.
142 2 : atom%potential%conf_type = barrier_conf
143 2 : atom%potential%acon = 200._dp
144 2 : atom%potential%rcon = 4._dp
145 2 : CALL section_vals_val_get(grb_section, "CONFINEMENT", r_val=scon)
146 2 : atom%potential%scon = scon
147 : ! generate main block geometrical exponents
148 2 : basis_ref => atom_ref%basis
149 38 : 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 2 : CALL allocate_grid_atom(basis%grid)
155 2 : CALL section_vals_val_get(grb_section, "QUADRATURE", i_val=quadtype)
156 2 : CALL section_vals_val_get(grb_section, "GRID_POINTS", i_val=ngp)
157 2 : IF (ngp <= 0) THEN
158 0 : CPABORT("# point radial grid < 0")
159 : END IF
160 2 : CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
161 2 : basis%grid%nr = ngp
162 : !
163 2 : maxl = atom%state%maxl_occ
164 2 : basis%basis_type = GTO_BASIS
165 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_CORE", i_val=num_gto)
166 14 : basis%nbas = 0
167 6 : basis%nbas(0:maxl) = num_gto
168 14 : basis%nprim = basis%nbas
169 2 : CALL section_vals_val_get(grb_section, "GEOMETRICAL_FACTOR", r_val=cval)
170 2 : CALL section_vals_val_get(grb_section, "GEO_START_VALUE", r_val=aval)
171 14 : m = MAXVAL(basis%nbas)
172 6 : ALLOCATE (basis%am(m, 0:lmat))
173 86 : basis%am = 0._dp
174 14 : DO l = 0, lmat
175 38 : DO i = 1, basis%nbas(l)
176 24 : ll = i - 1
177 36 : basis%am(i, l) = aval*cval**(ll)
178 : END DO
179 : END DO
180 :
181 2 : basis%eps_eig = basis_ref%eps_eig
182 2 : basis%geometrical = .TRUE.
183 2 : basis%aval = aval
184 2 : basis%cval = cval
185 14 : basis%start = 0
186 :
187 : ! initialize basis function on a radial grid
188 2 : nr = basis%grid%nr
189 14 : m = MAXVAL(basis%nbas)
190 10 : ALLOCATE (basis%bf(nr, m, 0:lmat))
191 6 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
192 6 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
193 28886 : basis%bf = 0._dp
194 28886 : basis%dbf = 0._dp
195 28886 : basis%ddbf = 0._dp
196 14 : DO l = 0, lmat
197 38 : DO i = 1, basis%nbas(l)
198 24 : al = basis%am(i, l)
199 9636 : DO k = 1, nr
200 9600 : rk = basis%grid%rad(k)
201 9600 : ear = EXP(-al*basis%grid%rad(k)**2)
202 9600 : basis%bf(k, i, l) = rk**l*ear
203 9600 : 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 9624 : 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 2 : NULLIFY (orbitals)
211 14 : mo = MAXVAL(atom%state%maxn_calc)
212 14 : mb = MAXVAL(basis%nbas)
213 2 : CALL create_atom_orbs(orbitals, mb, mo)
214 2 : CALL set_atom(atom, orbitals=orbitals)
215 :
216 2 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
217 2 : CALL atom_fit_grb(atom, basis, iw, powell_section)
218 2 : CALL set_atom(atom, basis=basis)
219 :
220 : ! generate response contractions
221 2 : CALL section_vals_val_get(grb_section, "DELTA_CHARGE", r_val=delta)
222 2 : CALL section_vals_val_get(grb_section, "DERIVATIVES", i_val=nder)
223 2 : IF (iw > 0) THEN
224 2 : WRITE (iw, '(/,A,T76,I5)') " Generate Response Basis Sets with Order ", nder
225 : END IF
226 :
227 2 : state => atom%state
228 : ! find HOMO
229 2 : lhomo = -1
230 2 : nhomo = -1
231 2 : emax = -HUGE(1._dp)
232 6 : DO l = 0, state%maxl_occ
233 10 : DO i = 1, state%maxn_occ(l)
234 8 : IF (atom%orbitals%ener(i, l) > emax) THEN
235 4 : lhomo = l
236 4 : nhomo = i
237 4 : emax = atom%orbitals%ener(i, l)
238 4 : fhomo = state%occupation(l, i)
239 : END IF
240 : END DO
241 : END DO
242 :
243 2 : s1 = SIZE(atom%orbitals%wfn, 1)
244 2 : s2 = SIZE(atom%orbitals%wfn, 2)
245 12 : ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
246 14 : s2 = MAXVAL(state%maxn_occ) + nder
247 14 : ALLOCATE (rbasis(s1, s2, 0:lmat), qbasis(s1, s2, 0:lmat))
248 2 : rbasis = 0._dp
249 2 : qbasis = 0._dp
250 :
251 : ! calculate integrals
252 426 : ALLOCATE (atint)
253 2 : CALL atom_int_setup(atint, basis, potential=atom%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
254 2 : CALL atom_ppint_setup(atint, basis, potential=atom%potential)
255 2 : IF (atom%pp_calc) THEN
256 2 : NULLIFY (atint%tzora, atint%hdkh)
257 : ELSE
258 : ! relativistic correction terms
259 0 : CALL atom_relint_setup(atint, basis, atom%relativistic, zcore=REAL(atom%z, dp))
260 : END IF
261 2 : CALL set_atom(atom, integrals=atint)
262 :
263 2 : CALL calculate_atom(atom, iw=0)
264 16 : DO ider = -nder, nder
265 14 : dene = REAL(ider, KIND=dp)*delta
266 14 : CPASSERT(fhomo > ABS(dene))
267 14 : state%occupation(lhomo, nhomo) = fhomo + dene
268 14 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
269 686 : wfn(:, :, :, ider) = atom%orbitals%wfn
270 16 : state%occupation(lhomo, nhomo) = fhomo
271 : END DO
272 2 : IF (iw > 0) THEN
273 2 : WRITE (iw, '(A,T76,I5)') " Total number of electronic structure calculations ", 2*nder + 1
274 : END IF
275 :
276 2 : ovlp => atom%integrals%ovlp
277 :
278 6 : DO l = 0, state%maxl_occ
279 4 : IF (iw > 0) THEN
280 4 : WRITE (iw, '(A,T76,I5)') " Response derivatives for l quantum number ", l
281 : END IF
282 : ! occupied states
283 8 : DO i = 1, MAX(state%maxn_occ(l), 1)
284 32 : rbasis(:, i, l) = wfn(:, i, l, 0)
285 : END DO
286 : ! differentiation
287 16 : DO ider = 1, nder
288 12 : i = MAX(state%maxn_occ(l), 1)
289 4 : SELECT CASE (ider)
290 : CASE (1)
291 28 : rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
292 : CASE (2)
293 28 : 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 28 : + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
297 : CASE DEFAULT
298 12 : 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 4 : n = state%maxn_occ(l) + nder
304 4 : m = atom%basis%nbas(l)
305 20 : DO i = 1, n
306 40 : DO j = 1, i - 1
307 2136 : o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
308 184 : rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
309 : END DO
310 1424 : o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), [m]))
311 116 : rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
312 : END DO
313 :
314 : ! check
315 16 : ALLOCATE (amat(n, n))
316 2208 : 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 20 : DO i = 1, n
318 20 : amat(i, i) = amat(i, i) - 1._dp
319 : END DO
320 84 : IF (MAXVAL(ABS(amat)) > error_threshold) THEN
321 0 : IF (iw > 0) WRITE (iw, '(A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
322 : END IF
323 4 : DEALLOCATE (amat)
324 :
325 : ! Quickstep normalization
326 4 : expzet = 0.25_dp*REAL(2*l + 3, dp)
327 4 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
328 30 : DO i = 1, m
329 24 : zeta = (2._dp*atom%basis%am(i, l))**expzet
330 124 : 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 2 : IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Valence Response Basis Sets"
337 2 : CALL init_orbital_pointers(lmat)
338 2 : CALL init_spherical_harmonics(lmat, 0)
339 10 : DO ider = 0, nder
340 8 : NULLIFY (basis_vrb)
341 152 : 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 8 : CALL allocate_grid_atom(basis_vrb%grid)
347 8 : CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
348 8 : basis_vrb%grid%nr = ngp
349 : !
350 8 : basis_vrb%eps_eig = basis_ref%eps_eig
351 8 : basis_vrb%geometrical = .FALSE.
352 8 : basis_vrb%basis_type = CGTO_BASIS
353 104 : basis_vrb%nprim = basis%nprim
354 56 : basis_vrb%nbas = 0
355 24 : DO l = 0, state%maxl_occ
356 24 : basis_vrb%nbas(l) = state%maxn_occ(l) + ider
357 : END DO
358 56 : m = MAXVAL(basis_vrb%nprim)
359 56 : n = MAXVAL(basis_vrb%nbas)
360 24 : ALLOCATE (basis_vrb%am(m, 0:lmat))
361 680 : basis_vrb%am = basis%am
362 : ! contractions
363 40 : ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
364 24 : DO l = 0, state%maxl_occ
365 16 : m = basis_vrb%nprim(l)
366 16 : n = basis_vrb%nbas(l)
367 304 : 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 8 : nr = basis_vrb%grid%nr
372 56 : m = MAXVAL(basis_vrb%nbas)
373 40 : ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
374 24 : ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
375 24 : ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
376 48176 : basis_vrb%bf = 0._dp
377 48176 : basis_vrb%dbf = 0._dp
378 48176 : basis_vrb%ddbf = 0._dp
379 56 : DO l = 0, lmat
380 152 : DO i = 1, basis_vrb%nprim(l)
381 96 : al = basis_vrb%am(i, l)
382 38544 : DO k = 1, nr
383 38400 : rk = basis_vrb%grid%rad(k)
384 38400 : ear = EXP(-al*basis_vrb%grid%rad(k)**2)
385 134496 : DO j = 1, basis_vrb%nbas(l)
386 96000 : 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 96000 : + (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 134400 : 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 8 : IF (iw > 0) THEN
398 8 : CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
399 8 : WRITE (iw, '(A,A)') " Basis set ", TRIM(abas)
400 : END IF
401 8 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
402 8 : cradx = crad*1.00_dp
403 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
404 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
405 8 : cradx = crad*1.10_dp
406 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
407 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
408 8 : cradx = crad*1.20_dp
409 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
410 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
411 34 : vbasis(ider)%basis => basis_vrb
412 : END DO
413 2 : CALL deallocate_orbital_pointers
414 2 : CALL deallocate_spherical_harmonics
415 :
416 : ! get density maximum
417 6 : ALLOCATE (rho(basis%grid%nr))
418 2 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
419 2 : CALL atom_density(rho(:), atom%orbitals%pmat, atom%basis, maxl, typ="RHO")
420 804 : n = SUM(MAXLOC(rho(:)))
421 2 : rmax = basis%grid%rad(n)
422 2 : IF (rmax < 0.1_dp) rmax = 1.0_dp
423 2 : DEALLOCATE (rho)
424 :
425 : ! generate polarization sets
426 2 : maxl = atom%state%maxl_occ
427 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto)
428 2 : num_pol = num_gto
429 2 : IF (num_gto > 0) THEN
430 1 : IF (iw > 0) THEN
431 1 : WRITE (iw, '(/,A)') " Polarization basis set "
432 : END IF
433 7 : ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto))
434 1 : pbasis = 0.0_dp
435 : ! optimize exponents
436 1 : lval = maxl + 1
437 1 : zval = SQRT(REAL(2*lval + 2, dp))*REAL(lval + 1, dp)/(2._dp*rmax)
438 1 : aval = atom%basis%am(1, 0)
439 1 : cval = 2.5_dp
440 1 : rconf = atom%potential%scon
441 1 : CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section)
442 : ! calculate contractions
443 5 : DO i = 1, num_gto
444 5 : alp(i) = aval*cval**(i - 1)
445 : END DO
446 3 : ALLOCATE (rho(num_gto))
447 5 : DO l = maxl + 1, MIN(maxl + num_gto, 7)
448 4 : zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
449 4 : CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l))
450 4 : IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
451 5 : " Polarization basis set contraction for lval=", l, "zval=", zval
452 : END DO
453 1 : DEALLOCATE (rho)
454 : END IF
455 :
456 : ! generate valence expansion sets
457 2 : maxl = atom%state%maxl_occ
458 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_EXTENDED", i_val=num_gto)
459 2 : CALL section_vals_val_get(grb_section, "EXTENSION_BASIS", i_vals=num_bas)
460 2 : next_bas(0:lmat) = 0
461 2 : IF (num_bas(1) == -1) THEN
462 0 : DO l = 0, maxl
463 0 : next_bas(l) = maxl - l + 1
464 : END DO
465 : ELSE
466 2 : n = MIN(SIZE(num_bas, 1), 4)
467 6 : next_bas(0:n - 1) = num_bas(1:n)
468 : END IF
469 2 : next_prim = 0
470 14 : DO l = 0, lmat
471 14 : IF (next_bas(l) > 0) next_prim(l) = num_gto
472 : END DO
473 2 : IF (iw > 0) THEN
474 2 : CALL basis_label(abas, next_prim, next_bas)
475 2 : WRITE (iw, '(/,A,A)') " Extension basis set ", TRIM(abas)
476 : END IF
477 14 : n = MAXVAL(next_prim)
478 14 : m = MAXVAL(next_bas)
479 11 : ALLOCATE (ebasis(n, n, 0:lmat), ale(n))
480 2 : basis_vrb => vbasis(0)%basis
481 2 : amin = atom%basis%aval/atom%basis%cval**1.5_dp
482 6 : DO i = 1, n
483 6 : ale(i) = amin*atom%basis%cval**(i - 1)
484 : END DO
485 2 : ebasis = 0._dp
486 5 : ALLOCATE (rho(n))
487 2 : rconf = 2.0_dp*atom%potential%scon
488 14 : DO l = 0, lmat
489 12 : IF (next_bas(l) < 1) CYCLE
490 2 : zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
491 2 : CALL hydrogenic(zval, rconf, l, ale, n, rho, ebasis(:, :, l))
492 2 : IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
493 4 : " Extension basis set contraction for lval=", l, "zval=", zval
494 : END DO
495 2 : DEALLOCATE (rho)
496 : ! check for condition numbers
497 2 : IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Extended Basis Sets"
498 2 : CALL init_orbital_pointers(lmat)
499 2 : CALL init_spherical_harmonics(lmat, 0)
500 10 : DO ider = 0, nder
501 8 : NULLIFY (basis_vrb)
502 152 : 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 8 : CALL allocate_grid_atom(basis_vrb%grid)
508 8 : CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
509 8 : basis_vrb%grid%nr = ngp
510 : !
511 8 : basis_vrb%eps_eig = basis_ref%eps_eig
512 8 : basis_vrb%geometrical = .FALSE.
513 8 : basis_vrb%basis_type = CGTO_BASIS
514 104 : basis_vrb%nprim = basis%nprim + next_prim
515 56 : basis_vrb%nbas = 0
516 24 : DO l = 0, state%maxl_occ
517 24 : basis_vrb%nbas(l) = state%maxn_occ(l) + ider + next_bas(l)
518 : END DO
519 56 : m = MAXVAL(basis_vrb%nprim)
520 24 : ALLOCATE (basis_vrb%am(m, 0:lmat))
521 : ! exponents
522 8 : m = SIZE(basis%am, 1)
523 680 : basis_vrb%am(1:m, :) = basis%am(1:m, :)
524 8 : n = SIZE(ale, 1)
525 24 : DO l = 0, state%maxl_occ
526 56 : basis_vrb%am(m + 1:m + n, l) = ale(1:n)
527 : END DO
528 : ! contractions
529 56 : m = MAXVAL(basis_vrb%nprim)
530 56 : n = MAXVAL(basis_vrb%nbas)
531 40 : ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
532 1664 : basis_vrb%cm = 0.0_dp
533 24 : DO l = 0, state%maxl_occ
534 16 : m = basis%nprim(l)
535 16 : n = state%maxn_occ(l) + ider
536 296 : basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
537 84 : 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 8 : nr = basis_vrb%grid%nr
542 56 : m = MAXVAL(basis_vrb%nbas)
543 40 : ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
544 24 : ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
545 24 : ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
546 67424 : basis_vrb%bf = 0._dp
547 67424 : basis_vrb%dbf = 0._dp
548 67424 : basis_vrb%ddbf = 0._dp
549 56 : DO l = 0, lmat
550 184 : DO i = 1, basis_vrb%nprim(l)
551 128 : al = basis_vrb%am(i, l)
552 51376 : DO k = 1, nr
553 51200 : rk = basis_vrb%grid%rad(k)
554 51200 : ear = EXP(-al*basis_vrb%grid%rad(k)**2)
555 227328 : DO j = 1, basis_vrb%nbas(l)
556 176000 : 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 176000 : + (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 227200 : 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 8 : IF (iw > 0) THEN
568 8 : CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
569 8 : WRITE (iw, '(A,A)') " Basis set ", TRIM(abas)
570 : END IF
571 8 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
572 8 : cradx = crad*1.00_dp
573 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
574 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
575 8 : cradx = crad*1.10_dp
576 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
577 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
578 8 : cradx = crad*1.20_dp
579 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
580 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
581 34 : vbasis(nder + 1 + ider)%basis => basis_vrb
582 : END DO
583 2 : CALL deallocate_orbital_pointers
584 2 : CALL deallocate_spherical_harmonics
585 :
586 : ! Tests for energy
587 2 : energy_ref = atom_ref%energy%etot
588 2 : IF (iw > 0) WRITE (iw, '(/,A,A)') " Basis set tests "
589 2 : IF (iw > 0) WRITE (iw, '(T10,A,T59,F22.9)') " Reference Energy [a.u.] ", energy_ref
590 18 : DO ider = 0, 2*nder + 1
591 : ! generate an atom type
592 16 : NULLIFY (atom_test)
593 16 : CALL create_atom_type(atom_test)
594 : CALL copy_atom_basics(atom_ref, atom_test, state=.TRUE., potential=.TRUE., &
595 16 : optimization=.TRUE., xc=.TRUE.)
596 16 : basis_grb => vbasis(ider)%basis
597 16 : NULLIFY (orbitals)
598 112 : mo = MAXVAL(atom_test%state%maxn_calc)
599 112 : mb = MAXVAL(basis_grb%nbas)
600 16 : CALL create_atom_orbs(orbitals, mb, mo)
601 16 : CALL set_atom(atom_test, orbitals=orbitals, basis=basis_grb)
602 : ! calculate integrals
603 3408 : ALLOCATE (atint)
604 16 : CALL atom_int_setup(atint, basis_grb, potential=atom_test%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
605 16 : CALL atom_ppint_setup(atint, basis_grb, potential=atom_test%potential)
606 16 : IF (atom_test%pp_calc) THEN
607 16 : NULLIFY (atint%tzora, atint%hdkh)
608 : ELSE
609 : ! relativistic correction terms
610 0 : CALL atom_relint_setup(atint, basis_grb, atom_test%relativistic, zcore=REAL(atom_test%z, dp))
611 : END IF
612 16 : CALL set_atom(atom_test, integrals=atint)
613 : !
614 16 : CALL calculate_atom(atom_test, iw=0)
615 16 : IF (ider <= nder) THEN
616 8 : energy_vb(ider) = atom_test%energy%etot
617 16 : IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (VB)", ider, " Energy [a.u.] ", &
618 16 : energy_ref - energy_vb(ider), energy_vb(ider)
619 : ELSE
620 8 : i = ider - nder - 1
621 8 : energy_ex(i) = atom_test%energy%etot
622 16 : IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (EX)", i, " Energy [a.u.] ", &
623 16 : energy_ref - energy_ex(i), energy_ex(i)
624 : END IF
625 16 : CALL atom_int_release(atint)
626 16 : CALL atom_ppint_release(atint)
627 16 : CALL atom_relint_release(atint)
628 16 : DEALLOCATE (atom_test%state, atom_test%potential, atint)
629 18 : CALL release_atom_type(atom_test)
630 : END DO
631 :
632 : ! Quickstep normalization polarization basis
633 18 : DO l = 0, 7
634 16 : expzet = 0.25_dp*REAL(2*l + 3, dp)
635 16 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
636 50 : DO i = 1, num_pol
637 32 : zeta = (2._dp*alp(i))**expzet
638 176 : 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 14 : DO l = 0, lmat
643 12 : expzet = 0.25_dp*REAL(2*l + 3, dp)
644 12 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
645 22 : DO i = 1, next_prim(l)
646 8 : zeta = (2._dp*ale(i))**expzet
647 32 : 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 2 : CALL section_vals_val_get(grb_section, "NAME_BODY", c_val=basname)
653 2 : CALL open_file(file_name="GRB_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iunit)
654 : ! header info
655 8 : headline = ""
656 2 : headline(1) = "#"
657 2 : headline(2) = "# Generated with CP2K Atom Code"
658 2 : headline(3) = "#"
659 2 : CALL grb_print_basis(header=headline, iunit=iunit)
660 : ! valence basis
661 2 : basline(1) = ""
662 2 : WRITE (basline(1), "(T2,A)") ADJUSTL(ptable(atom_ref%z)%symbol)
663 10 : DO ider = 0, nder
664 8 : basline(1) = ""
665 8 : 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 10 : al=vbasis(ider)%basis%am(:, 0), gcc=qbasis, iunit=iunit)
668 : END DO
669 : ! polarization basis
670 2 : maxl = atom_ref%state%maxl_occ
671 6 : DO l = maxl + 1, MIN(maxl + num_pol, 7)
672 4 : nbas = 0
673 14 : DO i = maxl + 1, l
674 14 : nbas(i) = l - i + 1
675 : END DO
676 4 : i = l - maxl
677 4 : basline(1) = ""
678 4 : WRITE (basline(1), "(T2,A,T5,A,I1)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-POL-", i
679 6 : CALL grb_print_basis(header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit)
680 : END DO
681 : ! extension set
682 14 : IF (SUM(next_bas) > 0) THEN
683 1 : basline(1) = ""
684 1 : WRITE (basline(1), "(T2,A,T5,A)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-EXT"
685 1 : CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit)
686 : END IF
687 : !
688 2 : CALL close_file(unit_number=iunit)
689 :
690 : ! clean up
691 2 : IF (ALLOCATED(pbasis)) DEALLOCATE (pbasis)
692 2 : IF (ALLOCATED(alp)) DEALLOCATE (alp)
693 2 : IF (ALLOCATED(ebasis)) DEALLOCATE (ebasis)
694 2 : DEALLOCATE (wfn, rbasis, qbasis, ale)
695 :
696 24 : DO ider = 0, 10
697 24 : IF (ASSOCIATED(vbasis(ider)%basis)) THEN
698 16 : CALL release_atom_basis(vbasis(ider)%basis)
699 16 : DEALLOCATE (vbasis(ider)%basis)
700 : END IF
701 : END DO
702 :
703 2 : CALL atom_int_release(atom%integrals)
704 2 : CALL atom_ppint_release(atom%integrals)
705 2 : CALL atom_relint_release(atom%integrals)
706 2 : CALL release_atom_basis(basis)
707 2 : DEALLOCATE (atom%potential, atom%state, atom%integrals, basis)
708 2 : CALL release_atom_type(atom)
709 :
710 2 : IF (iw > 0) WRITE (iw, '(" ",79("*"))')
711 :
712 18 : 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 15 : 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 15 : IF (PRESENT(header)) THEN
739 34 : DO i = 1, SIZE(header, 1)
740 34 : IF (header(i) /= "") THEN
741 19 : WRITE (iunit, "(A)") TRIM(header(i))
742 : END IF
743 : END DO
744 : END IF
745 :
746 15 : IF (PRESENT(nprim)) THEN
747 13 : IF (nprim > 0) THEN
748 13 : CPASSERT(PRESENT(nbas))
749 13 : CPASSERT(PRESENT(al))
750 13 : CPASSERT(PRESENT(gcc))
751 :
752 34 : DO i = LBOUND(nbas, 1), UBOUND(nbas, 1)
753 21 : IF (nbas(i) > 0) THEN
754 13 : lmin = i
755 13 : EXIT
756 : END IF
757 : END DO
758 63 : DO i = UBOUND(nbas, 1), LBOUND(nbas, 1), -1
759 63 : IF (nbas(i) > 0) THEN
760 13 : lmax = i
761 13 : EXIT
762 : END IF
763 : END DO
764 :
765 13 : nval = lmax
766 13 : WRITE (iunit, *) " 1"
767 13 : WRITE (iunit, "(40I3)") nval, lmin, lmax, nprim, (nbas(l), l=lmin, lmax)
768 81 : DO i = nprim, 1, -1
769 68 : WRITE (iunit, "(G20.12)", advance="no") al(i)
770 212 : DO l = lmin, lmax
771 544 : DO j = 1, nbas(l)
772 476 : WRITE (iunit, "(F16.10)", advance="no") gcc(i, j, l)
773 : END DO
774 : END DO
775 81 : WRITE (iunit, *)
776 : END DO
777 13 : WRITE (iunit, *)
778 : END IF
779 : END IF
780 :
781 15 : 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 18 : 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 18 : label = ""
800 18 : lmax = MIN(UBOUND(np, 1), UBOUND(nb, 1), 7)
801 18 : i = 1
802 18 : label(i:i) = "("
803 126 : DO l = 0, lmax
804 126 : IF (np(l) > 0) THEN
805 34 : i = i + 1
806 34 : IF (np(l) > 9) THEN
807 8 : WRITE (label(i:i + 1), "(I2)") np(l)
808 8 : i = i + 2
809 : ELSE
810 26 : WRITE (label(i:i), "(I1)") np(l)
811 26 : i = i + 1
812 : END IF
813 34 : label(i:i) = lq(l)
814 : END IF
815 : END DO
816 18 : i = i + 1
817 18 : label(i:i + 6) = ") --> ["
818 18 : i = i + 6
819 126 : DO l = 0, lmax
820 126 : IF (nb(l) > 0) THEN
821 34 : i = i + 1
822 34 : IF (nb(l) > 9) THEN
823 0 : WRITE (label(i:i + 1), "(I2)") nb(l)
824 0 : i = i + 2
825 : ELSE
826 34 : WRITE (label(i:i), "(I1)") nb(l)
827 34 : i = i + 1
828 : END IF
829 34 : label(i:i) = lq(l)
830 : END IF
831 : END DO
832 18 : i = i + 1
833 18 : label(i:i) = "]"
834 :
835 18 : 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 170 : 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 36890 : ALLOCATE (atint)
858 : ! calculate integrals
859 170 : NULLIFY (pot)
860 170 : eri_c = .FALSE.
861 170 : eri_e = .FALSE.
862 170 : pot => atom%potential
863 170 : zval = atom%z
864 170 : reltyp = atom%relativistic
865 170 : do_eric = atom%coulomb_integral_type
866 170 : do_erie = atom%exchange_integral_type
867 170 : IF (do_eric == do_analytic) eri_c = .TRUE.
868 170 : IF (do_erie == do_analytic) eri_e = .TRUE.
869 : ! general integrals
870 170 : CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
871 : ! potential
872 170 : CALL atom_ppint_setup(atint, basis, potential=pot)
873 170 : IF (atom%pp_calc) THEN
874 170 : NULLIFY (atint%tzora, atint%hdkh)
875 : ELSE
876 : ! relativistic correction terms
877 0 : CALL atom_relint_setup(atint, basis, reltyp, zcore=REAL(zval, dp))
878 : END IF
879 170 : CALL set_atom(atom, basis=basis)
880 170 : CALL set_atom(atom, integrals=atint)
881 170 : CALL calculate_atom(atom, iw)
882 170 : afun = atom%energy%etot
883 170 : CALL atom_int_release(atint)
884 170 : CALL atom_ppint_release(atint)
885 170 : CALL atom_relint_release(atint)
886 170 : DEALLOCATE (atint)
887 170 : 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 18 : 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 18 : atom%z = atom_ref%z
905 18 : atom%zcore = atom_ref%zcore
906 18 : atom%pp_calc = atom_ref%pp_calc
907 18 : atom%method_type = atom_ref%method_type
908 18 : atom%relativistic = atom_ref%relativistic
909 18 : atom%coulomb_integral_type = atom_ref%coulomb_integral_type
910 18 : atom%exchange_integral_type = atom_ref%exchange_integral_type
911 :
912 18 : NULLIFY (atom%potential, atom%state, atom%xc_section)
913 18 : NULLIFY (atom%basis, atom%integrals, atom%orbitals, atom%fmat)
914 :
915 18 : IF (PRESENT(state)) THEN
916 18 : IF (state) THEN
917 6660 : ALLOCATE (atom%state)
918 18 : atom%state = atom_ref%state
919 : END IF
920 : END IF
921 :
922 18 : IF (PRESENT(potential)) THEN
923 18 : IF (potential) THEN
924 97182 : ALLOCATE (atom%potential)
925 18 : atom%potential = atom_ref%potential
926 : END IF
927 : END IF
928 :
929 18 : IF (PRESENT(optimization)) THEN
930 18 : IF (optimization) THEN
931 18 : atom%optimization = atom_ref%optimization
932 : END IF
933 : END IF
934 :
935 18 : IF (PRESENT(xc)) THEN
936 18 : IF (xc) THEN
937 18 : atom%xc_section => atom_ref%xc_section
938 : END IF
939 : END IF
940 :
941 18 : 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 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x
961 : TYPE(opt_state_type) :: ostate
962 :
963 0 : CPASSERT(basis%geometrical)
964 :
965 2 : CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
966 2 : CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
967 2 : CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
968 :
969 2 : ostate%nvar = 2
970 2 : ALLOCATE (x(2))
971 2 : x(1) = SQRT(basis%aval)
972 2 : x(2) = SQRT(basis%cval)
973 :
974 2 : ostate%nf = 0
975 2 : ostate%iprint = 1
976 2 : ostate%unit = iunit
977 :
978 2 : ostate%state = 0
979 2 : IF (iunit > 0) THEN
980 2 : WRITE (iunit, '(/," POWELL| Start optimization procedure")')
981 2 : WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
982 : END IF
983 2 : n10 = MAX(ostate%maxfun/100, 1)
984 :
985 2 : fopt = HUGE(0._dp)
986 :
987 : DO
988 :
989 176 : IF (ostate%state == 2) THEN
990 7310 : basis%am = 0._dp
991 1190 : DO l = 0, lmat
992 3230 : DO i = 1, basis%nbas(l)
993 2040 : ll = i - 1 + basis%start(l)
994 3060 : basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
995 : END DO
996 : END DO
997 170 : basis%aval = x(1)*x(1)
998 170 : basis%cval = x(2)*x(2)
999 2455310 : basis%bf = 0._dp
1000 2455310 : basis%dbf = 0._dp
1001 2455310 : basis%ddbf = 0._dp
1002 170 : nr = basis%grid%nr
1003 1190 : DO l = 0, lmat
1004 3230 : DO i = 1, basis%nbas(l)
1005 2040 : al = basis%am(i, l)
1006 819060 : DO k = 1, nr
1007 816000 : rk = basis%grid%rad(k)
1008 816000 : ear = EXP(-al*basis%grid%rad(k)**2)
1009 816000 : basis%bf(k, i, l) = rk**l*ear
1010 816000 : 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 818040 : 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 170 : CALL grb_fit(atom, basis, ostate%f, 0)
1017 170 : fopt = MIN(fopt, ostate%f)
1018 : END IF
1019 :
1020 176 : IF (ostate%state == -1) EXIT
1021 :
1022 174 : CALL powell_optimize(ostate%nvar, x, ostate)
1023 :
1024 174 : IF (ostate%nf == 2 .AND. iunit > 0) THEN
1025 2 : WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1026 : END IF
1027 176 : IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1028 : WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1029 2 : INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
1030 : END IF
1031 :
1032 : END DO
1033 :
1034 2 : ostate%state = 8
1035 2 : CALL powell_optimize(ostate%nvar, x, ostate)
1036 :
1037 2 : IF (iunit > 0) THEN
1038 2 : WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1039 2 : WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1040 : END IF
1041 : ! x->basis
1042 86 : basis%am = 0._dp
1043 14 : DO l = 0, lmat
1044 38 : DO i = 1, basis%nbas(l)
1045 24 : ll = i - 1 + basis%start(l)
1046 36 : basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
1047 : END DO
1048 : END DO
1049 2 : basis%aval = x(1)*x(1)
1050 2 : basis%cval = x(2)*x(2)
1051 28886 : basis%bf = 0._dp
1052 28886 : basis%dbf = 0._dp
1053 28886 : basis%ddbf = 0._dp
1054 2 : nr = basis%grid%nr
1055 14 : DO l = 0, lmat
1056 38 : DO i = 1, basis%nbas(l)
1057 24 : al = basis%am(i, l)
1058 9636 : DO k = 1, nr
1059 9600 : rk = basis%grid%rad(k)
1060 9600 : ear = EXP(-al*basis%grid%rad(k)**2)
1061 9600 : basis%bf(k, i, l) = rk**l*ear
1062 9600 : 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 9624 : 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 2 : DEALLOCATE (x)
1070 :
1071 : ! final result
1072 2 : IF (iunit > 0) THEN
1073 2 : WRITE (iunit, '(/,A)') " Optimized Geometrical GTO basis set"
1074 2 : WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", basis%aval, &
1075 4 : " Proportionality factor: ", basis%cval
1076 14 : DO l = 0, lmat
1077 14 : WRITE (iunit, '(T41,A,I2,T76,I5)') " Number of exponents for l=", l, basis%nbas(l)
1078 : END DO
1079 : END IF
1080 :
1081 2 : IF (iunit > 0) WRITE (iunit, '(/,A)') " Condition number of uncontracted basis set"
1082 2 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
1083 2 : CALL init_orbital_pointers(lmat)
1084 2 : CALL init_spherical_harmonics(lmat, 0)
1085 2 : cradx = crad*1.00_dp
1086 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1087 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1088 2 : cradx = crad*1.10_dp
1089 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1090 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1091 2 : cradx = crad*1.20_dp
1092 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1093 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1094 2 : CALL deallocate_orbital_pointers
1095 2 : CALL deallocate_spherical_harmonics
1096 :
1097 8 : 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 1 : 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 7 : ALLOCATE (am(nbas), ener(nbas), orb(nbas, nbas))
1126 :
1127 1 : CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
1128 1 : CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
1129 1 : CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
1130 :
1131 1 : ostate%nvar = 2
1132 1 : x(1) = SQRT(aval)
1133 1 : x(2) = SQRT(cval)
1134 :
1135 1 : ostate%nf = 0
1136 1 : ostate%iprint = 1
1137 1 : ostate%unit = iunit
1138 :
1139 1 : ostate%state = 0
1140 1 : IF (iunit > 0) THEN
1141 1 : WRITE (iunit, '(/," POWELL| Start optimization procedure")')
1142 1 : WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
1143 : END IF
1144 1 : n10 = MAX(ostate%maxfun/100, 1)
1145 :
1146 1 : fopt = HUGE(0._dp)
1147 :
1148 : DO
1149 :
1150 81 : IF (ostate%state == 2) THEN
1151 78 : aval = x(1)*x(1)
1152 78 : cval = x(2)*x(2)
1153 390 : DO i = 1, nbas
1154 390 : am(i) = aval*cval**(i - 1)
1155 : END DO
1156 78 : CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1157 78 : ostate%f = ener(1)
1158 78 : fopt = MIN(fopt, ostate%f)
1159 : END IF
1160 :
1161 81 : IF (ostate%state == -1) EXIT
1162 :
1163 80 : CALL powell_optimize(ostate%nvar, x, ostate)
1164 :
1165 80 : IF (ostate%nf == 2 .AND. iunit > 0) THEN
1166 1 : WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1167 : END IF
1168 81 : IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1169 : WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1170 1 : INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
1171 : END IF
1172 :
1173 : END DO
1174 :
1175 1 : ostate%state = 8
1176 1 : CALL powell_optimize(ostate%nvar, x, ostate)
1177 :
1178 1 : IF (iunit > 0) THEN
1179 1 : WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1180 1 : WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1181 : END IF
1182 : ! x->basis
1183 1 : aval = x(1)*x(1)
1184 1 : cval = x(2)*x(2)
1185 :
1186 : ! final result
1187 1 : IF (iunit > 0) THEN
1188 1 : WRITE (iunit, '(/,A,T51,A,T76,I5)') " Optimized Polarization basis set", &
1189 2 : " Number of exponents:", nbas
1190 1 : WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", aval, &
1191 2 : " Proportionality factor: ", cval
1192 : END IF
1193 :
1194 1 : DEALLOCATE (am, ener, orb)
1195 :
1196 1 : 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 84 : 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 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
1221 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: confmat, hmat, potmat, smat, tmat
1222 :
1223 84 : n = nbas
1224 1008 : ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
1225 : ! calclulate overlap matrix
1226 84 : CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
1227 : ! calclulate kinetic energy matrix
1228 84 : CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
1229 : ! calclulate core potential matrix
1230 84 : CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
1231 : ! calclulate confinement potential matrix
1232 84 : cf = 0.1_dp
1233 84 : k = 10
1234 84 : CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
1235 : ! Hamiltionian
1236 1764 : 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 84 : lwork = 100*n
1239 420 : ALLOCATE (w(n), work(lwork))
1240 84 : CALL dsygv(1, "V", "U", n, hmat, n, smat, n, w, work, lwork, info)
1241 84 : CPASSERT(info == 0)
1242 1764 : orb(1:n, 1:n) = hmat(1:n, 1:n)
1243 420 : ener(1:n) = w(1:n)
1244 84 : DEALLOCATE (w, work)
1245 84 : DEALLOCATE (smat, tmat, potmat, confmat, hmat)
1246 :
1247 84 : END SUBROUTINE hydrogenic
1248 :
1249 44 : END MODULE atom_grb
|
