Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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("")
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) &
156 0 : CPABORT("# point radial grid < 0")
157 2 : CALL create_grid_atom(basis%grid, ngp, 1, 1, 0, quadtype)
158 2 : basis%grid%nr = ngp
159 : !
160 2 : maxl = atom%state%maxl_occ
161 2 : basis%basis_type = GTO_BASIS
162 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_CORE", i_val=num_gto)
163 14 : basis%nbas = 0
164 6 : basis%nbas(0:maxl) = num_gto
165 14 : basis%nprim = basis%nbas
166 2 : CALL section_vals_val_get(grb_section, "GEOMETRICAL_FACTOR", r_val=cval)
167 2 : CALL section_vals_val_get(grb_section, "GEO_START_VALUE", r_val=aval)
168 14 : m = MAXVAL(basis%nbas)
169 6 : ALLOCATE (basis%am(m, 0:lmat))
170 86 : basis%am = 0._dp
171 14 : DO l = 0, lmat
172 38 : DO i = 1, basis%nbas(l)
173 24 : ll = i - 1
174 36 : basis%am(i, l) = aval*cval**(ll)
175 : END DO
176 : END DO
177 :
178 2 : basis%eps_eig = basis_ref%eps_eig
179 2 : basis%geometrical = .TRUE.
180 2 : basis%aval = aval
181 2 : basis%cval = cval
182 14 : basis%start = 0
183 :
184 : ! initialize basis function on a radial grid
185 2 : nr = basis%grid%nr
186 14 : m = MAXVAL(basis%nbas)
187 10 : ALLOCATE (basis%bf(nr, m, 0:lmat))
188 6 : ALLOCATE (basis%dbf(nr, m, 0:lmat))
189 6 : ALLOCATE (basis%ddbf(nr, m, 0:lmat))
190 28886 : basis%bf = 0._dp
191 28886 : basis%dbf = 0._dp
192 28886 : basis%ddbf = 0._dp
193 14 : DO l = 0, lmat
194 38 : DO i = 1, basis%nbas(l)
195 24 : al = basis%am(i, l)
196 9636 : DO k = 1, nr
197 9600 : rk = basis%grid%rad(k)
198 9600 : ear = EXP(-al*basis%grid%rad(k)**2)
199 9600 : basis%bf(k, i, l) = rk**l*ear
200 9600 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
201 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
202 9624 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
203 : END DO
204 : END DO
205 : END DO
206 :
207 2 : NULLIFY (orbitals)
208 14 : mo = MAXVAL(atom%state%maxn_calc)
209 14 : mb = MAXVAL(basis%nbas)
210 2 : CALL create_atom_orbs(orbitals, mb, mo)
211 2 : CALL set_atom(atom, orbitals=orbitals)
212 :
213 2 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
214 2 : CALL atom_fit_grb(atom, basis, iw, powell_section)
215 2 : CALL set_atom(atom, basis=basis)
216 :
217 : ! generate response contractions
218 2 : CALL section_vals_val_get(grb_section, "DELTA_CHARGE", r_val=delta)
219 2 : CALL section_vals_val_get(grb_section, "DERIVATIVES", i_val=nder)
220 2 : IF (iw > 0) THEN
221 2 : WRITE (iw, '(/,A,T76,I5)') " Generate Response Basis Sets with Order ", nder
222 : END IF
223 :
224 2 : state => atom%state
225 : ! find HOMO
226 2 : lhomo = -1
227 2 : nhomo = -1
228 2 : emax = -HUGE(1._dp)
229 6 : DO l = 0, state%maxl_occ
230 10 : DO i = 1, state%maxn_occ(l)
231 8 : IF (atom%orbitals%ener(i, l) > emax) THEN
232 4 : lhomo = l
233 4 : nhomo = i
234 4 : emax = atom%orbitals%ener(i, l)
235 4 : fhomo = state%occupation(l, i)
236 : END IF
237 : END DO
238 : END DO
239 :
240 2 : s1 = SIZE(atom%orbitals%wfn, 1)
241 2 : s2 = SIZE(atom%orbitals%wfn, 2)
242 12 : ALLOCATE (wfn(s1, s2, 0:lmat, -nder:nder))
243 14 : s2 = MAXVAL(state%maxn_occ) + nder
244 14 : ALLOCATE (rbasis(s1, s2, 0:lmat), qbasis(s1, s2, 0:lmat))
245 350 : rbasis = 0._dp
246 350 : qbasis = 0._dp
247 :
248 : ! calculate integrals
249 426 : ALLOCATE (atint)
250 2 : CALL atom_int_setup(atint, basis, potential=atom%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
251 2 : CALL atom_ppint_setup(atint, basis, potential=atom%potential)
252 2 : IF (atom%pp_calc) THEN
253 2 : NULLIFY (atint%tzora, atint%hdkh)
254 : ELSE
255 : ! relativistic correction terms
256 0 : CALL atom_relint_setup(atint, basis, atom%relativistic, zcore=REAL(atom%z, dp))
257 : END IF
258 2 : CALL set_atom(atom, integrals=atint)
259 :
260 2 : CALL calculate_atom(atom, iw=0)
261 16 : DO ider = -nder, nder
262 14 : dene = REAL(ider, KIND=dp)*delta
263 14 : CPASSERT(fhomo > ABS(dene))
264 14 : state%occupation(lhomo, nhomo) = fhomo + dene
265 14 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
266 686 : wfn(:, :, :, ider) = atom%orbitals%wfn
267 16 : state%occupation(lhomo, nhomo) = fhomo
268 : END DO
269 2 : IF (iw > 0) THEN
270 2 : WRITE (iw, '(A,T76,I5)') " Total number of electronic structure calculations ", 2*nder + 1
271 : END IF
272 :
273 2 : ovlp => atom%integrals%ovlp
274 :
275 6 : DO l = 0, state%maxl_occ
276 4 : IF (iw > 0) THEN
277 4 : WRITE (iw, '(A,T76,I5)') " Response derivatives for l quantum number ", l
278 : END IF
279 : ! occupied states
280 8 : DO i = 1, MAX(state%maxn_occ(l), 1)
281 32 : rbasis(:, i, l) = wfn(:, i, l, 0)
282 : END DO
283 : ! differentiation
284 16 : DO ider = 1, nder
285 12 : i = MAX(state%maxn_occ(l), 1)
286 4 : SELECT CASE (ider)
287 : CASE (1)
288 28 : rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
289 : CASE (2)
290 28 : rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
291 : CASE (3)
292 : rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
293 28 : + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
294 : CASE DEFAULT
295 12 : CPABORT("")
296 : END SELECT
297 : END DO
298 :
299 : ! orthogonalization, use gram-schmidt in order to keep the natural order (semi-core, valence, response) of the wfn.
300 4 : n = state%maxn_occ(l) + nder
301 4 : m = atom%basis%nbas(l)
302 20 : DO i = 1, n
303 40 : DO j = 1, i - 1
304 2304 : o = DOT_PRODUCT(rbasis(1:m, j, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
305 184 : rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
306 : END DO
307 1536 : o = DOT_PRODUCT(rbasis(1:m, i, l), RESHAPE(MATMUL(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
308 116 : rbasis(1:m, i, l) = rbasis(1:m, i, l)/SQRT(o)
309 : END DO
310 :
311 : ! check
312 16 : ALLOCATE (amat(n, n))
313 2320 : 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)))
314 20 : DO i = 1, n
315 20 : amat(i, i) = amat(i, i) - 1._dp
316 : END DO
317 84 : IF (MAXVAL(ABS(amat)) > 1.e-12) THEN
318 0 : IF (iw > 0) WRITE (iw, '(A,G20.10)') " Orthogonality error ", MAXVAL(ABS(amat))
319 : END IF
320 4 : DEALLOCATE (amat)
321 :
322 : ! Quickstep normalization
323 4 : expzet = 0.25_dp*REAL(2*l + 3, dp)
324 4 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
325 30 : DO i = 1, m
326 24 : zeta = (2._dp*atom%basis%am(i, l))**expzet
327 124 : qbasis(i, 1:n, l) = rbasis(i, 1:n, l)*prefac/zeta
328 : END DO
329 :
330 : END DO
331 :
332 : ! check for condition numbers
333 2 : IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Valence Response Basis Sets"
334 2 : CALL init_orbital_pointers(lmat)
335 2 : CALL init_spherical_harmonics(lmat, 0)
336 10 : DO ider = 0, nder
337 8 : NULLIFY (basis_vrb)
338 152 : ALLOCATE (basis_vrb)
339 : NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
340 : basis_vrb%dbf, basis_vrb%ddbf)
341 : ! allocate and initialize the atomic grid
342 : NULLIFY (basis_vrb%grid)
343 8 : CALL allocate_grid_atom(basis_vrb%grid)
344 8 : CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
345 8 : basis_vrb%grid%nr = ngp
346 : !
347 8 : basis_vrb%eps_eig = basis_ref%eps_eig
348 8 : basis_vrb%geometrical = .FALSE.
349 8 : basis_vrb%basis_type = CGTO_BASIS
350 104 : basis_vrb%nprim = basis%nprim
351 56 : basis_vrb%nbas = 0
352 24 : DO l = 0, state%maxl_occ
353 24 : basis_vrb%nbas(l) = state%maxn_occ(l) + ider
354 : END DO
355 56 : m = MAXVAL(basis_vrb%nprim)
356 56 : n = MAXVAL(basis_vrb%nbas)
357 24 : ALLOCATE (basis_vrb%am(m, 0:lmat))
358 680 : basis_vrb%am = basis%am
359 : ! contractions
360 40 : ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
361 24 : DO l = 0, state%maxl_occ
362 16 : m = basis_vrb%nprim(l)
363 16 : n = basis_vrb%nbas(l)
364 304 : basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
365 : END DO
366 :
367 : ! initialize basis function on a radial grid
368 8 : nr = basis_vrb%grid%nr
369 56 : m = MAXVAL(basis_vrb%nbas)
370 40 : ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
371 24 : ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
372 24 : ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
373 48176 : basis_vrb%bf = 0._dp
374 48176 : basis_vrb%dbf = 0._dp
375 48176 : basis_vrb%ddbf = 0._dp
376 56 : DO l = 0, lmat
377 152 : DO i = 1, basis_vrb%nprim(l)
378 96 : al = basis_vrb%am(i, l)
379 38544 : DO k = 1, nr
380 38400 : rk = basis_vrb%grid%rad(k)
381 38400 : ear = EXP(-al*basis_vrb%grid%rad(k)**2)
382 134496 : DO j = 1, basis_vrb%nbas(l)
383 96000 : basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
384 : basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
385 96000 : + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
386 : basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
387 : (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
388 134400 : 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
389 : END DO
390 : END DO
391 : END DO
392 : END DO
393 :
394 8 : IF (iw > 0) THEN
395 8 : CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
396 8 : WRITE (iw, '(A,A)') " Basis set ", TRIM(abas)
397 : END IF
398 8 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
399 8 : cradx = crad*1.00_dp
400 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
401 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
402 8 : cradx = crad*1.10_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.20_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 34 : vbasis(ider)%basis => basis_vrb
409 : END DO
410 2 : CALL deallocate_orbital_pointers
411 2 : CALL deallocate_spherical_harmonics
412 :
413 : ! get density maximum
414 6 : ALLOCATE (rho(basis%grid%nr))
415 2 : CALL calculate_atom(atom, iw=0, noguess=.TRUE.)
416 2 : CALL atom_density(rho(:), atom%orbitals%pmat, atom%basis, maxl, typ="RHO")
417 802 : n = SUM(MAXLOC(rho(:)))
418 2 : rmax = basis%grid%rad(n)
419 2 : IF (rmax < 0.1_dp) rmax = 1.0_dp
420 2 : DEALLOCATE (rho)
421 :
422 : ! generate polarization sets
423 2 : maxl = atom%state%maxl_occ
424 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_POLARIZATION", i_val=num_gto)
425 2 : num_pol = num_gto
426 2 : IF (num_gto > 0) THEN
427 1 : IF (iw > 0) THEN
428 1 : WRITE (iw, '(/,A)') " Polarization basis set "
429 : END IF
430 7 : ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto))
431 169 : pbasis = 0.0_dp
432 : ! optimize exponents
433 1 : lval = maxl + 1
434 1 : zval = SQRT(REAL(2*lval + 2, dp))*REAL(lval + 1, dp)/(2._dp*rmax)
435 1 : aval = atom%basis%am(1, 0)
436 1 : cval = 2.5_dp
437 1 : rconf = atom%potential%scon
438 1 : CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section)
439 : ! calculate contractions
440 5 : DO i = 1, num_gto
441 5 : alp(i) = aval*cval**(i - 1)
442 : END DO
443 2 : ALLOCATE (rho(num_gto))
444 5 : DO l = maxl + 1, MIN(maxl + num_gto, 7)
445 4 : zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
446 4 : CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l))
447 4 : IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
448 5 : " Polarization basis set contraction for lval=", l, "zval=", zval
449 : END DO
450 1 : DEALLOCATE (rho)
451 : END IF
452 :
453 : ! generate valence expansion sets
454 2 : maxl = atom%state%maxl_occ
455 2 : CALL section_vals_val_get(grb_section, "NUM_GTO_EXTENDED", i_val=num_gto)
456 2 : CALL section_vals_val_get(grb_section, "EXTENSION_BASIS", i_vals=num_bas)
457 2 : next_bas(0:lmat) = 0
458 2 : IF (num_bas(1) == -1) THEN
459 0 : DO l = 0, maxl
460 0 : next_bas(l) = maxl - l + 1
461 : END DO
462 : ELSE
463 2 : n = MIN(SIZE(num_bas, 1), 4)
464 6 : next_bas(0:n - 1) = num_bas(1:n)
465 : END IF
466 2 : next_prim = 0
467 14 : DO l = 0, lmat
468 14 : IF (next_bas(l) > 0) next_prim(l) = num_gto
469 : END DO
470 2 : IF (iw > 0) THEN
471 2 : CALL basis_label(abas, next_prim, next_bas)
472 2 : WRITE (iw, '(/,A,A)') " Extension basis set ", TRIM(abas)
473 : END IF
474 14 : n = MAXVAL(next_prim)
475 14 : m = MAXVAL(next_bas)
476 11 : ALLOCATE (ebasis(n, n, 0:lmat), ale(n))
477 2 : basis_vrb => vbasis(0)%basis
478 2 : amin = atom%basis%aval/atom%basis%cval**1.5_dp
479 6 : DO i = 1, n
480 6 : ale(i) = amin*atom%basis%cval**(i - 1)
481 : END DO
482 134 : ebasis = 0._dp
483 3 : ALLOCATE (rho(n))
484 2 : rconf = 2.0_dp*atom%potential%scon
485 14 : DO l = 0, lmat
486 12 : IF (next_bas(l) < 1) CYCLE
487 2 : zval = SQRT(REAL(2*l + 2, dp))*REAL(l + 1, dp)/(2._dp*rmax)
488 2 : CALL hydrogenic(zval, rconf, l, ale, n, rho, ebasis(:, :, l))
489 2 : IF (iw > 0) WRITE (iw, '(T5,A,i5,T66,A,F10.4)') &
490 4 : " Extension basis set contraction for lval=", l, "zval=", zval
491 : END DO
492 2 : DEALLOCATE (rho)
493 : ! check for condition numbers
494 2 : IF (iw > 0) WRITE (iw, '(/,A)') " Condition Number of Extended Basis Sets"
495 2 : CALL init_orbital_pointers(lmat)
496 2 : CALL init_spherical_harmonics(lmat, 0)
497 10 : DO ider = 0, nder
498 8 : NULLIFY (basis_vrb)
499 152 : ALLOCATE (basis_vrb)
500 : NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
501 : basis_vrb%dbf, basis_vrb%ddbf)
502 : ! allocate and initialize the atomic grid
503 : NULLIFY (basis_vrb%grid)
504 8 : CALL allocate_grid_atom(basis_vrb%grid)
505 8 : CALL create_grid_atom(basis_vrb%grid, ngp, 1, 1, 0, quadtype)
506 8 : basis_vrb%grid%nr = ngp
507 : !
508 8 : basis_vrb%eps_eig = basis_ref%eps_eig
509 8 : basis_vrb%geometrical = .FALSE.
510 8 : basis_vrb%basis_type = CGTO_BASIS
511 104 : basis_vrb%nprim = basis%nprim + next_prim
512 56 : basis_vrb%nbas = 0
513 24 : DO l = 0, state%maxl_occ
514 24 : basis_vrb%nbas(l) = state%maxn_occ(l) + ider + next_bas(l)
515 : END DO
516 56 : m = MAXVAL(basis_vrb%nprim)
517 24 : ALLOCATE (basis_vrb%am(m, 0:lmat))
518 : ! exponents
519 8 : m = SIZE(basis%am, 1)
520 680 : basis_vrb%am(1:m, :) = basis%am(1:m, :)
521 : n = SIZE(ale, 1)
522 24 : DO l = 0, state%maxl_occ
523 56 : basis_vrb%am(m + 1:m + n, l) = ale(1:n)
524 : END DO
525 : ! contractions
526 56 : m = MAXVAL(basis_vrb%nprim)
527 56 : n = MAXVAL(basis_vrb%nbas)
528 40 : ALLOCATE (basis_vrb%cm(m, n, 0:lmat))
529 1664 : basis_vrb%cm = 0.0_dp
530 24 : DO l = 0, state%maxl_occ
531 16 : m = basis%nprim(l)
532 16 : n = state%maxn_occ(l) + ider
533 296 : basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
534 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)
535 : END DO
536 :
537 : ! initialize basis function on a radial grid
538 8 : nr = basis_vrb%grid%nr
539 56 : m = MAXVAL(basis_vrb%nbas)
540 40 : ALLOCATE (basis_vrb%bf(nr, m, 0:lmat))
541 24 : ALLOCATE (basis_vrb%dbf(nr, m, 0:lmat))
542 24 : ALLOCATE (basis_vrb%ddbf(nr, m, 0:lmat))
543 67424 : basis_vrb%bf = 0._dp
544 67424 : basis_vrb%dbf = 0._dp
545 67424 : basis_vrb%ddbf = 0._dp
546 56 : DO l = 0, lmat
547 184 : DO i = 1, basis_vrb%nprim(l)
548 128 : al = basis_vrb%am(i, l)
549 51376 : DO k = 1, nr
550 51200 : rk = basis_vrb%grid%rad(k)
551 51200 : ear = EXP(-al*basis_vrb%grid%rad(k)**2)
552 227328 : DO j = 1, basis_vrb%nbas(l)
553 176000 : basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
554 : basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
555 176000 : + (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
556 : basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
557 : (REAL(l*(l - 1), dp)*rk**(l - 2) - 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + &
558 227200 : 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
559 : END DO
560 : END DO
561 : END DO
562 : END DO
563 :
564 8 : IF (iw > 0) THEN
565 8 : CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
566 8 : WRITE (iw, '(A,A)') " Basis set ", TRIM(abas)
567 : END IF
568 8 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
569 8 : cradx = crad*1.00_dp
570 8 : CALL atom_basis_condnum(basis_vrb, cradx, cnum)
571 8 : IF (iw > 0) WRITE (iw, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
572 8 : cradx = crad*1.10_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.20_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 34 : vbasis(nder + 1 + ider)%basis => basis_vrb
579 : END DO
580 2 : CALL deallocate_orbital_pointers
581 2 : CALL deallocate_spherical_harmonics
582 :
583 : ! Tests for energy
584 2 : energy_ref = atom_ref%energy%etot
585 2 : IF (iw > 0) WRITE (iw, '(/,A,A)') " Basis set tests "
586 2 : IF (iw > 0) WRITE (iw, '(T10,A,T59,F22.9)') " Reference Energy [a.u.] ", energy_ref
587 18 : DO ider = 0, 2*nder + 1
588 : ! generate an atom type
589 16 : NULLIFY (atom_test)
590 16 : CALL create_atom_type(atom_test)
591 : CALL copy_atom_basics(atom_ref, atom_test, state=.TRUE., potential=.TRUE., &
592 16 : optimization=.TRUE., xc=.TRUE.)
593 16 : basis_grb => vbasis(ider)%basis
594 16 : NULLIFY (orbitals)
595 112 : mo = MAXVAL(atom_test%state%maxn_calc)
596 112 : mb = MAXVAL(basis_grb%nbas)
597 16 : CALL create_atom_orbs(orbitals, mb, mo)
598 16 : CALL set_atom(atom_test, orbitals=orbitals, basis=basis_grb)
599 : ! calculate integrals
600 3408 : ALLOCATE (atint)
601 16 : CALL atom_int_setup(atint, basis_grb, potential=atom_test%potential, eri_coulomb=.FALSE., eri_exchange=.FALSE.)
602 16 : CALL atom_ppint_setup(atint, basis_grb, potential=atom_test%potential)
603 16 : IF (atom_test%pp_calc) THEN
604 16 : NULLIFY (atint%tzora, atint%hdkh)
605 : ELSE
606 : ! relativistic correction terms
607 0 : CALL atom_relint_setup(atint, basis_grb, atom_test%relativistic, zcore=REAL(atom_test%z, dp))
608 : END IF
609 16 : CALL set_atom(atom_test, integrals=atint)
610 : !
611 16 : CALL calculate_atom(atom_test, iw=0)
612 16 : IF (ider <= nder) THEN
613 8 : energy_vb(ider) = atom_test%energy%etot
614 16 : IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (VB)", ider, " Energy [a.u.] ", &
615 16 : energy_ref - energy_vb(ider), energy_vb(ider)
616 : ELSE
617 8 : i = ider - nder - 1
618 8 : energy_ex(i) = atom_test%energy%etot
619 16 : IF (iw > 0) WRITE (iw, '(T10,A,i1,A,T40,F13.9,T59,F22.9)') " GRB (EX)", i, " Energy [a.u.] ", &
620 16 : energy_ref - energy_ex(i), energy_ex(i)
621 : END IF
622 16 : CALL atom_int_release(atint)
623 16 : CALL atom_ppint_release(atint)
624 16 : CALL atom_relint_release(atint)
625 16 : DEALLOCATE (atom_test%state, atom_test%potential, atint)
626 18 : CALL release_atom_type(atom_test)
627 : END DO
628 :
629 : ! Quickstep normalization polarization basis
630 18 : DO l = 0, 7
631 16 : expzet = 0.25_dp*REAL(2*l + 3, dp)
632 16 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
633 50 : DO i = 1, num_pol
634 32 : zeta = (2._dp*alp(i))**expzet
635 176 : pbasis(i, 1:num_pol, l) = pbasis(i, 1:num_pol, l)*prefac/zeta
636 : END DO
637 : END DO
638 : ! Quickstep normalization extended basis
639 14 : DO l = 0, lmat
640 12 : expzet = 0.25_dp*REAL(2*l + 3, dp)
641 12 : prefac = SQRT(rootpi/2._dp**(l + 2)*dfac(2*l + 1))
642 22 : DO i = 1, next_prim(l)
643 8 : zeta = (2._dp*ale(i))**expzet
644 32 : ebasis(i, 1:next_bas(l), l) = ebasis(i, 1:next_bas(l), l)*prefac/zeta
645 : END DO
646 : END DO
647 :
648 : ! Print basis sets
649 2 : CALL section_vals_val_get(grb_section, "NAME_BODY", c_val=basname)
650 2 : CALL open_file(file_name="GRB_BASIS", file_status="UNKNOWN", file_action="WRITE", unit_number=iunit)
651 : ! header info
652 8 : headline = ""
653 2 : headline(1) = "#"
654 2 : headline(2) = "# Generated with CP2K Atom Code"
655 2 : headline(3) = "#"
656 2 : CALL grb_print_basis(header=headline, iunit=iunit)
657 : ! valence basis
658 2 : basline(1) = ""
659 2 : WRITE (basline(1), "(T2,A)") ADJUSTL(ptable(atom_ref%z)%symbol)
660 10 : DO ider = 0, nder
661 8 : basline(1) = ""
662 8 : WRITE (basline(1), "(T2,A,T5,A,I1)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-VAL-", ider
663 : CALL grb_print_basis(header=basline, nprim=vbasis(ider)%basis%nprim(0), nbas=vbasis(ider)%basis%nbas, &
664 10 : al=vbasis(ider)%basis%am(:, 0), gcc=qbasis, iunit=iunit)
665 : END DO
666 : ! polarization basis
667 2 : maxl = atom_ref%state%maxl_occ
668 6 : DO l = maxl + 1, MIN(maxl + num_pol, 7)
669 4 : nbas = 0
670 14 : DO i = maxl + 1, l
671 14 : nbas(i) = l - i + 1
672 : END DO
673 4 : i = l - maxl
674 4 : basline(1) = ""
675 4 : WRITE (basline(1), "(T2,A,T5,A,I1)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-POL-", i
676 6 : CALL grb_print_basis(header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit)
677 : END DO
678 : ! extension set
679 14 : IF (SUM(next_bas) > 0) THEN
680 1 : basline(1) = ""
681 1 : WRITE (basline(1), "(T2,A,T5,A)") ADJUSTL(ptable(atom_ref%z)%symbol), TRIM(ADJUSTL(basname))//"-EXT"
682 1 : CALL grb_print_basis(header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit)
683 : END IF
684 : !
685 2 : CALL close_file(unit_number=iunit)
686 :
687 : ! clean up
688 2 : IF (ALLOCATED(pbasis)) DEALLOCATE (pbasis)
689 2 : IF (ALLOCATED(alp)) DEALLOCATE (alp)
690 2 : IF (ALLOCATED(ebasis)) DEALLOCATE (ebasis)
691 2 : DEALLOCATE (wfn, rbasis, qbasis, ale)
692 :
693 24 : DO ider = 0, 10
694 24 : IF (ASSOCIATED(vbasis(ider)%basis)) THEN
695 16 : CALL release_atom_basis(vbasis(ider)%basis)
696 16 : DEALLOCATE (vbasis(ider)%basis)
697 : END IF
698 : END DO
699 :
700 2 : CALL atom_int_release(atom%integrals)
701 2 : CALL atom_ppint_release(atom%integrals)
702 2 : CALL atom_relint_release(atom%integrals)
703 2 : CALL release_atom_basis(basis)
704 2 : DEALLOCATE (atom%potential, atom%state, atom%integrals, basis)
705 2 : CALL release_atom_type(atom)
706 :
707 2 : IF (iw > 0) WRITE (iw, '(" ",79("*"))')
708 :
709 18 : END SUBROUTINE atom_grb_construction
710 :
711 : ! **************************************************************************************************
712 : !> \brief Print geometrical response basis set.
713 : !> \param header banner to print on top of the basis set
714 : !> \param nprim number of primitive exponents
715 : !> \param nbas number of basis functions for the given angular momentum
716 : !> \param al list of the primitive exponents
717 : !> \param gcc array of contraction coefficients of size
718 : !> (index-of-the-primitive-exponent, index-of-the-contraction-set, angular-momentum)
719 : !> \param iunit output file unit
720 : !> \par History
721 : !> * 11.2016 created [Juerg Hutter]
722 : ! **************************************************************************************************
723 15 : SUBROUTINE grb_print_basis(header, nprim, nbas, al, gcc, iunit)
724 : CHARACTER(len=*), DIMENSION(:), INTENT(IN), &
725 : OPTIONAL :: header
726 : INTEGER, INTENT(IN), OPTIONAL :: nprim
727 : INTEGER, DIMENSION(0:), INTENT(IN), OPTIONAL :: nbas
728 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: al
729 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN), &
730 : OPTIONAL :: gcc
731 : INTEGER, INTENT(IN) :: iunit
732 :
733 : INTEGER :: i, j, l, lmax, lmin, nval
734 :
735 15 : IF (PRESENT(header)) THEN
736 34 : DO i = 1, SIZE(header, 1)
737 34 : IF (header(i) .NE. "") THEN
738 19 : WRITE (iunit, "(A)") TRIM(header(i))
739 : END IF
740 : END DO
741 : END IF
742 :
743 15 : IF (PRESENT(nprim)) THEN
744 13 : IF (nprim > 0) THEN
745 13 : CPASSERT(PRESENT(nbas))
746 13 : CPASSERT(PRESENT(al))
747 13 : CPASSERT(PRESENT(gcc))
748 :
749 34 : DO i = LBOUND(nbas, 1), UBOUND(nbas, 1)
750 21 : IF (nbas(i) > 0) THEN
751 13 : lmin = i
752 13 : EXIT
753 : END IF
754 : END DO
755 63 : DO i = UBOUND(nbas, 1), LBOUND(nbas, 1), -1
756 63 : IF (nbas(i) > 0) THEN
757 13 : lmax = i
758 13 : EXIT
759 : END IF
760 : END DO
761 :
762 13 : nval = lmax
763 13 : WRITE (iunit, *) " 1"
764 13 : WRITE (iunit, "(40I3)") nval, lmin, lmax, nprim, (nbas(l), l=lmin, lmax)
765 81 : DO i = nprim, 1, -1
766 68 : WRITE (iunit, "(G20.12)", advance="no") al(i)
767 212 : DO l = lmin, lmax
768 544 : DO j = 1, nbas(l)
769 476 : WRITE (iunit, "(F16.10)", advance="no") gcc(i, j, l)
770 : END DO
771 : END DO
772 81 : WRITE (iunit, *)
773 : END DO
774 13 : WRITE (iunit, *)
775 : END IF
776 : END IF
777 :
778 15 : END SUBROUTINE grb_print_basis
779 :
780 : ! **************************************************************************************************
781 : !> \brief Compose the basis set label:
782 : !> (np(0)'s'np(1)'p'...) -> [nb(0)'s'nb(1)'p'...] .
783 : !> \param label basis set label
784 : !> \param np number of primitive basis functions per angular momentum
785 : !> \param nb number of contracted basis functions per angular momentum
786 : !> \par History
787 : !> * 11.2016 created [Juerg Hutter]
788 : ! **************************************************************************************************
789 18 : SUBROUTINE basis_label(label, np, nb)
790 : CHARACTER(len=*), INTENT(out) :: label
791 : INTEGER, DIMENSION(0:), INTENT(in) :: np, nb
792 :
793 : INTEGER :: i, l, lmax
794 : CHARACTER(len=1), DIMENSION(0:7), PARAMETER :: &
795 : 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 164 : 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 35588 : ALLOCATE (atint)
856 : ! calculate integrals
857 164 : NULLIFY (pot)
858 164 : eri_c = .FALSE.
859 164 : eri_e = .FALSE.
860 164 : pot => atom%potential
861 164 : zval = atom%z
862 164 : reltyp = atom%relativistic
863 164 : do_eric = atom%coulomb_integral_type
864 164 : do_erie = atom%exchange_integral_type
865 164 : IF (do_eric == do_analytic) eri_c = .TRUE.
866 164 : IF (do_erie == do_analytic) eri_e = .TRUE.
867 : ! general integrals
868 164 : CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
869 : ! potential
870 164 : CALL atom_ppint_setup(atint, basis, potential=pot)
871 164 : IF (atom%pp_calc) THEN
872 164 : 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 164 : CALL set_atom(atom, basis=basis)
878 164 : CALL set_atom(atom, integrals=atint)
879 164 : CALL calculate_atom(atom, iw)
880 164 : afun = atom%energy%etot
881 164 : CALL atom_int_release(atint)
882 164 : CALL atom_ppint_release(atint)
883 164 : CALL atom_relint_release(atint)
884 164 : DEALLOCATE (atint)
885 164 : 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 94914 : 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 170 : IF (ostate%state == 2) THEN
988 7052 : basis%am = 0._dp
989 1148 : DO l = 0, lmat
990 3116 : DO i = 1, basis%nbas(l)
991 1968 : ll = i - 1 + basis%start(l)
992 2952 : basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
993 : END DO
994 : END DO
995 164 : basis%aval = x(1)*x(1)
996 164 : basis%cval = x(2)*x(2)
997 2368652 : basis%bf = 0._dp
998 2368652 : basis%dbf = 0._dp
999 2368652 : basis%ddbf = 0._dp
1000 164 : nr = basis%grid%nr
1001 1148 : DO l = 0, lmat
1002 3116 : DO i = 1, basis%nbas(l)
1003 1968 : al = basis%am(i, l)
1004 790152 : DO k = 1, nr
1005 787200 : rk = basis%grid%rad(k)
1006 787200 : ear = EXP(-al*basis%grid%rad(k)**2)
1007 787200 : basis%bf(k, i, l) = rk**l*ear
1008 787200 : 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 789168 : 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 164 : CALL grb_fit(atom, basis, ostate%f, 0)
1015 164 : fopt = MIN(fopt, ostate%f)
1016 : END IF
1017 :
1018 170 : IF (ostate%state == -1) EXIT
1019 :
1020 168 : CALL powell_optimize(ostate%nvar, x, ostate)
1021 :
1022 168 : 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 170 : 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 76 : IF (ostate%state == 2) THEN
1149 73 : aval = x(1)*x(1)
1150 73 : cval = x(2)*x(2)
1151 365 : DO i = 1, nbas
1152 365 : am(i) = aval*cval**(i - 1)
1153 : END DO
1154 73 : CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1155 73 : ostate%f = ener(1)
1156 73 : fopt = MIN(fopt, ostate%f)
1157 : END IF
1158 :
1159 76 : IF (ostate%state == -1) EXIT
1160 :
1161 75 : CALL powell_optimize(ostate%nvar, x, ostate)
1162 :
1163 75 : 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 76 : 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 79 : 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 79 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
1219 79 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: confmat, hmat, potmat, smat, tmat
1220 :
1221 79 : n = nbas
1222 948 : ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
1223 : ! calclulate overlap matrix
1224 79 : CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
1225 : ! calclulate kinetic energy matrix
1226 79 : CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
1227 : ! calclulate core potential matrix
1228 79 : CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
1229 : ! calclulate confinement potential matrix
1230 79 : cf = 0.1_dp
1231 79 : k = 10
1232 79 : CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
1233 : ! Hamiltionian
1234 1659 : 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 79 : lwork = 100*n
1237 395 : ALLOCATE (w(n), work(lwork))
1238 79 : CALL dsygv(1, "V", "U", n, hmat, n, smat, n, w, work, lwork, info)
1239 79 : CPASSERT(info == 0)
1240 1659 : orb(1:n, 1:n) = hmat(1:n, 1:n)
1241 395 : ener(1:n) = w(1:n)
1242 79 : DEALLOCATE (w, work)
1243 79 : DEALLOCATE (smat, tmat, potmat, confmat, hmat)
1244 :
1245 79 : END SUBROUTINE hydrogenic
1246 :
1247 44 : END MODULE atom_grb
|
