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 : ! **************************************************************************************************
9 : MODULE atom_pseudo
10 : USE atom_electronic_structure, ONLY: calculate_atom
11 : USE atom_fit, ONLY: atom_fit_pseudo
12 : USE atom_operators, ONLY: atom_int_release,&
13 : atom_int_setup,&
14 : atom_ppint_release,&
15 : atom_ppint_setup,&
16 : atom_relint_release,&
17 : atom_relint_setup
18 : USE atom_output, ONLY: atom_print_basis,&
19 : atom_print_info,&
20 : atom_print_method,&
21 : atom_print_orbitals,&
22 : atom_print_potential
23 : USE atom_types, ONLY: &
24 : atom_basis_type, atom_integrals, atom_optimization_type, atom_orbitals, atom_p_type, &
25 : atom_potential_type, atom_state, create_atom_orbs, create_atom_type, init_atom_basis, &
26 : init_atom_potential, lmat, read_atom_opt_section, release_atom_basis, &
27 : release_atom_potential, release_atom_type, set_atom
28 : USE atom_utils, ONLY: atom_consistent_method,&
29 : atom_set_occupation,&
30 : get_maxl_occ,&
31 : get_maxn_occ
32 : USE cp_log_handling, ONLY: cp_get_default_logger,&
33 : cp_logger_type
34 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
35 : cp_print_key_unit_nr
36 : USE input_constants, ONLY: do_analytic,&
37 : poly_conf
38 : USE input_section_types, ONLY: section_vals_get,&
39 : section_vals_get_subs_vals,&
40 : section_vals_type,&
41 : section_vals_val_get
42 : USE kinds, ONLY: default_string_length,&
43 : dp
44 : USE periodic_table, ONLY: nelem,&
45 : ptable
46 : USE physcon, ONLY: bohr
47 : #include "./base/base_uses.f90"
48 :
49 : IMPLICIT NONE
50 : PRIVATE
51 : PUBLIC :: atom_pseudo_opt
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_pseudo'
54 :
55 : ! **************************************************************************************************
56 :
57 : CONTAINS
58 :
59 : ! **************************************************************************************************
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param atom_section ...
64 : ! **************************************************************************************************
65 26 : SUBROUTINE atom_pseudo_opt(atom_section)
66 : TYPE(section_vals_type), POINTER :: atom_section
67 :
68 : CHARACTER(len=*), PARAMETER :: routineN = 'atom_pseudo_opt'
69 :
70 : CHARACTER(LEN=2) :: elem
71 : CHARACTER(LEN=default_string_length), &
72 26 : DIMENSION(:), POINTER :: tmpstringlist
73 : INTEGER :: ads, do_eric, do_erie, handle, i, im, &
74 : in, iw, k, l, maxl, mb, method, mo, &
75 : n_meth, n_rep, nr_gh, reltyp, zcore, &
76 : zval, zz
77 : INTEGER, DIMENSION(0:lmat) :: maxn
78 26 : INTEGER, DIMENSION(:), POINTER :: cn
79 : LOGICAL :: do_gh, eri_c, eri_e, graph, pp_calc
80 : REAL(KIND=dp) :: ne, nm
81 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: pocc
82 : TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
83 : TYPE(atom_integrals), POINTER :: ae_int, pp_int
84 : TYPE(atom_optimization_type) :: optimization
85 : TYPE(atom_orbitals), POINTER :: orbitals
86 26 : TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info, atom_refs
87 : TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
88 : TYPE(atom_state), POINTER :: state, statepp
89 : TYPE(cp_logger_type), POINTER :: logger
90 : TYPE(section_vals_type), POINTER :: basis_section, method_section, &
91 : opt_section, potential_section, &
92 : powell_section, xc_section
93 :
94 26 : CALL timeset(routineN, handle)
95 :
96 : ! What atom do we calculate
97 26 : CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
98 26 : CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
99 26 : zz = 0
100 140 : DO i = 1, nelem
101 140 : IF (ptable(i)%symbol == elem) THEN
102 : zz = i
103 : EXIT
104 : END IF
105 : END DO
106 26 : IF (zz /= 1) zval = zz
107 :
108 : ! read and set up information on the basis sets
109 962 : ALLOCATE (ae_basis, pp_basis)
110 26 : basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
111 26 : NULLIFY (ae_basis%grid)
112 26 : CALL init_atom_basis(ae_basis, basis_section, zval, "AA")
113 26 : NULLIFY (pp_basis%grid)
114 26 : basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
115 26 : CALL init_atom_basis(pp_basis, basis_section, zval, "AP")
116 :
117 : ! print general and basis set information
118 26 : logger => cp_get_default_logger()
119 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
120 26 : IF (iw > 0) CALL atom_print_info(zval, "Atomic Energy Calculation", iw)
121 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
122 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
123 26 : IF (iw > 0) THEN
124 8 : CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
125 8 : CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
126 : END IF
127 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
128 :
129 : ! read and setup information on the pseudopotential
130 26 : NULLIFY (potential_section)
131 26 : potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
132 273130 : ALLOCATE (ae_pot, p_pot)
133 26 : CALL init_atom_potential(p_pot, potential_section, zval)
134 26 : CALL init_atom_potential(ae_pot, potential_section, -1)
135 26 : IF (.NOT. p_pot%confinement .AND. .NOT. ae_pot%confinement) THEN
136 : !set default confinement potential
137 24 : p_pot%confinement = .TRUE.
138 24 : p_pot%conf_type = poly_conf
139 24 : p_pot%scon = 2.0_dp
140 24 : p_pot%acon = 0.5_dp
141 : ! this seems to be the default in the old code
142 24 : p_pot%rcon = (2._dp*ptable(zval)%covalent_radius*bohr)**2
143 24 : ae_pot%confinement = .TRUE.
144 24 : ae_pot%conf_type = poly_conf
145 24 : ae_pot%scon = 2.0_dp
146 24 : ae_pot%acon = 0.5_dp
147 : ! this seems to be the default in the old code
148 24 : ae_pot%rcon = (2._dp*ptable(zval)%covalent_radius*bohr)**2
149 : END IF
150 :
151 : ! if the ERI's are calculated analytically, we have to precalculate them
152 26 : eri_c = .FALSE.
153 26 : CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
154 26 : IF (do_eric == do_analytic) eri_c = .TRUE.
155 26 : eri_e = .FALSE.
156 26 : CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
157 26 : IF (do_erie == do_analytic) eri_e = .TRUE.
158 26 : CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
159 26 : CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
160 :
161 : ! information on the states to be calculated
162 26 : CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
163 26 : maxn = 0
164 26 : CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
165 52 : DO in = 1, MIN(SIZE(cn), 4)
166 52 : maxn(in - 1) = cn(in)
167 : END DO
168 182 : DO in = 0, lmat
169 182 : maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
170 : END DO
171 :
172 : ! read optimization section
173 26 : opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
174 26 : CALL read_atom_opt_section(optimization, opt_section)
175 :
176 : ! Check for the total number of electron configurations to be calculated
177 26 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
178 : ! Check for the total number of method types to be calculated
179 26 : method_section => section_vals_get_subs_vals(atom_section, "METHOD")
180 26 : CALL section_vals_get(method_section, n_repetition=n_meth)
181 :
182 : ! integrals
183 11050 : ALLOCATE (ae_int, pp_int)
184 :
185 260 : ALLOCATE (atom_info(n_rep, n_meth), atom_refs(n_rep, n_meth))
186 :
187 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
188 26 : IF (iw > 0) THEN
189 13 : WRITE (iw, '(/," ",79("*"))')
190 13 : WRITE (iw, '(" ",26("*"),A,25("*"))') " Calculate Reference States "
191 13 : WRITE (iw, '(" ",79("*"))')
192 : END IF
193 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
194 :
195 52 : DO in = 1, n_rep
196 78 : DO im = 1, n_meth
197 :
198 26 : NULLIFY (atom_info(in, im)%atom, atom_refs(in, im)%atom)
199 26 : CALL create_atom_type(atom_info(in, im)%atom)
200 26 : CALL create_atom_type(atom_refs(in, im)%atom)
201 :
202 26 : atom_info(in, im)%atom%optimization = optimization
203 26 : atom_refs(in, im)%atom%optimization = optimization
204 :
205 26 : atom_info(in, im)%atom%z = zval
206 26 : atom_refs(in, im)%atom%z = zval
207 26 : xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
208 26 : atom_info(in, im)%atom%xc_section => xc_section
209 26 : atom_refs(in, im)%atom%xc_section => xc_section
210 :
211 18850 : ALLOCATE (state, statepp)
212 :
213 : ! get the electronic configuration
214 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
215 26 : c_vals=tmpstringlist)
216 : ! all electron configurations have to be with full core
217 26 : pp_calc = INDEX(tmpstringlist(1), "CORE") /= 0
218 26 : CPASSERT(.NOT. pp_calc)
219 :
220 : ! set occupations
221 26 : CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
222 26 : state%maxl_occ = get_maxl_occ(state%occ)
223 182 : state%maxn_occ = get_maxn_occ(state%occ)
224 : ! set number of states to be calculated
225 26 : state%maxl_calc = MAX(maxl, state%maxl_occ)
226 26 : state%maxl_calc = MIN(lmat, state%maxl_calc)
227 182 : state%maxn_calc = 0
228 96 : DO k = 0, state%maxl_calc
229 70 : ads = 2
230 70 : IF (state%maxn_occ(k) == 0) ads = 1
231 70 : state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k) + ads)
232 96 : state%maxn_calc(k) = MIN(state%maxn_calc(k), ae_basis%nbas(k))
233 : END DO
234 1846 : state%core = 0._dp
235 26 : CALL set_atom(atom_refs(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
236 :
237 26 : IF (state%multiplicity /= -1) THEN
238 : ! set alpha and beta occupations
239 142 : state%occa = 0._dp
240 142 : state%occb = 0._dp
241 14 : DO l = 0, lmat
242 12 : nm = REAL((2*l + 1), KIND=dp)
243 20 : DO k = 1, 10
244 18 : ne = state%occupation(l, k)
245 18 : IF (ne == 0._dp) THEN !empty shell
246 : EXIT !assume there are no holes
247 6 : ELSEIF (ne == 2._dp*nm) THEN !closed shell
248 4 : state%occa(l, k) = nm
249 4 : state%occb(l, k) = nm
250 2 : ELSEIF (state%multiplicity == -2) THEN !High spin case
251 0 : state%occa(l, k) = MIN(ne, nm)
252 0 : state%occb(l, k) = MAX(0._dp, ne - nm)
253 : ELSE
254 2 : state%occa(l, k) = 0.5_dp*(ne + state%multiplicity - 1._dp)
255 2 : state%occb(l, k) = ne - state%occa(l, k)
256 : END IF
257 : END DO
258 : END DO
259 : END IF
260 :
261 : ! set occupations for pseudopotential calculation
262 26 : CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
263 26 : CALL atom_set_occupation(tmpstringlist, statepp%core, pocc)
264 1846 : zcore = zval - NINT(SUM(statepp%core))
265 26 : CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
266 :
267 3666 : statepp%occ = state%occ - statepp%core
268 1846 : statepp%occupation = 0._dp
269 182 : DO l = 0, lmat
270 : k = 0
271 1742 : DO i = 1, 10
272 1716 : IF (statepp%occ(l, i) /= 0._dp) THEN
273 46 : k = k + 1
274 46 : statepp%occupation(l, k) = state%occ(l, i)
275 46 : IF (state%multiplicity /= -1) THEN
276 4 : statepp%occa(l, k) = state%occa(l, i) - statepp%core(l, i)/2
277 4 : statepp%occb(l, k) = state%occb(l, i) - statepp%core(l, i)/2
278 : END IF
279 : END IF
280 : END DO
281 : END DO
282 :
283 26 : statepp%maxl_occ = get_maxl_occ(statepp%occ)
284 182 : statepp%maxn_occ = get_maxn_occ(statepp%occ)
285 26 : statepp%maxl_calc = state%maxl_calc
286 182 : statepp%maxn_calc = 0
287 26 : maxn = get_maxn_occ(statepp%core)
288 96 : DO k = 0, statepp%maxl_calc
289 70 : statepp%maxn_calc(k) = state%maxn_calc(k) - maxn(k)
290 96 : statepp%maxn_calc(k) = MIN(statepp%maxn_calc(k), pp_basis%nbas(k))
291 : END DO
292 26 : statepp%multiplicity = state%multiplicity
293 :
294 26 : CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_section=im)
295 26 : CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
296 26 : CALL set_atom(atom_info(in, im)%atom, method_type=method)
297 26 : CALL set_atom(atom_refs(in, im)%atom, method_type=method, relativistic=reltyp)
298 :
299 : ! calculate integrals: pseudopotential basis
300 : ! general integrals
301 26 : CALL atom_int_setup(pp_int, pp_basis, potential=p_pot, eri_coulomb=eri_c, eri_exchange=eri_e)
302 : !
303 26 : NULLIFY (pp_int%tzora, pp_int%hdkh)
304 : ! potential
305 26 : CALL atom_ppint_setup(pp_int, pp_basis, potential=p_pot)
306 : !
307 26 : CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, integrals=pp_int, potential=p_pot)
308 338 : statepp%maxn_calc(:) = MIN(statepp%maxn_calc(:), pp_basis%nbas(:))
309 182 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
310 :
311 : ! calculate integrals: all electron basis
312 : ! general integrals
313 : CALL atom_int_setup(ae_int, ae_basis, potential=ae_pot, &
314 26 : eri_coulomb=eri_c, eri_exchange=eri_e)
315 : ! potential
316 26 : CALL atom_ppint_setup(ae_int, ae_basis, potential=ae_pot)
317 : ! relativistic correction terms
318 26 : CALL atom_relint_setup(ae_int, ae_basis, reltyp, zcore=REAL(zval, dp))
319 : !
320 26 : CALL set_atom(atom_refs(in, im)%atom, basis=ae_basis, integrals=ae_int, potential=ae_pot)
321 338 : state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
322 182 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
323 :
324 : CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
325 26 : exchange_integral_type=do_erie)
326 : CALL set_atom(atom_refs(in, im)%atom, coulomb_integral_type=do_eric, &
327 26 : exchange_integral_type=do_erie)
328 26 : atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
329 26 : atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
330 26 : atom_refs(in, im)%atom%hfx_pot%do_gh = do_gh
331 26 : atom_refs(in, im)%atom%hfx_pot%nr_gh = nr_gh
332 :
333 26 : CALL set_atom(atom_info(in, im)%atom, state=statepp)
334 26 : NULLIFY (orbitals)
335 182 : mo = MAXVAL(statepp%maxn_calc)
336 182 : mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
337 26 : CALL create_atom_orbs(orbitals, mb, mo)
338 26 : CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
339 :
340 26 : CALL set_atom(atom_refs(in, im)%atom, state=state)
341 26 : NULLIFY (orbitals)
342 182 : mo = MAXVAL(state%maxn_calc)
343 182 : mb = MAXVAL(atom_refs(in, im)%atom%basis%nbas)
344 26 : CALL create_atom_orbs(orbitals, mb, mo)
345 26 : CALL set_atom(atom_refs(in, im)%atom, orbitals=orbitals)
346 :
347 78 : IF (atom_consistent_method(atom_refs(in, im)%atom%method_type, atom_refs(in, im)%atom%state%multiplicity)) THEN
348 : !Print method info
349 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
350 26 : CALL atom_print_method(atom_refs(in, im)%atom, iw)
351 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
352 : !Calculate the electronic structure
353 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
354 26 : CALL calculate_atom(atom_refs(in, im)%atom, iw)
355 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
356 : END IF
357 : END DO
358 : END DO
359 :
360 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_PSEUDO", extension=".log")
361 26 : IF (iw > 0) THEN
362 13 : WRITE (iw, '(/," ",79("*"))')
363 13 : WRITE (iw, '(" ",21("*"),A,21("*"))') " Optimize Pseudopotential Parameters "
364 13 : WRITE (iw, '(" ",79("*"))')
365 : END IF
366 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_PSEUDO")
367 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
368 26 : IF (iw > 0) THEN
369 0 : CALL atom_print_potential(p_pot, iw)
370 : END IF
371 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
372 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_PSEUDO", extension=".log")
373 26 : IF (iw > 0) THEN
374 13 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
375 13 : CALL atom_fit_pseudo(atom_info, atom_refs, p_pot, iw, powell_section)
376 : END IF
377 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%FIT_PSEUDO")
378 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
379 26 : IF (iw > 0) THEN
380 0 : CALL atom_print_potential(p_pot, iw)
381 : END IF
382 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
383 :
384 : ! Print out the orbitals if requested
385 26 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%ORBITALS", extension=".log")
386 26 : CALL section_vals_val_get(atom_section, "PRINT%ORBITALS%XMGRACE", l_val=graph)
387 26 : IF (iw > 0) THEN
388 0 : DO in = 1, n_rep
389 0 : DO im = 1, n_meth
390 0 : CALL atom_print_orbitals(atom_info(in, im)%atom, iw, xmgrace=graph)
391 : END DO
392 : END DO
393 : END IF
394 26 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%ORBITALS")
395 :
396 : ! clean up
397 26 : CALL atom_int_release(ae_int)
398 26 : CALL atom_ppint_release(ae_int)
399 26 : CALL atom_relint_release(ae_int)
400 :
401 26 : CALL atom_int_release(pp_int)
402 26 : CALL atom_ppint_release(pp_int)
403 26 : CALL atom_relint_release(pp_int)
404 :
405 26 : CALL release_atom_basis(ae_basis)
406 26 : CALL release_atom_basis(pp_basis)
407 :
408 26 : CALL release_atom_potential(p_pot)
409 26 : CALL release_atom_potential(ae_pot)
410 :
411 52 : DO in = 1, n_rep
412 78 : DO im = 1, n_meth
413 26 : CALL release_atom_type(atom_info(in, im)%atom)
414 52 : CALL release_atom_type(atom_refs(in, im)%atom)
415 : END DO
416 : END DO
417 26 : DEALLOCATE (atom_info, atom_refs)
418 :
419 26 : DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
420 :
421 26 : CALL timestop(handle)
422 :
423 182 : END SUBROUTINE atom_pseudo_opt
424 :
425 : ! **************************************************************************************************
426 :
427 : END MODULE atom_pseudo
|
