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 : !> \brief Interface to the SIRIUS Library
10 : !> \par History
11 : !> 07.2018 initial create
12 : !> \author JHU
13 : !***************************************************************************************************
14 : #if defined(__SIRIUS)
15 : MODULE sirius_interface
16 : USE ISO_C_BINDING, ONLY: C_DOUBLE,&
17 : C_INT,&
18 : C_LOC
19 : USE atom_kind_orbitals, ONLY: calculate_atomic_orbitals,&
20 : gth_potential_conversion
21 : USE atom_types, ONLY: atom_gthpot_type
22 : USE atom_upf, ONLY: atom_upfpot_type
23 : USE atom_utils, ONLY: atom_local_potential
24 : USE atomic_kind_types, ONLY: atomic_kind_type,&
25 : get_atomic_kind
26 : USE cell_types, ONLY: cell_type,&
27 : real_to_scaled
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_p_file,&
32 : cp_print_key_finished_output,&
33 : cp_print_key_should_output,&
34 : cp_print_key_unit_nr
35 : USE external_potential_types, ONLY: gth_potential_type
36 : USE input_constants, ONLY: do_gapw_log
37 : USE input_cp2k_pwdft, ONLY: SIRIUS_FUNC_VDWDF,&
38 : SIRIUS_FUNC_VDWDF2,&
39 : SIRIUS_FUNC_VDWDFCX
40 : USE input_section_types, ONLY: section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_get_subs_vals2,&
43 : section_vals_type,&
44 : section_vals_val_get
45 : USE kinds, ONLY: default_string_length,&
46 : dp
47 : USE machine, ONLY: m_flush
48 : USE mathconstants, ONLY: fourpi,&
49 : gamma1
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE particle_types, ONLY: particle_type
52 : USE physcon, ONLY: massunit
53 : USE pwdft_environment_types, ONLY: pwdft_energy_type,&
54 : pwdft_env_get,&
55 : pwdft_env_set,&
56 : pwdft_environment_type
57 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
58 : create_grid_atom,&
59 : deallocate_grid_atom,&
60 : grid_atom_type
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_subsys_types, ONLY: qs_subsys_get,&
64 : qs_subsys_type
65 : USE sirius, ONLY: &
66 : SIRIUS_INTEGER_ARRAY_TYPE, SIRIUS_INTEGER_TYPE, SIRIUS_LOGICAL_ARRAY_TYPE, &
67 : SIRIUS_LOGICAL_TYPE, SIRIUS_NUMBER_ARRAY_TYPE, SIRIUS_NUMBER_TYPE, &
68 : SIRIUS_STRING_ARRAY_TYPE, SIRIUS_STRING_TYPE, sirius_add_atom, sirius_add_atom_type, &
69 : sirius_add_atom_type_radial_function, sirius_add_xc_functional, sirius_context_handler, &
70 : sirius_create_context, sirius_create_ground_state, sirius_create_kset_from_grid, &
71 : sirius_finalize, sirius_find_ground_state, sirius_get_band_energies, &
72 : sirius_get_band_occupancies, sirius_get_energy, sirius_get_forces, &
73 : sirius_get_kpoint_properties, sirius_get_num_kpoints, sirius_get_parameters, &
74 : sirius_get_stress_tensor, sirius_ground_state_handler, sirius_import_parameters, &
75 : sirius_initialize, sirius_initialize_context, sirius_is_initialized, &
76 : sirius_kpoint_set_handler, sirius_option_get_info, sirius_option_get_section_length, &
77 : sirius_option_set, sirius_set_atom_position, sirius_set_atom_type_dion, &
78 : sirius_set_atom_type_hubbard, sirius_set_atom_type_radial_grid, &
79 : sirius_set_lattice_vectors, sirius_set_mpi_grid_dims, sirius_update_ground_state
80 : #include "./base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 :
84 : PRIVATE
85 :
86 : ! Global parameters
87 :
88 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'sirius_interface'
89 :
90 : ! Public subroutines
91 :
92 : PUBLIC :: cp_sirius_create_env, &
93 : cp_sirius_energy_force, &
94 : cp_sirius_finalize, &
95 : cp_sirius_init, &
96 : cp_sirius_is_initialized, &
97 : cp_sirius_update_context
98 :
99 : CONTAINS
100 :
101 : !***************************************************************************************************
102 : !> \brief Initialize the Sirius library
103 : !> \par Creation (07.2018, JHU)
104 : !> \author JHU
105 : ! **************************************************************************************************
106 20 : SUBROUTINE cp_sirius_init()
107 20 : CALL sirius_initialize(.FALSE.)
108 20 : END SUBROUTINE cp_sirius_init
109 :
110 : !***************************************************************************************************
111 : !> \brief Check the initialisation status of the Sirius library
112 : !> \return Return the initialisation status of the Sirius library as boolean
113 : !> \par Creation (03.12.2025, MK)
114 : !> \author Matthias Krack
115 : ! **************************************************************************************************
116 9284 : LOGICAL FUNCTION cp_sirius_is_initialized()
117 9284 : CALL sirius_is_initialized(cp_sirius_is_initialized)
118 9284 : END FUNCTION cp_sirius_is_initialized
119 :
120 : !***************************************************************************************************
121 : !> \brief Finalize the Sirius library
122 : !> \par Creation (07.2018, JHU)
123 : !> \author JHU
124 : ! **************************************************************************************************
125 20 : SUBROUTINE cp_sirius_finalize()
126 20 : CALL sirius_finalize(.FALSE., .FALSE., .FALSE.)
127 20 : END SUBROUTINE cp_sirius_finalize
128 :
129 : !***************************************************************************************************
130 : !> \brief ...
131 : !> \param pwdft_env ...
132 : !> \param
133 : !> \par History
134 : !> 07.2018 Create the Sirius environment
135 : !> \author JHU
136 : ! **************************************************************************************************
137 20 : SUBROUTINE cp_sirius_create_env(pwdft_env)
138 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
139 : #if defined(__SIRIUS)
140 :
141 : CHARACTER(len=2) :: element_symbol
142 : CHARACTER(len=default_string_length) :: label
143 : INTEGER :: i, ii, jj, iatom, ibeta, ifun, ikind, iwf, j, l, &
144 : n, ns, natom, nbeta, nbs, nkind, nmesh, &
145 : num_mag_dims, sirius_mpi_comm, vdw_func, nu, lu, output_unit
146 20 : INTEGER, DIMENSION(:), POINTER :: mpi_grid_dims
147 : INTEGER(KIND=C_INT), DIMENSION(3) :: k_grid, k_shift
148 20 : INTEGER, DIMENSION(:), POINTER :: kk
149 : LOGICAL :: up, use_ref_cell
150 : LOGICAL(4) :: use_so, use_symmetry, dft_plus_u_atom
151 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: fun
152 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: dion
153 : REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v1, v2
154 : REAL(KIND=dp) :: al, angle1, angle2, cval, focc, &
155 : magnetization, mass, pf, rl, zeff, alpha_u, beta_u, &
156 : J0_u, J_u, U_u, occ_u, u_minus_J, vnlp, vnlm
157 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: beta, corden, ef, fe, locpot, rc, rp
158 : REAL(KIND=dp), DIMENSION(3) :: vr, vs, j_t
159 20 : REAL(KIND=dp), DIMENSION(:), POINTER :: density
160 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: wavefunction, wfninfo
161 : TYPE(atom_gthpot_type), POINTER :: gth_atompot
162 : TYPE(atom_upfpot_type), POINTER :: upf_pot
163 20 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
164 : TYPE(atomic_kind_type), POINTER :: atomic_kind
165 : TYPE(cell_type), POINTER :: my_cell
166 : TYPE(mp_para_env_type), POINTER :: para_env
167 : TYPE(grid_atom_type), POINTER :: atom_grid
168 : TYPE(gth_potential_type), POINTER :: gth_potential
169 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
170 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
171 : TYPE(qs_subsys_type), POINTER :: qs_subsys
172 : TYPE(section_vals_type), POINTER :: pwdft_section, pwdft_sub_section, &
173 : xc_fun, xc_section
174 : TYPE(sirius_context_handler) :: sctx
175 : TYPE(sirius_ground_state_handler) :: gs_handler
176 : TYPE(sirius_kpoint_set_handler) :: ks_handler
177 :
178 0 : CPASSERT(ASSOCIATED(pwdft_env))
179 :
180 20 : output_unit = cp_logger_get_default_io_unit()
181 : ! create context of simulation
182 20 : CALL pwdft_env_get(pwdft_env, para_env=para_env)
183 20 : sirius_mpi_comm = para_env%get_handle()
184 20 : CALL sirius_create_context(sirius_mpi_comm, sctx)
185 :
186 : ! the "fun" starts.
187 :
188 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_section, xc_input=xc_section)
189 :
190 : CALL section_vals_val_get(pwdft_section, "ignore_convergence_failure", &
191 20 : l_val=pwdft_env%ignore_convergence_failure)
192 : ! cp2k should *have* a function that return all xc_functionals. Doing
193 : ! manually is prone to errors
194 :
195 20 : IF (ASSOCIATED(xc_section)) THEN
196 20 : ifun = 0
197 38 : DO
198 58 : ifun = ifun + 1
199 58 : xc_fun => section_vals_get_subs_vals2(xc_section, i_section=ifun)
200 58 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
201 : ! Here, we do not have to check whether the functional name starts with XC_
202 : ! because we only allow the shorter form w/o XC_
203 58 : CALL sirius_add_xc_functional(sctx, "XC_"//TRIM(xc_fun%section%name))
204 : END DO
205 : END IF
206 :
207 : ! import control section
208 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "control")
209 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
210 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "control")
211 20 : CALL section_vals_val_get(pwdft_sub_section, "mpi_grid_dims", i_vals=mpi_grid_dims)
212 : END IF
213 :
214 : ! import parameters section
215 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "parameters")
216 :
217 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
218 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "parameters")
219 20 : CALL section_vals_val_get(pwdft_sub_section, "ngridk", i_vals=kk)
220 20 : k_grid(1) = kk(1)
221 20 : k_grid(2) = kk(2)
222 20 : k_grid(3) = kk(3)
223 :
224 20 : CALL section_vals_val_get(pwdft_sub_section, "shiftk", i_vals=kk)
225 20 : k_shift(1) = kk(1)
226 20 : k_shift(2) = kk(2)
227 20 : k_shift(3) = kk(3)
228 20 : CALL section_vals_val_get(pwdft_sub_section, "num_mag_dims", i_val=num_mag_dims)
229 20 : CALL section_vals_val_get(pwdft_sub_section, "use_symmetry", l_val=use_symmetry)
230 20 : CALL section_vals_val_get(pwdft_sub_section, "so_correction", l_val=use_so)
231 :
232 : ! now check if van der walls corrections are needed
233 : vdw_func = -1
234 : #ifdef __LIBVDWXC
235 20 : CALL section_vals_val_get(pwdft_sub_section, "vdw_functional", i_val=vdw_func)
236 0 : SELECT CASE (vdw_func)
237 : CASE (SIRIUS_FUNC_VDWDF)
238 0 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF")
239 : CASE (SIRIUS_FUNC_VDWDF2)
240 0 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF2")
241 : CASE (SIRIUS_FUNC_VDWDFCX)
242 20 : CALL sirius_add_xc_functional(sctx, "XC_FUNC_VDWDF2")
243 : CASE default
244 : END SELECT
245 : #endif
246 :
247 : END IF
248 :
249 : ! import mixer section
250 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "mixer")
251 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
252 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "mixer")
253 : END IF
254 :
255 : ! import settings section
256 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "settings")
257 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
258 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "settings")
259 : END IF
260 :
261 : ! import solver section
262 20 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "iterative_solver")
263 20 : IF (ASSOCIATED(pwdft_sub_section)) THEN
264 20 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "iterative_solver")
265 : END IF
266 :
267 : #if defined(__SIRIUS_DFTD4)
268 : ! import dftd3 and dftd4 section
269 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd4")
270 : IF (ASSOCIATED(pwdft_sub_section)) THEN
271 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd4")
272 : END IF
273 :
274 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "dftd3")
275 : IF (ASSOCIATED(pwdft_sub_section)) THEN
276 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "dftd3")
277 : END IF
278 : #endif
279 :
280 : !
281 : ! uncomment these lines when nlcg is officially supported
282 : !
283 : #if defined(__SIRIUS_NLCG)
284 : ! import nlcg section
285 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "nlcg")
286 : IF (ASSOCIATED(pwdft_sub_section)) THEN
287 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "nlcg")
288 : END IF
289 : #endif
290 :
291 : #if defined(__SIRIUS_VCSQNM)
292 : pwdft_sub_section => section_vals_get_subs_vals(pwdft_section, "vcsqnm")
293 : IF (ASSOCIATED(pwdft_sub_section)) THEN
294 : CALL cp_sirius_fill_in_section(sctx, pwdft_sub_section, "vcsqnm")
295 : END IF
296 : #endif
297 :
298 : !CALL sirius_dump_runtime_setup(sctx, "runtime.json")
299 20 : CALL sirius_import_parameters(sctx, '{}')
300 :
301 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
302 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
303 20 : CALL qs_subsys_get(qs_subsys, cell=my_cell, use_ref_cell=use_ref_cell)
304 80 : a1(:) = my_cell%hmat(:, 1)
305 80 : a2(:) = my_cell%hmat(:, 2)
306 80 : a3(:) = my_cell%hmat(:, 3)
307 20 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
308 :
309 20 : IF (use_ref_cell) THEN
310 0 : CPWARN("SIRIUS| The specified CELL_REF will be ignored for PW_DFT calculations")
311 : END IF
312 :
313 : ! set up the atomic type definitions
314 : CALL qs_subsys_get(qs_subsys, &
315 : atomic_kind_set=atomic_kind_set, &
316 : qs_kind_set=qs_kind_set, &
317 20 : particle_set=particle_set)
318 20 : nkind = SIZE(atomic_kind_set)
319 46 : DO ikind = 1, nkind
320 : CALL get_atomic_kind(atomic_kind_set(ikind), &
321 26 : name=label, element_symbol=element_symbol, mass=mass)
322 26 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
323 26 : NULLIFY (upf_pot, gth_potential)
324 26 : CALL get_qs_kind(qs_kind_set(ikind), upf_potential=upf_pot, gth_potential=gth_potential)
325 :
326 26 : IF (ASSOCIATED(upf_pot)) THEN
327 : CALL sirius_add_atom_type(sctx, label, fname=upf_pot%filename, &
328 : symbol=element_symbol, &
329 20 : mass=REAL(mass/massunit, KIND=C_DOUBLE))
330 :
331 6 : ELSEIF (ASSOCIATED(gth_potential)) THEN
332 : !
333 6 : NULLIFY (atom_grid)
334 6 : CALL allocate_grid_atom(atom_grid)
335 6 : nmesh = 929
336 6 : atom_grid%nr = nmesh
337 6 : CALL create_grid_atom(atom_grid, nmesh, 1, 1, 0, do_gapw_log)
338 24 : ALLOCATE (rp(nmesh), fun(nmesh))
339 6 : IF (atom_grid%rad(1) < atom_grid%rad(nmesh)) THEN
340 : up = .TRUE.
341 : ELSE
342 : up = .FALSE.
343 : END IF
344 : IF (up) THEN
345 0 : rp(1:nmesh) = atom_grid%rad(1:nmesh)
346 : ELSE
347 5580 : DO i = 1, nmesh
348 5580 : rp(i) = atom_grid%rad(nmesh - i + 1)
349 : END DO
350 : END IF
351 : ! add new atom type
352 : CALL sirius_add_atom_type(sctx, label, &
353 : zn=NINT(zeff + 0.001d0), &
354 : symbol=element_symbol, &
355 : mass=REAL(mass/massunit, KIND=C_DOUBLE), &
356 6 : spin_orbit=gth_potential%soc)
357 : !
358 2916 : ALLOCATE (gth_atompot)
359 6 : CALL gth_potential_conversion(gth_potential, gth_atompot)
360 : ! set radial grid
361 5580 : fun(1:nmesh) = rp(1:nmesh)
362 6 : CALL sirius_set_atom_type_radial_grid(sctx, label, nmesh, fun(1))
363 : ! set beta-projectors
364 : ! GTH SOC uses the same projectors, SIRIUS can use the same or different projectors for l+1/2, l-1/2 (l > 0 l+1/2 l < 0 l-/2 )
365 24 : ALLOCATE (ef(nmesh), beta(nmesh))
366 6 : ibeta = 0
367 30 : DO l = 0, 3
368 24 : IF (gth_atompot%nl(l) == 0) CYCLE
369 14 : rl = gth_atompot%rcnl(l)
370 : ! we need to multiply by r so that data transferred to sirius are r \beta(r) not beta(r)
371 13020 : ef(1:nmesh) = EXP(-0.5_dp*rp(1:nmesh)*rp(1:nmesh)/(rl*rl))
372 50 : DO i = 1, gth_atompot%nl(l)
373 30 : pf = rl**(l + 0.5_dp*(4._dp*i - 1._dp))
374 30 : j = l + 2*i - 1
375 30 : pf = SQRT(2._dp)/(pf*SQRT(gamma1(j)))
376 27900 : beta(:) = pf*rp**(l + 2*i - 2)*ef
377 30 : ibeta = ibeta + 1
378 27900 : fun(1:nmesh) = beta(1:nmesh)*rp(1:nmesh)
379 : CALL sirius_add_atom_type_radial_function(sctx, label, &
380 30 : "beta", fun(1), nmesh, l=l)
381 : ! we double the number of beta projectors for SO and l>0
382 54 : IF (gth_atompot%soc .AND. l /= 0) THEN
383 : CALL sirius_add_atom_type_radial_function(sctx, label, &
384 10 : "beta", fun(1), nmesh, l=-l)
385 : END IF
386 : END DO
387 : END DO
388 6 : DEALLOCATE (ef, beta)
389 6 : nbeta = ibeta
390 :
391 : ! nonlocal PP matrix elements
392 6 : IF (gth_atompot%soc) THEN
393 2 : nbs = 2*nbeta - gth_atompot%nl(0)
394 8 : ALLOCATE (dion(nbs, nbs))
395 : ELSE
396 16 : ALLOCATE (dion(nbeta, nbeta))
397 : END IF
398 458 : dion = 0.0_dp
399 6 : IF (gth_atompot%soc) THEN
400 2 : ns = gth_atompot%nl(0)
401 2 : IF (ns /= 0) THEN
402 26 : dion(1:ns, 1:ns) = gth_atompot%hnl(1:ns, 1:ns, 0)
403 : END IF
404 8 : DO l = 1, 3
405 6 : IF (gth_atompot%nl(l) == 0) CYCLE
406 16 : DO i = 1, gth_atompot%nl(l)
407 14 : ii = ns + 2*SUM(gth_atompot%nl(1:l - 1))
408 10 : ii = ii + 2*(i - 1) + 1
409 42 : DO j = 1, gth_atompot%nl(l)
410 34 : jj = ns + 2*SUM(gth_atompot%nl(1:l - 1))
411 26 : jj = jj + 2*(j - 1) + 1
412 26 : vnlp = gth_atompot%hnl(i, j, l) + 0.5_dp*l*gth_atompot%knl(i, j, l)
413 26 : vnlm = gth_atompot%hnl(i, j, l) - 0.5_dp*(l + 1)*gth_atompot%knl(i, j, l)
414 26 : dion(ii, jj) = vnlp
415 36 : dion(ii + 1, jj + 1) = vnlm
416 : END DO
417 : END DO
418 : END DO
419 2 : CALL sirius_set_atom_type_dion(sctx, label, nbs, dion(1, 1))
420 : ELSE
421 20 : DO l = 0, 3
422 16 : IF (gth_atompot%nl(l) == 0) CYCLE
423 14 : ibeta = SUM(gth_atompot%nl(0:l - 1)) + 1
424 8 : i = ibeta + gth_atompot%nl(l) - 1
425 52 : dion(ibeta:i, ibeta:i) = gth_atompot%hnl(1:gth_atompot%nl(l), 1:gth_atompot%nl(l), l)
426 : END DO
427 4 : CALL sirius_set_atom_type_dion(sctx, label, nbeta, dion(1, 1))
428 : END IF
429 :
430 6 : DEALLOCATE (dion)
431 :
432 : ! set non-linear core correction
433 6 : IF (gth_atompot%nlcc) THEN
434 0 : ALLOCATE (corden(nmesh), fe(nmesh), rc(nmesh))
435 0 : corden(:) = 0.0_dp
436 0 : n = gth_atompot%nexp_nlcc
437 0 : DO i = 1, n
438 0 : al = gth_atompot%alpha_nlcc(i)
439 0 : rc(:) = rp(:)/al
440 0 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
441 0 : DO j = 1, gth_atompot%nct_nlcc(i)
442 0 : cval = gth_atompot%cval_nlcc(j, i)
443 0 : corden(:) = corden(:) + fe(:)*rc(:)**(2*j - 2)*cval
444 : END DO
445 : END DO
446 0 : fun(1:nmesh) = corden(1:nmesh)*rp(1:nmesh)
447 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_core", &
448 0 : fun(1), nmesh)
449 0 : DEALLOCATE (corden, fe, rc)
450 : END IF
451 :
452 : ! local potential
453 18 : ALLOCATE (locpot(nmesh))
454 5580 : locpot(:) = 0.0_dp
455 6 : CALL atom_local_potential(locpot, gth_atompot, rp)
456 5580 : fun(1:nmesh) = locpot(1:nmesh)
457 : CALL sirius_add_atom_type_radial_function(sctx, label, "vloc", &
458 6 : fun(1), nmesh)
459 6 : DEALLOCATE (locpot)
460 : !
461 6 : NULLIFY (density, wavefunction, wfninfo)
462 : CALL calculate_atomic_orbitals(atomic_kind_set(ikind), qs_kind_set(ikind), &
463 : density=density, wavefunction=wavefunction, &
464 6 : wfninfo=wfninfo, agrid=atom_grid)
465 :
466 : ! set the atomic radial functions
467 22 : DO iwf = 1, SIZE(wavefunction, 2)
468 16 : focc = wfninfo(1, iwf)
469 16 : l = NINT(wfninfo(2, iwf))
470 : ! we can not easily get the principal quantum number
471 16 : nu = -1
472 16 : IF (up) THEN
473 0 : fun(1:nmesh) = wavefunction(1:nmesh, iwf)*rp(i)
474 : ELSE
475 14880 : DO i = 1, nmesh
476 14880 : fun(i) = wavefunction(nmesh - i + 1, iwf)*rp(i)
477 : END DO
478 : END IF
479 : CALL sirius_add_atom_type_radial_function(sctx, &
480 : label, "ps_atomic_wf", &
481 22 : fun(1), nmesh, l=l, occ=REAL(focc, KIND=C_DOUBLE), n=nu)
482 : END DO
483 :
484 : ! set total charge density of a free atom (to compute initial rho(r))
485 6 : IF (up) THEN
486 0 : fun(1:nmesh) = fourpi*density(1:nmesh)*atom_grid%rad(1:nmesh)**2
487 : ELSE
488 5580 : DO i = 1, nmesh
489 5580 : fun(i) = fourpi*density(nmesh - i + 1)*atom_grid%rad(nmesh - i + 1)**2
490 : END DO
491 : END IF
492 : CALL sirius_add_atom_type_radial_function(sctx, label, "ps_rho_total", &
493 6 : fun(1), nmesh)
494 :
495 6 : IF (ASSOCIATED(density)) DEALLOCATE (density)
496 6 : IF (ASSOCIATED(wavefunction)) DEALLOCATE (wavefunction)
497 6 : IF (ASSOCIATED(wfninfo)) DEALLOCATE (wfninfo)
498 :
499 6 : CALL deallocate_grid_atom(atom_grid)
500 6 : DEALLOCATE (rp, fun)
501 6 : DEALLOCATE (gth_atompot)
502 : !
503 : ELSE
504 : CALL cp_abort(__LOCATION__, &
505 0 : "CP2K/SIRIUS: atomic kind needs UPF or GTH potential definition")
506 : END IF
507 :
508 : CALL get_qs_kind(qs_kind_set(ikind), &
509 : dft_plus_u_atom=dft_plus_u_atom, &
510 : l_of_dft_plus_u=lu, &
511 : n_of_dft_plus_u=nu, &
512 : u_minus_j_target=u_minus_j, &
513 : U_of_dft_plus_u=U_u, &
514 : J_of_dft_plus_u=J_u, &
515 : alpha_of_dft_plus_u=alpha_u, &
516 : beta_of_dft_plus_u=beta_u, &
517 : J0_of_dft_plus_u=J0_u, &
518 26 : occupation_of_dft_plus_u=occ_u)
519 :
520 72 : IF (dft_plus_u_atom) THEN
521 0 : IF (nu < 1) THEN
522 0 : CPABORT("CP2K/SIRIUS (hubbard): principal quantum number not specified")
523 : END IF
524 :
525 0 : IF (lu < 0) THEN
526 0 : CPABORT("CP2K/SIRIUS (hubbard): l can not be negative.")
527 : END IF
528 :
529 0 : IF (occ_u < 0.0) THEN
530 0 : CPABORT("CP2K/SIRIUS (hubbard): the occupation number can not be negative.")
531 : END IF
532 :
533 0 : j_t(:) = 0.0
534 0 : IF (ABS(u_minus_j) < 1e-8) THEN
535 0 : j_t(1) = J_u
536 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
537 0 : occ_u, U_u, j_t, alpha_u, beta_u, J0_u)
538 : ELSE
539 : CALL sirius_set_atom_type_hubbard(sctx, label, lu, nu, &
540 0 : occ_u, u_minus_j, j_t, alpha_u, beta_u, J0_u)
541 : END IF
542 : END IF
543 :
544 : END DO
545 :
546 : ! add atoms to the unit cell
547 : ! WARNING: sirius accepts only fractional coordinates;
548 20 : natom = SIZE(particle_set)
549 58 : DO iatom = 1, natom
550 152 : vr(1:3) = particle_set(iatom)%r(1:3)
551 38 : CALL real_to_scaled(vs, vr, my_cell)
552 38 : atomic_kind => particle_set(iatom)%atomic_kind
553 38 : ikind = atomic_kind%kind_number
554 38 : CALL get_atomic_kind(atomic_kind, name=label)
555 38 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, magnetization=magnetization)
556 : ! angle of magnetization might come from input Atom x y z mx my mz
557 : ! or as an angle?
558 : ! Answer : SIRIUS only accept the magnetization as mx, my, mz
559 38 : IF (num_mag_dims == 3) THEN
560 4 : angle1 = 0.0_dp
561 4 : angle2 = 0.0_dp
562 4 : v1(1) = magnetization*SIN(angle1)*COS(angle2)
563 4 : v1(2) = magnetization*SIN(angle1)*SIN(angle2)
564 4 : v1(3) = magnetization*COS(angle1)
565 : ELSE
566 34 : v1 = 0._dp
567 34 : v1(3) = magnetization
568 : END IF
569 38 : v2(1:3) = vs(1:3)
570 96 : CALL sirius_add_atom(sctx, label, v2(1), v1(1))
571 : END DO
572 :
573 20 : CALL sirius_set_mpi_grid_dims(sctx, 2, mpi_grid_dims)
574 :
575 : ! initialize global variables/indices/arrays/etc. of the simulation
576 20 : CALL sirius_initialize_context(sctx)
577 :
578 : ! strictly speaking the parameter use_symmetry is initialized at the
579 : ! beginning but it does no harm to do it that way
580 20 : IF (use_symmetry) THEN
581 16 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.TRUE., kset_handler=ks_handler)
582 : ELSE
583 4 : CALL sirius_create_kset_from_grid(sctx, k_grid(1), k_shift(1), use_symmetry=.FALSE., kset_handler=ks_handler)
584 : END IF
585 : ! create ground-state class
586 20 : CALL sirius_create_ground_state(ks_handler, gs_handler)
587 :
588 20 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler, ks_handler=ks_handler)
589 : #endif
590 40 : END SUBROUTINE cp_sirius_create_env
591 :
592 : !***************************************************************************************************
593 : !> \brief ...
594 : !> \param pwdft_env ...
595 : !> \param
596 : !> \par History
597 : !> 07.2018 Update the Sirius environment
598 : !> \author JHU
599 : ! **************************************************************************************************
600 20 : SUBROUTINE cp_sirius_update_context(pwdft_env)
601 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
602 :
603 : INTEGER :: iatom, natom
604 : REAL(KIND=C_DOUBLE), DIMENSION(3) :: a1, a2, a3, v2
605 : REAL(KIND=dp), DIMENSION(3) :: vr, vs
606 : TYPE(cell_type), POINTER :: my_cell
607 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
608 : TYPE(qs_subsys_type), POINTER :: qs_subsys
609 : TYPE(sirius_context_handler) :: sctx
610 : TYPE(sirius_ground_state_handler) :: gs_handler
611 :
612 0 : CPASSERT(ASSOCIATED(pwdft_env))
613 20 : CALL pwdft_env_get(pwdft_env, sctx=sctx, gs_handler=gs_handler)
614 :
615 : ! get current positions and lattice vectors
616 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
617 :
618 : ! lattice vectors of the unit cell should be in [a.u.] (length is in [a.u.])
619 20 : CALL qs_subsys_get(qs_subsys, cell=my_cell)
620 80 : a1(:) = my_cell%hmat(:, 1)
621 80 : a2(:) = my_cell%hmat(:, 2)
622 80 : a3(:) = my_cell%hmat(:, 3)
623 20 : CALL sirius_set_lattice_vectors(sctx, a1(1), a2(1), a3(1))
624 :
625 : ! new atomic positions
626 20 : CALL qs_subsys_get(qs_subsys, particle_set=particle_set)
627 20 : natom = SIZE(particle_set)
628 58 : DO iatom = 1, natom
629 152 : vr(1:3) = particle_set(iatom)%r(1:3)
630 38 : CALL real_to_scaled(vs, vr, my_cell)
631 38 : v2(1:3) = vs(1:3)
632 58 : CALL sirius_set_atom_position(sctx, iatom, v2(1))
633 : END DO
634 :
635 : ! update ground-state class
636 20 : CALL sirius_update_ground_state(gs_handler)
637 :
638 20 : CALL pwdft_env_set(pwdft_env, sctx=sctx, gs_handler=gs_handler)
639 :
640 20 : END SUBROUTINE cp_sirius_update_context
641 :
642 : ! **************************************************************************************************
643 : !> \brief ...
644 : !> \param sctx ...
645 : !> \param section ...
646 : !> \param section_name ...
647 : ! **************************************************************************************************
648 100 : SUBROUTINE cp_sirius_fill_in_section(sctx, section, section_name)
649 : TYPE(sirius_context_handler), INTENT(INOUT) :: sctx
650 : TYPE(section_vals_type), POINTER :: section
651 : CHARACTER(*), INTENT(in) :: section_name
652 :
653 : CHARACTER(len=256), TARGET :: option_name
654 : CHARACTER(len=4096) :: description, usage
655 100 : CHARACTER(len=80), DIMENSION(:), POINTER :: tmp
656 : CHARACTER(len=80), TARGET :: str
657 : INTEGER :: ctype, elem, ic, j
658 100 : INTEGER, DIMENSION(:), POINTER :: ivals
659 : INTEGER, TARGET :: enum_length, ival, length, &
660 : num_possible_values, number_of_options
661 : LOGICAL :: explicit
662 100 : LOGICAL, DIMENSION(:), POINTER :: lvals
663 : LOGICAL, TARGET :: found, lval
664 100 : REAL(kind=dp), DIMENSION(:), POINTER :: rvals
665 : REAL(kind=dp), TARGET :: rval
666 :
667 100 : NULLIFY (rvals)
668 100 : NULLIFY (ivals)
669 100 : CALL sirius_option_get_section_length(section_name, number_of_options)
670 :
671 2140 : DO elem = 1, number_of_options
672 2040 : option_name = ''
673 : CALL sirius_option_get_info(section_name, &
674 : elem, &
675 : option_name, &
676 : 256, &
677 : ctype, &
678 : num_possible_values, &
679 : enum_length, &
680 : description, &
681 : 4096, &
682 : usage, &
683 2040 : 4096)
684 2140 : IF ((option_name /= 'memory_usage') .AND. (option_name /= 'xc_functionals') .AND. (option_name /= 'vk')) THEN
685 2000 : CALL section_vals_val_get(section, option_name, explicit=found)
686 2000 : IF (found) THEN
687 160 : SELECT CASE (ctype)
688 : CASE (SIRIUS_INTEGER_TYPE)
689 160 : CALL section_vals_val_get(section, option_name, i_val=ival)
690 160 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ival))
691 : CASE (SIRIUS_NUMBER_TYPE)
692 124 : CALL section_vals_val_get(section, option_name, r_val=rval)
693 124 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rval))
694 : CASE (SIRIUS_LOGICAL_TYPE)
695 30 : CALL section_vals_val_get(section, option_name, l_val=lval)
696 30 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lval))
697 : CASE (SIRIUS_STRING_TYPE) ! string nightmare
698 124 : str = ''
699 124 : CALL section_vals_val_get(section, option_name, explicit=explicit, c_val=str)
700 124 : str = TRIM(ADJUSTL(str))
701 10044 : DO j = 1, LEN(str)
702 9920 : ic = ICHAR(str(j:j))
703 10044 : IF (ic >= 65 .AND. ic < 90) str(j:j) = CHAR(ic + 32)
704 : END DO
705 :
706 124 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), max_length=LEN_TRIM(str))
707 : CASE (SIRIUS_INTEGER_ARRAY_TYPE)
708 20 : CALL section_vals_val_get(section, option_name, i_vals=ivals)
709 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(ivals(1)), &
710 20 : max_length=num_possible_values)
711 : CASE (SIRIUS_NUMBER_ARRAY_TYPE)
712 0 : CALL section_vals_val_get(section, option_name, r_vals=rvals)
713 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(rvals(1)), &
714 0 : max_length=num_possible_values)
715 : CASE (SIRIUS_LOGICAL_ARRAY_TYPE)
716 0 : CALL section_vals_val_get(section, option_name, l_vals=lvals)
717 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(lvals(1)), &
718 0 : max_length=num_possible_values)
719 : CASE (SIRIUS_STRING_ARRAY_TYPE)
720 0 : CALL section_vals_val_get(section, option_name, explicit=explicit, n_rep_val=length)
721 458 : DO j = 1, length
722 0 : str = ''
723 0 : CALL section_vals_val_get(section, option_name, i_rep_val=j, explicit=explicit, c_vals=tmp)
724 0 : str = TRIM(ADJUSTL(tmp(j)))
725 : CALL sirius_option_set(sctx, section_name, option_name, ctype, C_LOC(str), &
726 0 : max_length=LEN_TRIM(str), append=.TRUE.)
727 : END DO
728 : CASE DEFAULT
729 : END SELECT
730 : END IF
731 : END IF
732 : END DO
733 100 : END SUBROUTINE cp_sirius_fill_in_section
734 :
735 : !***************************************************************************************************
736 : !> \brief ...
737 : !> \param pwdft_env ...
738 : !> \param calculate_forces ...
739 : !> \param calculate_stress_tensor ...
740 : !> \param
741 : !> \par History
742 : !> 07.2018 start the Sirius library
743 : !> \author JHU
744 : ! **************************************************************************************************
745 20 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress_tensor)
746 : TYPE(pwdft_environment_type), INTENT(INOUT), &
747 : POINTER :: pwdft_env
748 : LOGICAL, INTENT(IN) :: calculate_forces, calculate_stress_tensor
749 :
750 : INTEGER :: iw, n1, n2
751 : LOGICAL :: do_print, gs_converged
752 : REAL(KIND=C_DOUBLE) :: etotal
753 20 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:, :) :: cforces
754 : REAL(KIND=C_DOUBLE), DIMENSION(3, 3) :: cstress
755 : REAL(KIND=dp), DIMENSION(3, 3) :: stress
756 20 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
757 : TYPE(cp_logger_type), POINTER :: logger
758 : TYPE(pwdft_energy_type), POINTER :: energy
759 : TYPE(section_vals_type), POINTER :: print_section, pwdft_input
760 : TYPE(sirius_ground_state_handler) :: gs_handler
761 :
762 0 : CPASSERT(ASSOCIATED(pwdft_env))
763 :
764 20 : NULLIFY (logger)
765 20 : logger => cp_get_default_logger()
766 20 : iw = cp_logger_get_default_io_unit(logger)
767 :
768 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, gs_handler=gs_handler)
769 20 : CALL sirius_find_ground_state(gs_handler, converged=gs_converged)
770 :
771 20 : IF (gs_converged) THEN
772 20 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS: ground state is converged"
773 : ELSE
774 0 : IF (pwdft_env%ignore_convergence_failure) THEN
775 0 : IF (iw > 0) WRITE (iw, '(A)') "CP2K/SIRIUS Warning: ground state is not converged"
776 : ELSE
777 0 : CPABORT("CP2K/SIRIUS (ground state): SIRIUS did not converge.")
778 : END IF
779 : END IF
780 20 : IF (iw > 0) CALL m_flush(iw)
781 :
782 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, energy=energy)
783 : etotal = 0.0_C_DOUBLE
784 :
785 20 : CALL sirius_get_energy(gs_handler, 'band-gap', etotal)
786 20 : energy%band_gap = etotal
787 :
788 : etotal = 0.0_C_DOUBLE
789 20 : CALL sirius_get_energy(gs_handler, 'total', etotal)
790 20 : energy%etotal = etotal
791 :
792 : ! extract entropy (TS returned by sirius is always negative, sign
793 : ! convention in QE)
794 : etotal = 0.0_C_DOUBLE
795 20 : CALL sirius_get_energy(gs_handler, 'demet', etotal)
796 20 : energy%entropy = -etotal
797 :
798 20 : IF (calculate_forces) THEN
799 14 : CALL pwdft_env_get(pwdft_env=pwdft_env, forces=forces)
800 14 : n1 = SIZE(forces, 1)
801 14 : n2 = SIZE(forces, 2)
802 :
803 56 : ALLOCATE (cforces(n2, n1))
804 142 : cforces = 0.0_C_DOUBLE
805 14 : CALL sirius_get_forces(gs_handler, 'total', cforces)
806 : ! Sirius computes the forces but cp2k use the gradient everywhere
807 : ! so a minus sign is needed.
808 : ! note also that sirius and cp2k store the forces transpose to each other
809 : ! sirius : forces(coordinates, atoms)
810 : ! cp2k : forces(atoms, coordinates)
811 152 : forces = -TRANSPOSE(cforces(:, :))
812 14 : DEALLOCATE (cforces)
813 : END IF
814 :
815 20 : IF (calculate_stress_tensor) THEN
816 0 : cstress = 0.0_C_DOUBLE
817 0 : CALL sirius_get_stress_tensor(gs_handler, 'total', cstress)
818 0 : stress(1:3, 1:3) = cstress(1:3, 1:3)
819 0 : CALL pwdft_env_set(pwdft_env=pwdft_env, stress=stress)
820 : END IF
821 :
822 20 : CALL pwdft_env_get(pwdft_env=pwdft_env, pwdft_input=pwdft_input)
823 20 : print_section => section_vals_get_subs_vals(pwdft_input, "PRINT")
824 20 : CALL section_vals_get(print_section, explicit=do_print)
825 20 : IF (do_print) THEN
826 2 : CALL cp_sirius_print_results(pwdft_env, print_section)
827 : END IF
828 40 : END SUBROUTINE cp_sirius_energy_force
829 :
830 : !***************************************************************************************************
831 : !> \brief ...
832 : !> \param pwdft_env ...
833 : !> \param print_section ...
834 : !> \param
835 : !> \par History
836 : !> 12.2019 init
837 : !> \author JHU
838 : ! **************************************************************************************************
839 2 : SUBROUTINE cp_sirius_print_results(pwdft_env, print_section)
840 : TYPE(pwdft_environment_type), INTENT(INOUT), &
841 : POINTER :: pwdft_env
842 : TYPE(section_vals_type), POINTER :: print_section
843 :
844 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
845 : INTEGER :: i, ik, iounit, ispn, iterstep, iv, iw, &
846 : nbands, nhist, nkpts, nspins
847 : INTEGER(KIND=C_INT) :: cint
848 : LOGICAL :: append, dos, ionode
849 : REAL(KIND=C_DOUBLE) :: creal
850 2 : REAL(KIND=C_DOUBLE), ALLOCATABLE, DIMENSION(:) :: slist
851 : REAL(KIND=dp) :: de, e_fermi(2), emax, emin, eval
852 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: wkpt
853 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
854 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: energies, occupations
855 : TYPE(cp_logger_type), POINTER :: logger
856 : TYPE(sirius_context_handler) :: sctx
857 : TYPE(sirius_ground_state_handler) :: gs_handler
858 : TYPE(sirius_kpoint_set_handler) :: ks_handler
859 :
860 2 : NULLIFY (logger)
861 4 : logger => cp_get_default_logger()
862 : ionode = logger%para_env%is_source()
863 2 : iounit = cp_logger_get_default_io_unit(logger)
864 :
865 : ! Density of States
866 2 : dos = BTEST(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
867 2 : IF (dos) THEN
868 2 : CALL pwdft_env_get(pwdft_env, ks_handler=ks_handler)
869 2 : CALL pwdft_env_get(pwdft_env, gs_handler=gs_handler)
870 2 : CALL pwdft_env_get(pwdft_env, sctx=sctx)
871 :
872 2 : CALL section_vals_val_get(print_section, "DOS%DELTA_E", r_val=de)
873 2 : CALL section_vals_val_get(print_section, "DOS%APPEND", l_val=append)
874 :
875 2 : CALL sirius_get_num_kpoints(ks_handler, cint)
876 2 : nkpts = cint
877 2 : CALL sirius_get_parameters(sctx, num_bands=cint)
878 2 : nbands = cint
879 2 : CALL sirius_get_parameters(sctx, num_spins=cint)
880 2 : nspins = cint
881 2 : e_fermi(:) = 0.0_dp
882 10 : ALLOCATE (energies(nbands, nspins, nkpts))
883 442 : energies = 0.0_dp
884 8 : ALLOCATE (occupations(nbands, nspins, nkpts))
885 442 : occupations = 0.0_dp
886 6 : ALLOCATE (wkpt(nkpts))
887 6 : ALLOCATE (slist(nbands))
888 10 : DO ik = 1, nkpts
889 8 : CALL sirius_get_kpoint_properties(ks_handler, ik, creal)
890 10 : wkpt(ik) = creal
891 : END DO
892 10 : DO ik = 1, nkpts
893 26 : DO ispn = 1, nspins
894 16 : CALL sirius_get_band_energies(ks_handler, ik, ispn, slist)
895 432 : energies(1:nbands, ispn, ik) = slist(1:nbands)
896 16 : CALL sirius_get_band_occupancies(ks_handler, ik, ispn, slist)
897 440 : occupations(1:nbands, ispn, ik) = slist(1:nbands)
898 : END DO
899 : END DO
900 442 : emin = MINVAL(energies)
901 442 : emax = MAXVAL(energies)
902 2 : nhist = NINT((emax - emin)/de) + 1
903 16 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
904 16110 : hist = 0.0_dp
905 16110 : occval = 0.0_dp
906 16110 : ehist = 0.0_dp
907 :
908 10 : DO ik = 1, nkpts
909 26 : DO ispn = 1, nspins
910 440 : DO i = 1, nbands
911 416 : eval = energies(i, ispn, ik) - emin
912 416 : iv = NINT(eval/de) + 1
913 416 : CPASSERT((iv > 0) .AND. (iv <= nhist))
914 416 : hist(iv, ispn) = hist(iv, ispn) + wkpt(ik)
915 432 : occval(iv, ispn) = occval(iv, ispn) + wkpt(ik)*occupations(i, ispn, ik)
916 : END DO
917 : END DO
918 : END DO
919 16110 : hist = hist/REAL(nbands, KIND=dp)
920 8054 : DO i = 1, nhist
921 24158 : ehist(i, 1:nspins) = emin + (i - 1)*de
922 : END DO
923 :
924 2 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
925 2 : my_act = "WRITE"
926 2 : IF (append .AND. iterstep > 1) THEN
927 0 : my_pos = "APPEND"
928 : ELSE
929 2 : my_pos = "REWIND"
930 : END IF
931 :
932 : iw = cp_print_key_unit_nr(logger, print_section, "DOS", &
933 : extension=".dos", file_position=my_pos, file_action=my_act, &
934 2 : file_form="FORMATTED")
935 2 : IF (iw > 0) THEN
936 1 : IF (nspins == 2) THEN
937 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
938 1 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
939 1 : WRITE (UNIT=iw, FMT="(T2,A, A)") " Energy[a.u.] Alpha_Density Occupation", &
940 2 : " Beta_Density Occupation"
941 : ELSE
942 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
943 0 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
944 0 : WRITE (UNIT=iw, FMT="(T2,A)") " Energy[a.u.] Density Occupation"
945 : END IF
946 4027 : DO i = 1, nhist
947 4026 : eval = emin + (i - 1)*de
948 4027 : IF (nspins == 2) THEN
949 4026 : WRITE (UNIT=iw, FMT="(F15.8,4F15.4)") eval, hist(i, 1), occval(i, 1), &
950 8052 : hist(i, 2), occval(i, 2)
951 : ELSE
952 0 : WRITE (UNIT=iw, FMT="(F15.8,2F15.4)") eval, hist(i, 1), occval(i, 1)
953 : END IF
954 : END DO
955 : END IF
956 2 : CALL cp_print_key_finished_output(iw, logger, print_section, "DOS")
957 :
958 2 : DEALLOCATE (energies, occupations, wkpt, slist)
959 4 : DEALLOCATE (hist, occval, ehist)
960 : END IF
961 4 : END SUBROUTINE cp_sirius_print_results
962 :
963 : END MODULE sirius_interface
964 :
965 : #else
966 :
967 : !***************************************************************************************************
968 : !> \brief Empty implementation in case SIRIUS is not compiled in.
969 : !***************************************************************************************************
970 : MODULE sirius_interface
971 : USE pwdft_environment_types, ONLY: pwdft_environment_type
972 : #include "./base/base_uses.f90"
973 :
974 : IMPLICIT NONE
975 :
976 : PRIVATE
977 :
978 : ! Public subroutines
979 :
980 : PUBLIC :: cp_sirius_create_env, &
981 : cp_sirius_energy_force, &
982 : cp_sirius_finalize, &
983 : cp_sirius_init, &
984 : cp_sirius_is_initialized, &
985 : cp_sirius_update_context
986 :
987 : CONTAINS
988 :
989 : ! **************************************************************************************************
990 : !> \brief Empty implementation in case SIRIUS is not compiled in.
991 : ! **************************************************************************************************
992 : SUBROUTINE cp_sirius_init()
993 : END SUBROUTINE cp_sirius_init
994 :
995 : ! **************************************************************************************************
996 : !> \brief Return always .FALSE. because the Sirius library is not compiled in.
997 : !> \return Return the initialisation status of the Sirius library as boolean
998 : ! **************************************************************************************************
999 : LOGICAL FUNCTION cp_sirius_is_initialized()
1000 : cp_sirius_is_initialized = .FALSE.
1001 : END FUNCTION cp_sirius_is_initialized
1002 :
1003 : ! **************************************************************************************************
1004 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1005 : ! **************************************************************************************************
1006 : SUBROUTINE cp_sirius_finalize()
1007 : END SUBROUTINE cp_sirius_finalize
1008 :
1009 : ! **************************************************************************************************
1010 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1011 : !> \param pwdft_env ...
1012 : ! **************************************************************************************************
1013 : SUBROUTINE cp_sirius_create_env(pwdft_env)
1014 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1015 :
1016 : MARK_USED(pwdft_env)
1017 : CPABORT("Sirius library is missing")
1018 : END SUBROUTINE cp_sirius_create_env
1019 :
1020 : ! **************************************************************************************************
1021 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1022 : !> \param pwdft_env ...
1023 : !> \param calculate_forces ...
1024 : !> \param calculate_stress ...
1025 : ! **************************************************************************************************
1026 : SUBROUTINE cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
1027 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1028 : LOGICAL :: calculate_forces, calculate_stress
1029 :
1030 : MARK_USED(pwdft_env)
1031 : MARK_USED(calculate_forces)
1032 : MARK_USED(calculate_stress)
1033 : CPABORT("Sirius library is missing")
1034 : END SUBROUTINE cp_sirius_energy_force
1035 :
1036 : ! **************************************************************************************************
1037 : !> \brief Empty implementation in case SIRIUS is not compiled in.
1038 : !> \param pwdft_env ...
1039 : ! **************************************************************************************************
1040 : SUBROUTINE cp_sirius_update_context(pwdft_env)
1041 : TYPE(pwdft_environment_type), POINTER :: pwdft_env
1042 :
1043 : MARK_USED(pwdft_env)
1044 : CPABORT("Sirius library is missing")
1045 : END SUBROUTINE cp_sirius_update_context
1046 :
1047 : END MODULE sirius_interface
1048 :
1049 : #endif
|
