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_basis
9 : USE atom_fit, ONLY: atom_fit_basis
10 : USE atom_output, ONLY: atom_print_basis,&
11 : atom_print_info,&
12 : atom_print_method,&
13 : atom_print_potential
14 : USE atom_types, ONLY: &
15 : atom_basis_type, atom_integrals, atom_optimization_type, atom_orbitals, atom_p_type, &
16 : atom_potential_type, atom_state, create_atom_orbs, create_atom_type, init_atom_basis, &
17 : init_atom_potential, lmat, read_atom_opt_section, release_atom_basis, &
18 : release_atom_potential, release_atom_type, set_atom
19 : USE atom_utils, ONLY: atom_consistent_method,&
20 : atom_set_occupation,&
21 : get_maxl_occ,&
22 : get_maxn_occ
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
26 : cp_print_key_unit_nr
27 : USE input_constants, ONLY: do_analytic
28 : USE input_section_types, ONLY: section_vals_get,&
29 : section_vals_get_subs_vals,&
30 : section_vals_type,&
31 : section_vals_val_get
32 : USE kinds, ONLY: default_string_length,&
33 : dp
34 : USE periodic_table, ONLY: nelem,&
35 : ptable
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 : PUBLIC :: atom_basis_opt
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_basis'
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Optimize the atomic basis set.
48 : !> \param atom_section ATOM input section
49 : !> \par History
50 : !> * 04.2009 created starting from the subroutine atom_energy() [Juerg Hutter]
51 : ! **************************************************************************************************
52 10 : SUBROUTINE atom_basis_opt(atom_section)
53 : TYPE(section_vals_type), POINTER :: atom_section
54 :
55 : CHARACTER(len=*), PARAMETER :: routineN = 'atom_basis_opt'
56 :
57 : CHARACTER(LEN=2) :: elem
58 : CHARACTER(LEN=default_string_length), &
59 10 : DIMENSION(:), POINTER :: tmpstringlist
60 : INTEGER :: do_eric, do_erie, handle, i, im, in, &
61 : iunit, iw, k, maxl, mb, method, mo, &
62 : n_meth, n_rep, nr_gh, reltyp, zcore, &
63 : zval, zz
64 : INTEGER, DIMENSION(0:lmat) :: maxn
65 10 : INTEGER, DIMENSION(:), POINTER :: cn
66 : LOGICAL :: do_gh, eri_c, eri_e, had_ae, had_pp, &
67 : pp_calc
68 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: pocc
69 : TYPE(atom_basis_type), POINTER :: ae_basis, pp_basis
70 : TYPE(atom_integrals), POINTER :: ae_int, pp_int
71 : TYPE(atom_optimization_type) :: optimization
72 : TYPE(atom_orbitals), POINTER :: orbitals
73 10 : TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
74 : TYPE(atom_potential_type), POINTER :: ae_pot, p_pot
75 : TYPE(atom_state), POINTER :: state
76 : TYPE(cp_logger_type), POINTER :: logger
77 : TYPE(section_vals_type), POINTER :: basis_section, method_section, &
78 : opt_section, potential_section, &
79 : powell_section, xc_section
80 :
81 10 : CALL timeset(routineN, handle)
82 :
83 : ! What atom do we calculate
84 10 : CALL section_vals_val_get(atom_section, "ATOMIC_NUMBER", i_val=zval)
85 10 : CALL section_vals_val_get(atom_section, "ELEMENT", c_val=elem)
86 10 : zz = 0
87 222 : DO i = 1, nelem
88 222 : IF (ptable(i)%symbol == elem) THEN
89 : zz = i
90 : EXIT
91 : END IF
92 : END DO
93 10 : IF (zz /= 1) zval = zz
94 :
95 : ! read and set up inofrmation on the basis sets
96 370 : ALLOCATE (ae_basis, pp_basis)
97 10 : basis_section => section_vals_get_subs_vals(atom_section, "AE_BASIS")
98 10 : NULLIFY (ae_basis%grid)
99 10 : CALL init_atom_basis(ae_basis, basis_section, zval, "AE")
100 10 : NULLIFY (pp_basis%grid)
101 10 : basis_section => section_vals_get_subs_vals(atom_section, "PP_BASIS")
102 10 : CALL init_atom_basis(pp_basis, basis_section, zval, "PP")
103 :
104 : ! print general and basis set information
105 10 : logger => cp_get_default_logger()
106 10 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%PROGRAM_BANNER", extension=".log")
107 10 : IF (iw > 0) CALL atom_print_info(zval, "Atomic Basis Optimization", iw)
108 10 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%PROGRAM_BANNER")
109 :
110 : ! read and setup information on the pseudopotential
111 10 : NULLIFY (potential_section)
112 10 : potential_section => section_vals_get_subs_vals(atom_section, "POTENTIAL")
113 107570 : ALLOCATE (ae_pot, p_pot)
114 10 : CALL init_atom_potential(p_pot, potential_section, zval)
115 10 : CALL init_atom_potential(ae_pot, potential_section, -1)
116 :
117 : ! if the ERI's are calculated analytically, we have to precalculate them
118 10 : eri_c = .FALSE.
119 10 : CALL section_vals_val_get(atom_section, "COULOMB_INTEGRALS", i_val=do_eric)
120 10 : IF (do_eric == do_analytic) eri_c = .TRUE.
121 10 : eri_e = .FALSE.
122 10 : CALL section_vals_val_get(atom_section, "EXCHANGE_INTEGRALS", i_val=do_erie)
123 10 : IF (do_erie == do_analytic) eri_e = .TRUE.
124 10 : CALL section_vals_val_get(atom_section, "USE_GAUSS_HERMITE", l_val=do_gh)
125 10 : CALL section_vals_val_get(atom_section, "GRID_POINTS_GH", i_val=nr_gh)
126 :
127 : ! information on the states to be calculated
128 10 : CALL section_vals_val_get(atom_section, "MAX_ANGULAR_MOMENTUM", i_val=maxl)
129 10 : maxn = 0
130 10 : CALL section_vals_val_get(atom_section, "CALCULATE_STATES", i_vals=cn)
131 20 : DO in = 1, MIN(SIZE(cn), 4)
132 20 : maxn(in - 1) = cn(in)
133 : END DO
134 70 : DO in = 0, lmat
135 60 : maxn(in) = MIN(maxn(in), ae_basis%nbas(in))
136 70 : maxn(in) = MIN(maxn(in), pp_basis%nbas(in))
137 : END DO
138 :
139 : ! read optimization section
140 10 : opt_section => section_vals_get_subs_vals(atom_section, "OPTIMIZATION")
141 10 : CALL read_atom_opt_section(optimization, opt_section)
142 :
143 10 : had_ae = .FALSE.
144 10 : had_pp = .FALSE.
145 :
146 : ! Check for the total number of electron configurations to be calculated
147 10 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", n_rep_val=n_rep)
148 : ! Check for the total number of method types to be calculated
149 10 : method_section => section_vals_get_subs_vals(atom_section, "METHOD")
150 10 : CALL section_vals_get(method_section, n_repetition=n_meth)
151 :
152 : ! integrals
153 4250 : ALLOCATE (ae_int, pp_int)
154 :
155 60 : ALLOCATE (atom_info(n_rep, n_meth))
156 :
157 20 : DO in = 1, n_rep
158 30 : DO im = 1, n_meth
159 :
160 10 : NULLIFY (atom_info(in, im)%atom)
161 10 : CALL create_atom_type(atom_info(in, im)%atom)
162 :
163 10 : atom_info(in, im)%atom%optimization = optimization
164 :
165 10 : atom_info(in, im)%atom%z = zval
166 10 : xc_section => section_vals_get_subs_vals(method_section, "XC", i_rep_section=im)
167 10 : atom_info(in, im)%atom%xc_section => xc_section
168 :
169 3630 : ALLOCATE (state)
170 :
171 : ! get the electronic configuration
172 : CALL section_vals_val_get(atom_section, "ELECTRON_CONFIGURATION", i_rep_val=in, &
173 10 : c_vals=tmpstringlist)
174 :
175 : ! set occupations
176 10 : CALL atom_set_occupation(tmpstringlist, state%occ, state%occupation, state%multiplicity)
177 10 : state%maxl_occ = get_maxl_occ(state%occ)
178 70 : state%maxn_occ = get_maxn_occ(state%occ)
179 :
180 : ! set number of states to be calculated
181 10 : state%maxl_calc = MAX(maxl, state%maxl_occ)
182 10 : state%maxl_calc = MIN(lmat, state%maxl_calc)
183 70 : state%maxn_calc = 0
184 30 : DO k = 0, state%maxl_calc
185 30 : state%maxn_calc(k) = MAX(maxn(k), state%maxn_occ(k))
186 : END DO
187 :
188 : ! is there a pseudo potential
189 22 : pp_calc = ANY(INDEX(tmpstringlist(1:), "CORE") /= 0)
190 10 : IF (pp_calc) THEN
191 : ! get and set the core occupations
192 6 : CALL section_vals_val_get(atom_section, "CORE", c_vals=tmpstringlist)
193 6 : CALL atom_set_occupation(tmpstringlist, state%core, pocc)
194 426 : zcore = zval - NINT(SUM(state%core))
195 6 : CALL set_atom(atom_info(in, im)%atom, zcore=zcore, pp_calc=.TRUE.)
196 6 : had_pp = .TRUE.
197 6 : CALL set_atom(atom_info(in, im)%atom, basis=pp_basis, potential=p_pot)
198 78 : state%maxn_calc(:) = MIN(state%maxn_calc(:), pp_basis%nbas(:))
199 42 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
200 : ELSE
201 284 : state%core = 0._dp
202 4 : CALL set_atom(atom_info(in, im)%atom, zcore=zval, pp_calc=.FALSE.)
203 4 : had_ae = .TRUE.
204 4 : CALL set_atom(atom_info(in, im)%atom, basis=ae_basis, potential=ae_pot)
205 52 : state%maxn_calc(:) = MIN(state%maxn_calc(:), ae_basis%nbas(:))
206 28 : CPASSERT(ALL(state%maxn_calc(:) >= state%maxn_occ))
207 : END IF
208 :
209 10 : CALL section_vals_val_get(method_section, "METHOD_TYPE", i_val=method, i_rep_val=im)
210 10 : CALL section_vals_val_get(method_section, "RELATIVISTIC", i_val=reltyp, i_rep_section=im)
211 10 : CALL set_atom(atom_info(in, im)%atom, method_type=method, relativistic=reltyp)
212 10 : CALL set_atom(atom_info(in, im)%atom, state=state)
213 : CALL set_atom(atom_info(in, im)%atom, coulomb_integral_type=do_eric, &
214 10 : exchange_integral_type=do_erie)
215 10 : atom_info(in, im)%atom%hfx_pot%do_gh = do_gh
216 10 : atom_info(in, im)%atom%hfx_pot%nr_gh = nr_gh
217 :
218 10 : IF (atom_consistent_method(method, state%multiplicity)) THEN
219 10 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%METHOD_INFO", extension=".log")
220 10 : CALL atom_print_method(atom_info(in, im)%atom, iw)
221 10 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%METHOD_INFO")
222 10 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%POTENTIAL", extension=".log")
223 10 : IF (pp_calc) THEN
224 6 : IF (iw > 0) CALL atom_print_potential(p_pot, iw)
225 : ELSE
226 4 : IF (iw > 0) CALL atom_print_potential(ae_pot, iw)
227 : END IF
228 10 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%POTENTIAL")
229 : ELSE
230 0 : CPABORT("METHOD_TYPE and MULTIPLICITY are incompatible")
231 : END IF
232 :
233 10 : NULLIFY (orbitals)
234 70 : mo = MAXVAL(state%maxn_calc)
235 70 : mb = MAXVAL(atom_info(in, im)%atom%basis%nbas)
236 10 : CALL create_atom_orbs(orbitals, mb, mo)
237 30 : CALL set_atom(atom_info(in, im)%atom, orbitals=orbitals)
238 :
239 : END DO
240 : END DO
241 :
242 : ! Start the Optimization
243 10 : powell_section => section_vals_get_subs_vals(atom_section, "POWELL")
244 10 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%SCF_INFO", extension=".log")
245 10 : iunit = cp_print_key_unit_nr(logger, atom_section, "PRINT%FIT_BASIS", extension=".log")
246 10 : IF (had_ae) THEN
247 4 : pp_calc = .FALSE.
248 4 : CALL atom_fit_basis(atom_info, ae_basis, pp_calc, iunit, powell_section)
249 : END IF
250 10 : IF (had_pp) THEN
251 6 : pp_calc = .TRUE.
252 6 : CALL atom_fit_basis(atom_info, pp_basis, pp_calc, iunit, powell_section)
253 : END IF
254 10 : CALL cp_print_key_finished_output(iunit, logger, atom_section, "PRINT%FIT_BASIS")
255 10 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%SCF_INFO")
256 10 : iw = cp_print_key_unit_nr(logger, atom_section, "PRINT%BASIS_SET", extension=".log")
257 10 : IF (iw > 0) THEN
258 0 : CALL atom_print_basis(ae_basis, iw, " All Electron Basis")
259 0 : CALL atom_print_basis(pp_basis, iw, " Pseudopotential Basis")
260 : END IF
261 10 : CALL cp_print_key_finished_output(iw, logger, atom_section, "PRINT%BASIS_SET")
262 :
263 10 : CALL release_atom_basis(ae_basis)
264 10 : CALL release_atom_basis(pp_basis)
265 :
266 10 : CALL release_atom_potential(p_pot)
267 10 : CALL release_atom_potential(ae_pot)
268 :
269 20 : DO in = 1, n_rep
270 30 : DO im = 1, n_meth
271 20 : CALL release_atom_type(atom_info(in, im)%atom)
272 : END DO
273 : END DO
274 10 : DEALLOCATE (atom_info)
275 :
276 10 : DEALLOCATE (ae_pot, p_pot, ae_basis, pp_basis, ae_int, pp_int)
277 :
278 10 : CALL timestop(handle)
279 :
280 40 : END SUBROUTINE atom_basis_opt
281 :
282 : END MODULE atom_basis
|
