Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 lapack, ONLY: lapack_ssygv
41 : USE mathconstants, ONLY: dfac,&
42 : rootpi
43 : USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
44 : init_orbital_pointers
45 : USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
46 : init_spherical_harmonics
47 : USE periodic_table, ONLY: ptable
48 : USE physcon, ONLY: bohr
49 : USE powell, ONLY: opt_state_type,&
50 : powell_optimize
51 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
52 : create_grid_atom
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : TYPE basis_p_type
58 : TYPE(atom_basis_type), POINTER :: basis => NULL()
59 : END TYPE basis_p_type
60 :
61 : PRIVATE
62 : PUBLIC :: atom_grb_construction
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_grb'
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief Construct geometrical response basis set.
70 : !> \param atom_info information about the atomic kind. Two-dimensional array of size
71 : !> (electronic-configuration, electronic-structure-method)
72 : !> \param atom_section ATOM input section
73 : !> \param iw output file unit
74 : !> \par History
75 : !> * 11.2016 created [Juerg Hutter]
76 : ! **************************************************************************************************
77 2 : SUBROUTINE atom_grb_construction(atom_info, atom_section, iw)
78 :
79 : TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
80 : TYPE(section_vals_type), POINTER :: atom_section
81 : INTEGER, INTENT(IN) :: iw
82 :
83 : CHARACTER(len=default_string_length) :: abas, basname
84 : CHARACTER(len=default_string_length), DIMENSION(1) :: basline
85 : CHARACTER(len=default_string_length), DIMENSION(3) :: headline
86 : INTEGER :: i, ider, is, iunit, j, k, l, lhomo, ll, &
87 : lval, m, maxl, mb, method, mo, n, &
88 : nder, ngp, nhomo, nr, num_gto, &
89 : num_pol, quadtype, s1, s2
90 : INTEGER, DIMENSION(0:7) :: nbas
91 : INTEGER, DIMENSION(0:lmat) :: next_bas, next_prim
92 2 : INTEGER, DIMENSION(:), POINTER :: num_bas
93 : REAL(KIND=dp) :: al, amin, aval, cnum, crad, cradx, cval, delta, dene, ear, emax, &
94 : energy_ex(0:lmat), energy_ref, energy_vb(0:lmat), expzet, fhomo, o, prefac, rconf, rk, &
95 : rmax, scon, zeta, zval
96 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ale, alp, rho
97 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat
98 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ebasis, pbasis, qbasis, rbasis
99 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: wfn
100 2 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: ovlp
101 : TYPE(atom_basis_type), POINTER :: basis, basis_grb, basis_ref, basis_vrb
102 : TYPE(atom_integrals), POINTER :: atint
103 : TYPE(atom_orbitals), POINTER :: orbitals
104 : TYPE(atom_state), POINTER :: state
105 : TYPE(atom_type), POINTER :: atom, atom_ref, atom_test
106 24 : TYPE(basis_p_type), DIMENSION(0:10) :: vbasis
107 : TYPE(section_vals_type), POINTER :: grb_section, powell_section
108 :
109 2 : IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T28,A,/," ",79("*"))') "GEOMETRICAL RESPONSE BASIS"
110 :
111 24 : DO i = 0, 10
112 24 : NULLIFY (vbasis(i)%basis)
113 : END DO
114 : ! make some basic checks
115 6 : is = SIZE(atom_info)
116 2 : IF (iw > 0 .AND. is > 1) THEN
117 0 : WRITE (iw, '(/,A,/)') " WARNING: Only use first electronic structure/method for basis set generation"
118 : END IF
119 2 : atom_ref => atom_info(1, 1)%atom
120 :
121 : ! check method
122 2 : method = atom_ref%method_type
123 0 : SELECT CASE (method)
124 : CASE (do_rks_atom, do_rhf_atom)
125 : ! restricted methods are okay
126 : CASE (do_uks_atom, do_uhf_atom, do_rohf_atom)
127 0 : CPABORT("Unrestricted methods not allowed for GRB generation")
128 : CASE DEFAULT
129 2 : CPABORT("")
130 : END SELECT
131 :
132 : ! input for basis optimization
133 2 : grb_section => section_vals_get_subs_vals(atom_section, "PRINT%GEOMETRICAL_RESPONSE_BASIS")
134 :
135 : ! generate an atom type
136 2 : NULLIFY (atom)
137 2 : CALL create_atom_type(atom)
138 2 : CALL copy_atom_basics(atom_ref, atom, state=.TRUE., potential=.TRUE., optimization=.TRUE., xc=.TRUE.)
139 : ! set confinement potential
140 2 : atom%potential%confinement = .TRUE.
141 2 : atom%potential%conf_type = barrier_conf
142 2 : atom%potential%acon = 200._dp
143 2 : atom%potential%rcon = 4._dp
144 2 : CALL section_vals_val_get(grb_section, "CONFINEMENT", r_val=scon)
145 2 : atom%potential%scon = scon
146 : ! generate main block geometrical exponents
147 2 : basis_ref => atom_ref%basis
148 38 : ALLOCATE (basis)
149 : NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
150 : ! get information on quadrature type and number of grid points
151 : ! allocate and initialize the atomic grid
152 : NULLIFY (basis%grid)
153 2 : CALL allocate_grid_atom(basis%grid)
154 2 : CALL section_vals_val_get(grb_section, "QUADRATURE", i_val=quadtype)
155 2 : CALL section_vals_val_get(grb_section, "GRID_POINTS", i_val=ngp)
156 2 : IF (ngp <= 0) &
157 0 : CPABORT("# point radial grid < 0")
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 350 : rbasis = 0._dp
247 350 : 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("")
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 2304 : 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 1536 : 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 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)))
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 802 : 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 169 : 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 2 : 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 134 : ebasis = 0._dp
484 3 : 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 : 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) .NE. "") 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 :: &
796 : lq = (/"s", "p", "d", "f", "g", "h", "i", "k"/)
797 :
798 18 : label = ""
799 18 : lmax = MIN(UBOUND(np, 1), UBOUND(nb, 1), 7)
800 18 : i = 1
801 18 : label(i:i) = "("
802 126 : DO l = 0, lmax
803 126 : IF (np(l) > 0) THEN
804 34 : i = i + 1
805 34 : IF (np(l) > 9) THEN
806 8 : WRITE (label(i:i + 1), "(I2)") np(l)
807 8 : i = i + 2
808 : ELSE
809 26 : WRITE (label(i:i), "(I1)") np(l)
810 26 : i = i + 1
811 : END IF
812 34 : label(i:i) = lq(l)
813 : END IF
814 : END DO
815 18 : i = i + 1
816 18 : label(i:i + 6) = ") --> ["
817 18 : i = i + 6
818 126 : DO l = 0, lmax
819 126 : IF (nb(l) > 0) THEN
820 34 : i = i + 1
821 34 : IF (nb(l) > 9) THEN
822 0 : WRITE (label(i:i + 1), "(I2)") nb(l)
823 0 : i = i + 2
824 : ELSE
825 34 : WRITE (label(i:i), "(I1)") nb(l)
826 34 : i = i + 1
827 : END IF
828 34 : label(i:i) = lq(l)
829 : END IF
830 : END DO
831 18 : i = i + 1
832 18 : label(i:i) = "]"
833 :
834 18 : END SUBROUTINE basis_label
835 :
836 : ! **************************************************************************************************
837 : !> \brief Compute the total energy for the given atomic kind and basis set.
838 : !> \param atom information about the atomic kind
839 : !> \param basis basis set to fit
840 : !> \param afun (output) atomic total energy
841 : !> \param iw output file unit
842 : !> \par History
843 : !> * 11.2016 created [Juerg Hutter]
844 : ! **************************************************************************************************
845 164 : SUBROUTINE grb_fit(atom, basis, afun, iw)
846 : TYPE(atom_type), POINTER :: atom
847 : TYPE(atom_basis_type), POINTER :: basis
848 : REAL(dp), INTENT(OUT) :: afun
849 : INTEGER, INTENT(IN) :: iw
850 :
851 : INTEGER :: do_eric, do_erie, reltyp, zval
852 : LOGICAL :: eri_c, eri_e
853 : TYPE(atom_integrals), POINTER :: atint
854 : TYPE(atom_potential_type), POINTER :: pot
855 :
856 35588 : ALLOCATE (atint)
857 : ! calculate integrals
858 164 : NULLIFY (pot)
859 164 : eri_c = .FALSE.
860 164 : eri_e = .FALSE.
861 164 : pot => atom%potential
862 164 : zval = atom%z
863 164 : reltyp = atom%relativistic
864 164 : do_eric = atom%coulomb_integral_type
865 164 : do_erie = atom%exchange_integral_type
866 164 : IF (do_eric == do_analytic) eri_c = .TRUE.
867 164 : IF (do_erie == do_analytic) eri_e = .TRUE.
868 : ! general integrals
869 164 : CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
870 : ! potential
871 164 : CALL atom_ppint_setup(atint, basis, potential=pot)
872 164 : IF (atom%pp_calc) THEN
873 164 : NULLIFY (atint%tzora, atint%hdkh)
874 : ELSE
875 : ! relativistic correction terms
876 0 : CALL atom_relint_setup(atint, basis, reltyp, zcore=REAL(zval, dp))
877 : END IF
878 164 : CALL set_atom(atom, basis=basis)
879 164 : CALL set_atom(atom, integrals=atint)
880 164 : CALL calculate_atom(atom, iw)
881 164 : afun = atom%energy%etot
882 164 : CALL atom_int_release(atint)
883 164 : CALL atom_ppint_release(atint)
884 164 : CALL atom_relint_release(atint)
885 164 : DEALLOCATE (atint)
886 164 : END SUBROUTINE grb_fit
887 :
888 : ! **************************************************************************************************
889 : !> \brief Copy basic information about the atomic kind.
890 : !> \param atom_ref atom to copy
891 : !> \param atom new atom to create
892 : !> \param state also copy electronic state and occupation numbers
893 : !> \param potential also copy pseudo-potential
894 : !> \param optimization also copy optimization procedure
895 : !> \param xc also copy the XC input section
896 : !> \par History
897 : !> * 11.2016 created [Juerg Hutter]
898 : ! **************************************************************************************************
899 18 : SUBROUTINE copy_atom_basics(atom_ref, atom, state, potential, optimization, xc)
900 : TYPE(atom_type), POINTER :: atom_ref, atom
901 : LOGICAL, INTENT(IN), OPTIONAL :: state, potential, optimization, xc
902 :
903 18 : atom%z = atom_ref%z
904 18 : atom%zcore = atom_ref%zcore
905 18 : atom%pp_calc = atom_ref%pp_calc
906 18 : atom%method_type = atom_ref%method_type
907 18 : atom%relativistic = atom_ref%relativistic
908 18 : atom%coulomb_integral_type = atom_ref%coulomb_integral_type
909 18 : atom%exchange_integral_type = atom_ref%exchange_integral_type
910 :
911 18 : NULLIFY (atom%potential, atom%state, atom%xc_section)
912 18 : NULLIFY (atom%basis, atom%integrals, atom%orbitals, atom%fmat)
913 :
914 18 : IF (PRESENT(state)) THEN
915 18 : IF (state) THEN
916 6660 : ALLOCATE (atom%state)
917 18 : atom%state = atom_ref%state
918 : END IF
919 : END IF
920 :
921 18 : IF (PRESENT(potential)) THEN
922 18 : IF (potential) THEN
923 94914 : ALLOCATE (atom%potential)
924 18 : atom%potential = atom_ref%potential
925 : END IF
926 : END IF
927 :
928 18 : IF (PRESENT(optimization)) THEN
929 18 : IF (optimization) THEN
930 18 : atom%optimization = atom_ref%optimization
931 : END IF
932 : END IF
933 :
934 18 : IF (PRESENT(xc)) THEN
935 18 : IF (xc) THEN
936 18 : atom%xc_section => atom_ref%xc_section
937 : END IF
938 : END IF
939 :
940 18 : END SUBROUTINE copy_atom_basics
941 :
942 : ! **************************************************************************************************
943 : !> \brief Optimise a geometrical response basis set.
944 : !> \param atom information about the atomic kind
945 : !> \param basis basis set to fit
946 : !> \param iunit output file unit
947 : !> \param powell_section POWELL input section
948 : !> \par History
949 : !> * 11.2016 created [Juerg Hutter]
950 : ! **************************************************************************************************
951 2 : SUBROUTINE atom_fit_grb(atom, basis, iunit, powell_section)
952 : TYPE(atom_type), POINTER :: atom
953 : TYPE(atom_basis_type), POINTER :: basis
954 : INTEGER, INTENT(IN) :: iunit
955 : TYPE(section_vals_type), POINTER :: powell_section
956 :
957 : INTEGER :: i, k, l, ll, n10, nr
958 : REAL(KIND=dp) :: al, cnum, crad, cradx, ear, fopt, rk
959 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x
960 : TYPE(opt_state_type) :: ostate
961 :
962 0 : CPASSERT(basis%geometrical)
963 :
964 2 : CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
965 2 : CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
966 2 : CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
967 :
968 2 : ostate%nvar = 2
969 2 : ALLOCATE (x(2))
970 2 : x(1) = SQRT(basis%aval)
971 2 : x(2) = SQRT(basis%cval)
972 :
973 2 : ostate%nf = 0
974 2 : ostate%iprint = 1
975 2 : ostate%unit = iunit
976 :
977 2 : ostate%state = 0
978 2 : IF (iunit > 0) THEN
979 2 : WRITE (iunit, '(/," POWELL| Start optimization procedure")')
980 2 : WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
981 : END IF
982 2 : n10 = MAX(ostate%maxfun/100, 1)
983 :
984 2 : fopt = HUGE(0._dp)
985 :
986 : DO
987 :
988 170 : IF (ostate%state == 2) THEN
989 7052 : basis%am = 0._dp
990 1148 : DO l = 0, lmat
991 3116 : DO i = 1, basis%nbas(l)
992 1968 : ll = i - 1 + basis%start(l)
993 2952 : basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
994 : END DO
995 : END DO
996 164 : basis%aval = x(1)*x(1)
997 164 : basis%cval = x(2)*x(2)
998 2368652 : basis%bf = 0._dp
999 2368652 : basis%dbf = 0._dp
1000 2368652 : basis%ddbf = 0._dp
1001 164 : nr = basis%grid%nr
1002 1148 : DO l = 0, lmat
1003 3116 : DO i = 1, basis%nbas(l)
1004 1968 : al = basis%am(i, l)
1005 790152 : DO k = 1, nr
1006 787200 : rk = basis%grid%rad(k)
1007 787200 : ear = EXP(-al*basis%grid%rad(k)**2)
1008 787200 : basis%bf(k, i, l) = rk**l*ear
1009 787200 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1010 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
1011 789168 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1012 : END DO
1013 : END DO
1014 : END DO
1015 164 : CALL grb_fit(atom, basis, ostate%f, 0)
1016 164 : fopt = MIN(fopt, ostate%f)
1017 : END IF
1018 :
1019 170 : IF (ostate%state == -1) EXIT
1020 :
1021 168 : CALL powell_optimize(ostate%nvar, x, ostate)
1022 :
1023 168 : IF (ostate%nf == 2 .AND. iunit > 0) THEN
1024 2 : WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1025 : END IF
1026 170 : IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1027 : WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1028 2 : INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
1029 : END IF
1030 :
1031 : END DO
1032 :
1033 2 : ostate%state = 8
1034 2 : CALL powell_optimize(ostate%nvar, x, ostate)
1035 :
1036 2 : IF (iunit > 0) THEN
1037 2 : WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1038 2 : WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1039 : END IF
1040 : ! x->basis
1041 86 : basis%am = 0._dp
1042 14 : DO l = 0, lmat
1043 38 : DO i = 1, basis%nbas(l)
1044 24 : ll = i - 1 + basis%start(l)
1045 36 : basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
1046 : END DO
1047 : END DO
1048 2 : basis%aval = x(1)*x(1)
1049 2 : basis%cval = x(2)*x(2)
1050 28886 : basis%bf = 0._dp
1051 28886 : basis%dbf = 0._dp
1052 28886 : basis%ddbf = 0._dp
1053 2 : nr = basis%grid%nr
1054 14 : DO l = 0, lmat
1055 38 : DO i = 1, basis%nbas(l)
1056 24 : al = basis%am(i, l)
1057 9636 : DO k = 1, nr
1058 9600 : rk = basis%grid%rad(k)
1059 9600 : ear = EXP(-al*basis%grid%rad(k)**2)
1060 9600 : basis%bf(k, i, l) = rk**l*ear
1061 9600 : basis%dbf(k, i, l) = (REAL(l, dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1062 : basis%ddbf(k, i, l) = (REAL(l*(l - 1), dp)*rk**(l - 2) - &
1063 9624 : 2._dp*al*REAL(2*l + 1, dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1064 : END DO
1065 : END DO
1066 : END DO
1067 :
1068 2 : DEALLOCATE (x)
1069 :
1070 : ! final result
1071 2 : IF (iunit > 0) THEN
1072 2 : WRITE (iunit, '(/,A)') " Optimized Geometrical GTO basis set"
1073 2 : WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", basis%aval, &
1074 4 : " Proportionality factor: ", basis%cval
1075 14 : DO l = 0, lmat
1076 14 : WRITE (iunit, '(T41,A,I2,T76,I5)') " Number of exponents for l=", l, basis%nbas(l)
1077 : END DO
1078 : END IF
1079 :
1080 2 : IF (iunit > 0) WRITE (iunit, '(/,A)') " Condition number of uncontracted basis set"
1081 2 : crad = 2.0_dp*ptable(atom%z)%covalent_radius*bohr
1082 2 : CALL init_orbital_pointers(lmat)
1083 2 : CALL init_spherical_harmonics(lmat, 0)
1084 2 : cradx = crad*1.00_dp
1085 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1086 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1087 2 : cradx = crad*1.10_dp
1088 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1089 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1090 2 : cradx = crad*1.20_dp
1091 2 : CALL atom_basis_condnum(basis, cradx, cnum)
1092 2 : IF (iunit > 0) WRITE (iunit, '(T5,A,F15.3,T50,A,F14.4)') " Lattice constant:", cradx, "Condition number:", cnum
1093 2 : CALL deallocate_orbital_pointers
1094 2 : CALL deallocate_spherical_harmonics
1095 :
1096 8 : END SUBROUTINE atom_fit_grb
1097 :
1098 : ! **************************************************************************************************
1099 : !> \brief Optimize 'aval' and 'cval' parameters which define the geometrical response basis set.
1100 : !> \param zval nuclear charge
1101 : !> \param rconf confinement radius
1102 : !> \param lval angular momentum
1103 : !> \param aval (input/output) exponent of the first Gaussian basis function in the series
1104 : !> \param cval (input/output) factor of geometrical series
1105 : !> \param nbas number of basis functions
1106 : !> \param iunit output file unit
1107 : !> \param powell_section POWELL input section
1108 : !> \par History
1109 : !> * 11.2016 created [Juerg Hutter]
1110 : ! **************************************************************************************************
1111 1 : SUBROUTINE atom_fit_pol(zval, rconf, lval, aval, cval, nbas, iunit, powell_section)
1112 : REAL(KIND=dp), INTENT(IN) :: zval, rconf
1113 : INTEGER, INTENT(IN) :: lval
1114 : REAL(KIND=dp), INTENT(INOUT) :: aval, cval
1115 : INTEGER, INTENT(IN) :: nbas, iunit
1116 : TYPE(section_vals_type), POINTER :: powell_section
1117 :
1118 : INTEGER :: i, n10
1119 : REAL(KIND=dp) :: fopt, x(2)
1120 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: am, ener
1121 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1122 : TYPE(opt_state_type) :: ostate
1123 :
1124 7 : ALLOCATE (am(nbas), ener(nbas), orb(nbas, nbas))
1125 :
1126 1 : CALL section_vals_val_get(powell_section, "ACCURACY", r_val=ostate%rhoend)
1127 1 : CALL section_vals_val_get(powell_section, "STEP_SIZE", r_val=ostate%rhobeg)
1128 1 : CALL section_vals_val_get(powell_section, "MAX_FUN", i_val=ostate%maxfun)
1129 :
1130 1 : ostate%nvar = 2
1131 1 : x(1) = SQRT(aval)
1132 1 : x(2) = SQRT(cval)
1133 :
1134 1 : ostate%nf = 0
1135 1 : ostate%iprint = 1
1136 1 : ostate%unit = iunit
1137 :
1138 1 : ostate%state = 0
1139 1 : IF (iunit > 0) THEN
1140 1 : WRITE (iunit, '(/," POWELL| Start optimization procedure")')
1141 1 : WRITE (iunit, '(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
1142 : END IF
1143 1 : n10 = MAX(ostate%maxfun/100, 1)
1144 :
1145 1 : fopt = HUGE(0._dp)
1146 :
1147 : DO
1148 :
1149 76 : IF (ostate%state == 2) THEN
1150 73 : aval = x(1)*x(1)
1151 73 : cval = x(2)*x(2)
1152 365 : DO i = 1, nbas
1153 365 : am(i) = aval*cval**(i - 1)
1154 : END DO
1155 73 : CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1156 73 : ostate%f = ener(1)
1157 73 : fopt = MIN(fopt, ostate%f)
1158 : END IF
1159 :
1160 76 : IF (ostate%state == -1) EXIT
1161 :
1162 75 : CALL powell_optimize(ostate%nvar, x, ostate)
1163 :
1164 75 : IF (ostate%nf == 2 .AND. iunit > 0) THEN
1165 1 : WRITE (iunit, '(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1166 : END IF
1167 76 : IF (MOD(ostate%nf, n10) == 0 .AND. iunit > 0) THEN
1168 : WRITE (iunit, '(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1169 1 : INT(REAL(ostate%nf, dp)/REAL(ostate%maxfun, dp)*100._dp), fopt
1170 : END IF
1171 :
1172 : END DO
1173 :
1174 1 : ostate%state = 8
1175 1 : CALL powell_optimize(ostate%nvar, x, ostate)
1176 :
1177 1 : IF (iunit > 0) THEN
1178 1 : WRITE (iunit, '(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1179 1 : WRITE (iunit, '(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1180 : END IF
1181 : ! x->basis
1182 1 : aval = x(1)*x(1)
1183 1 : cval = x(2)*x(2)
1184 :
1185 : ! final result
1186 1 : IF (iunit > 0) THEN
1187 1 : WRITE (iunit, '(/,A,T51,A,T76,I5)') " Optimized Polarization basis set", &
1188 2 : " Number of exponents:", nbas
1189 1 : WRITE (iunit, '(A,F15.8,T41,A,F15.8)') " Initial exponent: ", aval, &
1190 2 : " Proportionality factor: ", cval
1191 : END IF
1192 :
1193 1 : DEALLOCATE (am, ener, orb)
1194 :
1195 1 : END SUBROUTINE atom_fit_pol
1196 :
1197 : ! **************************************************************************************************
1198 : !> \brief Calculate orbitals of a hydrogen-like atom.
1199 : !> \param zval nuclear charge
1200 : !> \param rconf confinement radius
1201 : !> \param lval angular momentum
1202 : !> \param am list of basis functions' exponents
1203 : !> \param nbas number of basis functions
1204 : !> \param ener orbital energies
1205 : !> \param orb expansion coefficients of atomic wavefunctions
1206 : !> \par History
1207 : !> * 11.2016 created [Juerg Hutter]
1208 : ! **************************************************************************************************
1209 79 : SUBROUTINE hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1210 : REAL(KIND=dp), INTENT(IN) :: zval, rconf
1211 : INTEGER, INTENT(IN) :: lval
1212 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: am
1213 : INTEGER, INTENT(IN) :: nbas
1214 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ener
1215 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: orb
1216 :
1217 : INTEGER :: info, k, lwork, n
1218 : REAL(KIND=dp) :: cf
1219 79 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
1220 79 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: confmat, hmat, potmat, smat, tmat
1221 :
1222 79 : n = nbas
1223 948 : ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
1224 : ! calclulate overlap matrix
1225 79 : CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
1226 : ! calclulate kinetic energy matrix
1227 79 : CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
1228 : ! calclulate core potential matrix
1229 79 : CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
1230 : ! calclulate confinement potential matrix
1231 79 : cf = 0.1_dp
1232 79 : k = 10
1233 79 : CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
1234 : ! Hamiltionian
1235 1659 : hmat(1:n, 1:n) = tmat(1:n, 1:n) - zval*potmat(1:n, 1:n) + cf*confmat(1:n, 1:n)
1236 : ! solve
1237 79 : lwork = 100*n
1238 395 : ALLOCATE (w(n), work(lwork))
1239 79 : CALL lapack_ssygv(1, "V", "U", n, hmat, n, smat, n, w, work, lwork, info)
1240 79 : CPASSERT(info == 0)
1241 1659 : orb(1:n, 1:n) = hmat(1:n, 1:n)
1242 395 : ener(1:n) = w(1:n)
1243 79 : DEALLOCATE (w, work)
1244 79 : DEALLOCATE (smat, tmat, potmat, confmat, hmat)
1245 :
1246 79 : END SUBROUTINE hydrogenic
1247 :
1248 44 : END MODULE atom_grb
|
