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