Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for a linear scaling quickstep SCF run based on the density
10 : !> matrix, with a focus on the interface between dm_ls_scf and qs
11 : !> \par History
12 : !> 2011.04 created [Joost VandeVondele]
13 : !> \author Joost VandeVondele
14 : ! **************************************************************************************************
15 : MODULE dm_ls_scf_qs
16 : USE atomic_kind_types, ONLY: atomic_kind_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
20 : dbcsr_distribution_get, dbcsr_distribution_hold, dbcsr_distribution_new, &
21 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_info, &
22 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
23 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
24 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_get_default_unit_nr,&
27 : cp_logger_type
28 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
29 : USE dm_ls_scf_types, ONLY: ls_cluster_atomic,&
30 : ls_cluster_molecular,&
31 : ls_mstruct_type,&
32 : ls_scf_env_type
33 : USE input_constants, ONLY: external_density_guess,&
34 : ls_cluster_atomic,&
35 : ls_cluster_molecular
36 : USE kinds, ONLY: default_string_length,&
37 : dp
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE particle_list_types, ONLY: particle_list_type
40 : USE particle_types, ONLY: particle_type
41 : USE pw_env_types, ONLY: pw_env_get,&
42 : pw_env_type
43 : USE pw_methods, ONLY: pw_zero
44 : USE pw_pool_types, ONLY: pw_pool_p_type,&
45 : pw_pool_type
46 : USE pw_types, ONLY: pw_c1d_gs_type,&
47 : pw_r3d_rs_type
48 : USE qs_atomic_block, ONLY: calculate_atomic_block_dm
49 : USE qs_collocate_density, ONLY: calculate_rho_elec
50 : USE qs_core_energies, ONLY: calculate_ptrace
51 : USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
52 : gspace_mixing_nr
53 : USE qs_energy_types, ONLY: qs_energy_type
54 : USE qs_environment_types, ONLY: get_qs_env,&
55 : qs_environment_type
56 : USE qs_external_density, ONLY: read_scf_guess_density
57 : USE qs_gspace_mixing, ONLY: gspace_mixing
58 : USE qs_harris_types, ONLY: harris_type
59 : USE qs_harris_utils, ONLY: harris_density_update
60 : USE qs_initial_guess, ONLY: calculate_mopac_dm
61 : USE qs_kind_types, ONLY: qs_kind_type
62 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
63 : USE qs_ks_types, ONLY: qs_ks_did_change,&
64 : qs_ks_env_type,&
65 : set_ks_env
66 : USE qs_mixing_utils, ONLY: charge_mixing_init,&
67 : mixing_allocate,&
68 : mixing_init
69 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
70 : USE qs_rho_atom_types, ONLY: rho_atom_type
71 : USE qs_rho_methods, ONLY: qs_rho_update_rho
72 : USE qs_rho_types, ONLY: qs_rho_get,&
73 : qs_rho_type
74 : USE qs_subsys_types, ONLY: qs_subsys_get,&
75 : qs_subsys_type
76 : USE scf_control_types, ONLY: scf_control_type
77 : USE tblite_interface, ONLY: tb_get_energy
78 : USE tblite_types, ONLY: tblite_type
79 : #include "./base/base_uses.f90"
80 :
81 : IMPLICIT NONE
82 :
83 : PRIVATE
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_qs'
86 :
87 : PUBLIC :: matrix_ls_create, matrix_qs_to_ls, matrix_ls_to_qs, ls_scf_init_qs, &
88 : ls_nonscf_ks, ls_nonscf_energy, ls_scf_dm_to_ks, ls_scf_qs_atomic_guess, &
89 : write_matrix_to_cube, rho_mixing_ls_init, matrix_decluster
90 :
91 : CONTAINS
92 :
93 : ! **************************************************************************************************
94 : !> \brief create a matrix for use (and as a template) in ls based on a qs template
95 : !> \param matrix_ls ...
96 : !> \param matrix_qs ...
97 : !> \param ls_mstruct ...
98 : !> \par History
99 : !> 2011.03 created [Joost VandeVondele]
100 : !> 2015.09 add support for PAO [Ole Schuett]
101 : !> \author Joost VandeVondele
102 : ! **************************************************************************************************
103 1010 : SUBROUTINE matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
104 : TYPE(dbcsr_type) :: matrix_ls, matrix_qs
105 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
106 :
107 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_ls_create'
108 :
109 : CHARACTER(len=default_string_length) :: name
110 : INTEGER :: handle, iatom, imol, jatom, natom, nmol
111 1010 : INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: atom_to_cluster, atom_to_cluster_primus, &
112 1010 : clustered_blk_sizes, primus_of_mol
113 1010 : INTEGER, DIMENSION(:), POINTER :: clustered_col_dist, clustered_row_dist, &
114 1010 : ls_blk_sizes, ls_col_dist, ls_row_dist
115 : TYPE(dbcsr_distribution_type) :: ls_dist, ls_dist_clustered
116 :
117 1010 : CALL timeset(routineN, handle)
118 :
119 : ! Defaults -----------------------------------------------------------------------------------
120 1010 : CALL dbcsr_get_info(matrix_qs, col_blk_size=ls_blk_sizes, distribution=ls_dist)
121 1010 : CALL dbcsr_distribution_hold(ls_dist)
122 1010 : CALL dbcsr_distribution_get(ls_dist, row_dist=ls_row_dist, col_dist=ls_col_dist)
123 :
124 : ! PAO ----------------------------------------------------------------------------------------
125 1010 : IF (ls_mstruct%do_pao) THEN
126 512 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=ls_blk_sizes)
127 : END IF
128 :
129 : ! Clustering ---------------------------------------------------------------------------------
130 1068 : SELECT CASE (ls_mstruct%cluster_type)
131 : CASE (ls_cluster_atomic)
132 : ! do nothing
133 : CASE (ls_cluster_molecular)
134 : ! create format of the clustered matrix
135 58 : CALL dbcsr_get_info(matrix_qs, nblkrows_total=natom)
136 526 : nmol = MAXVAL(ls_mstruct%atom_to_molecule)
137 174 : ALLOCATE (atom_to_cluster_primus(natom))
138 116 : ALLOCATE (atom_to_cluster(natom))
139 174 : ALLOCATE (primus_of_mol(nmol))
140 526 : DO iatom = 1, natom
141 468 : atom_to_cluster(iatom) = ls_mstruct%atom_to_molecule(iatom)
142 : ! the first atom of the molecule is the primus
143 : ! if the number of atoms per molecule is independent of system size, this is not a quadratic loop
144 : ! it assumes that all atoms of the molecule are consecutive.
145 1722 : DO jatom = iatom, 1, -1
146 1722 : IF (ls_mstruct%atom_to_molecule(jatom) == atom_to_cluster(iatom)) THEN
147 1254 : atom_to_cluster_primus(iatom) = jatom
148 : ELSE
149 : EXIT
150 : END IF
151 : END DO
152 526 : primus_of_mol(atom_to_cluster(iatom)) = atom_to_cluster_primus(iatom)
153 : END DO
154 :
155 : ! row
156 116 : ALLOCATE (clustered_row_dist(nmol))
157 188 : DO imol = 1, nmol
158 188 : clustered_row_dist(imol) = ls_row_dist(primus_of_mol(imol))
159 : END DO
160 :
161 : ! col
162 116 : ALLOCATE (clustered_col_dist(nmol))
163 188 : DO imol = 1, nmol
164 188 : clustered_col_dist(imol) = ls_col_dist(primus_of_mol(imol))
165 : END DO
166 :
167 116 : ALLOCATE (clustered_blk_sizes(nmol))
168 58 : clustered_blk_sizes = 0
169 526 : DO iatom = 1, natom
170 : clustered_blk_sizes(atom_to_cluster(iatom)) = clustered_blk_sizes(atom_to_cluster(iatom)) + &
171 526 : ls_blk_sizes(iatom)
172 : END DO
173 58 : ls_blk_sizes => clustered_blk_sizes ! redirect pointer
174 :
175 : ! create new distribution
176 : CALL dbcsr_distribution_new(ls_dist_clustered, &
177 : template=ls_dist, &
178 : row_dist=clustered_row_dist, &
179 : col_dist=clustered_col_dist, &
180 58 : reuse_arrays=.TRUE.)
181 58 : CALL dbcsr_distribution_release(ls_dist)
182 116 : ls_dist = ls_dist_clustered
183 :
184 : CASE DEFAULT
185 1010 : CPABORT("Unknown LS cluster type")
186 : END SELECT
187 :
188 : ! Create actual matrix -----------------------------------------------------------------------
189 1010 : CALL dbcsr_get_info(matrix_qs, name=name)
190 : CALL dbcsr_create(matrix_ls, &
191 : name=name, &
192 : dist=ls_dist, &
193 : matrix_type="S", &
194 : row_blk_size=ls_blk_sizes, &
195 1010 : col_blk_size=ls_blk_sizes)
196 1010 : CALL dbcsr_distribution_release(ls_dist)
197 1010 : CALL dbcsr_finalize(matrix_ls)
198 :
199 1010 : CALL timestop(handle)
200 :
201 2020 : END SUBROUTINE matrix_ls_create
202 :
203 : ! **************************************************************************************************
204 : !> \brief first link to QS, copy a QS matrix to LS matrix
205 : !> used to isolate QS style matrices from LS style
206 : !> will be useful for future features (e.g. precision, symmetry, blocking, ...)
207 : !> \param matrix_ls ...
208 : !> \param matrix_qs ...
209 : !> \param ls_mstruct ...
210 : !> \param covariant ...
211 : !> \par History
212 : !> 2010.10 created [Joost VandeVondele]
213 : !> 2015.09 add support for PAO [Ole Schuett]
214 : !> \author Joost VandeVondele
215 : ! **************************************************************************************************
216 28102 : SUBROUTINE matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
217 : TYPE(dbcsr_type) :: matrix_ls, matrix_qs
218 : TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
219 : LOGICAL, INTENT(IN) :: covariant
220 :
221 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_qs_to_ls'
222 :
223 : INTEGER :: handle
224 28102 : INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
225 : TYPE(dbcsr_type) :: matrix_pao, matrix_tmp
226 : TYPE(dbcsr_type), POINTER :: matrix_trafo
227 :
228 28102 : CALL timeset(routineN, handle)
229 :
230 28102 : IF (.NOT. ls_mstruct%do_pao) THEN
231 2658 : CALL matrix_cluster(matrix_ls, matrix_qs, ls_mstruct)
232 :
233 : ELSE ! using pao
234 25444 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
235 : CALL dbcsr_create(matrix_pao, &
236 : matrix_type="N", &
237 : template=matrix_qs, &
238 : row_blk_size=pao_blk_sizes, &
239 25444 : col_blk_size=pao_blk_sizes)
240 :
241 25444 : matrix_trafo => ls_mstruct%matrix_A ! contra-variant
242 25444 : IF (covariant) matrix_trafo => ls_mstruct%matrix_B ! co-variant
243 25444 : CALL dbcsr_create(matrix_tmp, template=matrix_trafo)
244 :
245 25444 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qs, matrix_trafo, 0.0_dp, matrix_tmp)
246 25444 : CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_trafo, matrix_tmp, 0.0_dp, matrix_pao)
247 25444 : CALL dbcsr_release(matrix_tmp)
248 :
249 25444 : CALL matrix_cluster(matrix_ls, matrix_pao, ls_mstruct)
250 25444 : CALL dbcsr_release(matrix_pao)
251 : END IF
252 :
253 28102 : CALL timestop(handle)
254 :
255 28102 : END SUBROUTINE matrix_qs_to_ls
256 :
257 : ! **************************************************************************************************
258 : !> \brief Performs molecular blocking and reduction to single precision if enabled
259 : !> \param matrix_out ...
260 : !> \param matrix_in ...
261 : !> \param ls_mstruct ...
262 : !> \author Ole Schuett
263 : ! **************************************************************************************************
264 28102 : SUBROUTINE matrix_cluster(matrix_out, matrix_in, ls_mstruct)
265 : TYPE(dbcsr_type) :: matrix_out, matrix_in
266 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
267 :
268 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_cluster'
269 :
270 : INTEGER :: handle
271 : TYPE(dbcsr_type) :: matrix_in_nosym
272 :
273 28102 : CALL timeset(routineN, handle)
274 :
275 54644 : SELECT CASE (ls_mstruct%cluster_type)
276 : CASE (ls_cluster_atomic)
277 26542 : CALL dbcsr_copy(matrix_out, matrix_in)
278 :
279 : CASE (ls_cluster_molecular)
280 : ! desymmetrize the qs matrix
281 1560 : CALL dbcsr_create(matrix_in_nosym, template=matrix_in, matrix_type="N")
282 1560 : CALL dbcsr_desymmetrize(matrix_in, matrix_in_nosym)
283 :
284 : ! perform the magic complete redistribute copy
285 1560 : CALL dbcsr_complete_redistribute(matrix_in_nosym, matrix_out)
286 1560 : CALL dbcsr_release(matrix_in_nosym)
287 :
288 : CASE DEFAULT
289 28102 : CPABORT("Unknown LS cluster type")
290 : END SELECT
291 :
292 28102 : CALL timestop(handle)
293 :
294 28102 : END SUBROUTINE matrix_cluster
295 :
296 : ! **************************************************************************************************
297 : !> \brief second link to QS, copy a LS matrix to QS matrix
298 : !> used to isolate QS style matrices from LS style
299 : !> will be useful for future features (e.g. precision, symmetry, blocking, ...)
300 : !> \param matrix_qs ...
301 : !> \param matrix_ls ...
302 : !> \param ls_mstruct ...
303 : !> \param covariant ...
304 : !> \param keep_sparsity will be passed on to dbcsr_copy, by default set to .TRUE.
305 : !> \par History
306 : !> 2010.10 created [Joost VandeVondele]
307 : !> 2015.09 add support for PAO [Ole Schuett]
308 : !> \author Joost VandeVondele
309 : ! **************************************************************************************************
310 4212 : SUBROUTINE matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
311 : TYPE(dbcsr_type) :: matrix_qs, matrix_ls
312 : TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
313 : LOGICAL :: covariant
314 : LOGICAL, OPTIONAL :: keep_sparsity
315 :
316 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_ls_to_qs'
317 :
318 : INTEGER :: handle
319 4212 : INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
320 : LOGICAL :: my_keep_sparsity
321 : TYPE(dbcsr_type) :: matrix_declustered, matrix_tmp1, &
322 : matrix_tmp2
323 : TYPE(dbcsr_type), POINTER :: matrix_trafo
324 :
325 4212 : CALL timeset(routineN, handle)
326 :
327 4212 : my_keep_sparsity = .TRUE.
328 4212 : IF (PRESENT(keep_sparsity)) THEN
329 294 : my_keep_sparsity = keep_sparsity
330 : END IF
331 :
332 4212 : IF (.NOT. ls_mstruct%do_pao) THEN
333 2544 : CALL dbcsr_create(matrix_declustered, template=matrix_qs)
334 2544 : CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
335 2544 : CALL dbcsr_copy(matrix_qs, matrix_declustered, keep_sparsity=my_keep_sparsity)
336 2544 : CALL dbcsr_release(matrix_declustered)
337 :
338 : ELSE ! using pao
339 1668 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
340 : CALL dbcsr_create(matrix_declustered, &
341 : template=matrix_qs, &
342 : row_blk_size=pao_blk_sizes, &
343 1668 : col_blk_size=pao_blk_sizes)
344 :
345 1668 : CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
346 :
347 1668 : matrix_trafo => ls_mstruct%matrix_B ! contra-variant
348 1668 : IF (covariant) matrix_trafo => ls_mstruct%matrix_A ! co-variant
349 1668 : CALL dbcsr_create(matrix_tmp1, template=matrix_trafo)
350 1668 : CALL dbcsr_create(matrix_tmp2, template=matrix_qs)
351 1668 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_declustered, 0.0_dp, matrix_tmp1)
352 1668 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp1, matrix_trafo, 0.0_dp, matrix_tmp2)
353 1668 : CALL dbcsr_copy(matrix_qs, matrix_tmp2, keep_sparsity=my_keep_sparsity)
354 1668 : CALL dbcsr_release(matrix_declustered)
355 1668 : CALL dbcsr_release(matrix_tmp1)
356 1668 : CALL dbcsr_release(matrix_tmp2)
357 : END IF
358 :
359 4212 : CALL timestop(handle)
360 :
361 4212 : END SUBROUTINE matrix_ls_to_qs
362 :
363 : ! **************************************************************************************************
364 : !> \brief Reverses molecular blocking and reduction to single precision if enabled
365 : !> \param matrix_out ...
366 : !> \param matrix_in ...
367 : !> \param ls_mstruct ...
368 : !> \author Ole Schuett
369 : ! **************************************************************************************************
370 12192 : SUBROUTINE matrix_decluster(matrix_out, matrix_in, ls_mstruct)
371 : TYPE(dbcsr_type) :: matrix_out, matrix_in
372 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
373 :
374 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_decluster'
375 :
376 : INTEGER :: handle
377 :
378 12192 : CALL timeset(routineN, handle)
379 :
380 23498 : SELECT CASE (ls_mstruct%cluster_type)
381 : CASE (ls_cluster_atomic)
382 11306 : CALL dbcsr_copy(matrix_out, matrix_in)
383 :
384 : CASE (ls_cluster_molecular)
385 : ! perform the magic complete redistribute copy
386 886 : CALL dbcsr_complete_redistribute(matrix_in, matrix_out)
387 :
388 : CASE DEFAULT
389 12192 : CPABORT("Unknown LS cluster type")
390 : END SELECT
391 :
392 12192 : CALL timestop(handle)
393 :
394 12192 : END SUBROUTINE matrix_decluster
395 :
396 : ! **************************************************************************************************
397 : !> \brief further required initialization of QS.
398 : !> Might be factored-out since this seems common code with the other SCF.
399 : !> \param qs_env ...
400 : !> \par History
401 : !> 2010.10 created [Joost VandeVondele]
402 : !> \author Joost VandeVondele
403 : ! **************************************************************************************************
404 980 : SUBROUTINE ls_scf_init_qs(qs_env)
405 : TYPE(qs_environment_type), POINTER :: qs_env
406 :
407 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_qs'
408 :
409 : INTEGER :: handle, ispin, nspin, unit_nr
410 : TYPE(cp_logger_type), POINTER :: logger
411 980 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
412 : TYPE(dft_control_type), POINTER :: dft_control
413 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
414 980 : POINTER :: sab_orb
415 : TYPE(qs_ks_env_type), POINTER :: ks_env
416 :
417 980 : NULLIFY (sab_orb)
418 980 : CALL timeset(routineN, handle)
419 :
420 : ! get a useful output_unit
421 980 : logger => cp_get_default_logger()
422 980 : IF (logger%para_env%is_source()) THEN
423 490 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
424 : ELSE
425 : unit_nr = -1
426 : END IF
427 :
428 : ! get basic quantities from the qs_env
429 : CALL get_qs_env(qs_env, dft_control=dft_control, &
430 : matrix_s=matrix_s, &
431 : matrix_ks=matrix_ks, &
432 : ks_env=ks_env, &
433 980 : sab_orb=sab_orb)
434 :
435 980 : nspin = dft_control%nspins
436 :
437 : ! we might have to create matrix_ks
438 980 : IF (.NOT. ASSOCIATED(matrix_ks)) THEN
439 0 : CALL dbcsr_allocate_matrix_set(matrix_ks, nspin)
440 0 : DO ispin = 1, nspin
441 0 : ALLOCATE (matrix_ks(ispin)%matrix)
442 0 : CALL dbcsr_create(matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
443 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix, sab_orb)
444 0 : CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
445 : END DO
446 0 : CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
447 : END IF
448 :
449 980 : CALL timestop(handle)
450 :
451 980 : END SUBROUTINE ls_scf_init_qs
452 :
453 : ! **************************************************************************************************
454 : !> \brief get an atomic initial guess
455 : !> \param qs_env ...
456 : !> \param ls_scf_env ...
457 : !> \param energy ...
458 : !> \param nonscf ...
459 : !> \par History
460 : !> 2012.11 created [Joost VandeVondele]
461 : !> \author Joost VandeVondele
462 : ! **************************************************************************************************
463 358 : SUBROUTINE ls_scf_qs_atomic_guess(qs_env, ls_scf_env, energy, nonscf)
464 : TYPE(qs_environment_type), POINTER :: qs_env
465 : TYPE(ls_scf_env_type) :: ls_scf_env
466 : REAL(KIND=dp) :: energy
467 : LOGICAL, INTENT(IN), OPTIONAL :: nonscf
468 :
469 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_qs_atomic_guess'
470 :
471 : INTEGER :: handle, nspin, unit_nr
472 : INTEGER, DIMENSION(2) :: nelectron_spin
473 : LOGICAL :: do_scf, has_unit_metric
474 358 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
475 : TYPE(cp_logger_type), POINTER :: logger
476 358 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
477 : TYPE(dft_control_type), POINTER :: dft_control
478 : TYPE(mp_para_env_type), POINTER :: para_env
479 358 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
480 : TYPE(qs_energy_type), POINTER :: qs_energy
481 358 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
482 : TYPE(qs_ks_env_type), POINTER :: ks_env
483 : TYPE(qs_rho_type), POINTER :: rho
484 : TYPE(scf_control_type), POINTER :: scf_control
485 :
486 358 : CALL timeset(routineN, handle)
487 358 : NULLIFY (rho, rho_ao, scf_control)
488 :
489 : ! get a useful output_unit
490 358 : logger => cp_get_default_logger()
491 358 : IF (logger%para_env%is_source()) THEN
492 179 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
493 : ELSE
494 179 : unit_nr = -1
495 : END IF
496 :
497 : ! get basic quantities from the qs_env
498 : CALL get_qs_env(qs_env, dft_control=dft_control, &
499 : matrix_s=matrix_s, &
500 : matrix_ks=matrix_ks, &
501 : ks_env=ks_env, &
502 : energy=qs_energy, &
503 : atomic_kind_set=atomic_kind_set, &
504 : qs_kind_set=qs_kind_set, &
505 : particle_set=particle_set, &
506 : has_unit_metric=has_unit_metric, &
507 : para_env=para_env, &
508 : nelectron_spin=nelectron_spin, &
509 : rho=rho, &
510 358 : scf_control=scf_control)
511 :
512 358 : CALL qs_rho_get(rho, rho_ao=rho_ao)
513 :
514 358 : nspin = dft_control%nspins
515 358 : IF (qs_env%harris_method .AND. scf_control%density_guess == external_density_guess) THEN
516 0 : CPABORT("EXTERNAL_DENSITY SCF guess is incompatible with HARRIS_METHOD")
517 : END IF
518 :
519 : ! create an initial atomic guess
520 358 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
521 : dft_control%qs_control%xtb) THEN
522 : CALL calculate_mopac_dm(rho_ao, matrix_s(1)%matrix, has_unit_metric, &
523 : dft_control, particle_set, atomic_kind_set, qs_kind_set, &
524 94 : nspin, nelectron_spin, para_env)
525 : ELSE
526 : CALL calculate_atomic_block_dm(rho_ao, matrix_s(1)%matrix, atomic_kind_set, qs_kind_set, &
527 264 : nspin, nelectron_spin, unit_nr, para_env)
528 : END IF
529 :
530 358 : do_scf = .TRUE.
531 358 : IF (PRESENT(nonscf)) do_scf = .NOT. nonscf
532 288 : IF (do_scf) THEN
533 346 : IF (scf_control%density_guess == external_density_guess) THEN
534 2 : CALL read_scf_guess_density(qs_env, scf_control%external_density_file_name)
535 : ELSE
536 344 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
537 : END IF
538 346 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
539 346 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
540 346 : CALL ls_scf_tblite_energy(qs_env, qs_energy)
541 346 : energy = qs_energy%total
542 : ELSE
543 12 : CALL ls_nonscf_ks(qs_env, ls_scf_env, energy)
544 : END IF
545 :
546 358 : CALL timestop(handle)
547 :
548 358 : END SUBROUTINE ls_scf_qs_atomic_guess
549 :
550 : ! **************************************************************************************************
551 : !> \brief use the density matrix in ls_scf_env to compute the new energy and KS matrix
552 : !> \param qs_env ...
553 : !> \param ls_scf_env ...
554 : !> \param energy_new ...
555 : !> \param iscf ...
556 : !> \par History
557 : !> 2011.04 created [Joost VandeVondele]
558 : !> 2015.02 added gspace density mixing [Patrick Seewald]
559 : !> \author Joost VandeVondele
560 : ! **************************************************************************************************
561 3368 : SUBROUTINE ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
562 : TYPE(qs_environment_type), POINTER :: qs_env
563 : TYPE(ls_scf_env_type) :: ls_scf_env
564 : REAL(KIND=dp) :: energy_new
565 : INTEGER, INTENT(IN) :: iscf
566 :
567 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_dm_to_ks'
568 :
569 : INTEGER :: handle, ispin, nspin, unit_nr
570 : TYPE(cp_logger_type), POINTER :: logger
571 3368 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
572 : TYPE(mp_para_env_type), POINTER :: para_env
573 : TYPE(qs_energy_type), POINTER :: energy
574 : TYPE(qs_rho_type), POINTER :: rho
575 :
576 3368 : NULLIFY (energy, rho, rho_ao)
577 3368 : CALL timeset(routineN, handle)
578 :
579 3368 : logger => cp_get_default_logger()
580 3368 : IF (logger%para_env%is_source()) THEN
581 1684 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
582 : ELSE
583 : unit_nr = -1
584 : END IF
585 :
586 3368 : nspin = ls_scf_env%nspins
587 3368 : CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
588 3368 : CALL qs_rho_get(rho, rho_ao=rho_ao)
589 :
590 : ! set the new density matrix
591 6878 : DO ispin = 1, nspin
592 : CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
593 6878 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
594 : END DO
595 :
596 : ! compute the corresponding KS matrix and new energy, mix density if requested
597 3368 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
598 3368 : IF (ls_scf_env%do_rho_mixing) THEN
599 0 : IF (ls_scf_env%density_mixing_method == direct_mixing_nr) THEN
600 0 : CPABORT("Direct P mixing not implemented in linear scaling SCF. ")
601 : END IF
602 0 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
603 0 : IF (iscf > MAX(ls_scf_env%mixing_store%nskip_mixing, 1)) THEN
604 : CALL gspace_mixing(qs_env, ls_scf_env%density_mixing_method, &
605 : ls_scf_env%mixing_store, rho, para_env, &
606 0 : iscf - 1)
607 0 : IF (unit_nr > 0) THEN
608 : WRITE (unit_nr, '(A57)') &
609 0 : "*********************************************************"
610 : WRITE (unit_nr, '(A13,F5.3,A20,A6,A7,I3)') &
611 0 : " Using ALPHA=", ls_scf_env%mixing_store%alpha, &
612 0 : " to mix rho: method=", ls_scf_env%mixing_store%iter_method, ", iscf=", iscf
613 : WRITE (unit_nr, '(A8,F5.3,A6,F5.3,A8)') &
614 0 : " rho_nw=", ls_scf_env%mixing_store%alpha, "*rho + ", &
615 0 : 1.0_dp - ls_scf_env%mixing_store%alpha, "*rho_old"
616 : WRITE (unit_nr, '(A57)') &
617 0 : "*********************************************************"
618 : END IF
619 : END IF
620 : END IF
621 : END IF
622 :
623 3368 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
624 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
625 3368 : just_energy=.FALSE., print_active=.TRUE.)
626 3368 : CALL ls_scf_tblite_energy(qs_env, energy)
627 3368 : energy_new = energy%total
628 :
629 3368 : CALL timestop(handle)
630 :
631 3368 : END SUBROUTINE ls_scf_dm_to_ks
632 :
633 : ! **************************************************************************************************
634 : !> \brief use the external density in ls_scf_env to compute the new KS matrix
635 : !> \param qs_env ...
636 : !> \param ls_scf_env ...
637 : !> \param energy_new ...
638 : ! **************************************************************************************************
639 72 : SUBROUTINE ls_nonscf_ks(qs_env, ls_scf_env, energy_new)
640 : TYPE(qs_environment_type), POINTER :: qs_env
641 : TYPE(ls_scf_env_type) :: ls_scf_env
642 : REAL(KIND=dp) :: energy_new
643 :
644 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_nonscf_ks'
645 :
646 : INTEGER :: handle, ispin, nspin
647 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
648 : TYPE(harris_type), POINTER :: harris_env
649 : TYPE(mp_para_env_type), POINTER :: para_env
650 : TYPE(qs_energy_type), POINTER :: energy
651 : TYPE(qs_rho_type), POINTER :: rho
652 :
653 72 : NULLIFY (energy, rho, rho_ao)
654 72 : CALL timeset(routineN, handle)
655 :
656 72 : nspin = ls_scf_env%nspins
657 72 : CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
658 72 : CALL qs_rho_get(rho, rho_ao=rho_ao)
659 :
660 : ! set the new density matrix
661 144 : DO ispin = 1, nspin
662 : CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
663 144 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
664 : END DO
665 :
666 72 : IF (qs_env%harris_method) THEN
667 20 : CALL get_qs_env(qs_env, harris_env=harris_env)
668 20 : CALL harris_density_update(qs_env, harris_env)
669 : END IF
670 : ! compute the corresponding KS matrix and new energy
671 72 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
672 72 : IF (ls_scf_env%do_rho_mixing) THEN
673 0 : CPABORT("P mixing not implemented in linear scaling NONSCF. ")
674 : END IF
675 :
676 72 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
677 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
678 72 : just_energy=.FALSE., print_active=.TRUE.)
679 72 : CALL ls_scf_tblite_energy(qs_env, energy)
680 72 : energy_new = energy%total
681 :
682 72 : CALL timestop(handle)
683 :
684 72 : END SUBROUTINE ls_nonscf_ks
685 :
686 : ! **************************************************************************************************
687 : !> \brief use the new density matrix in ls_scf_env to compute the new energy
688 : !> \param qs_env ...
689 : !> \param ls_scf_env ...
690 : ! **************************************************************************************************
691 72 : SUBROUTINE ls_nonscf_energy(qs_env, ls_scf_env)
692 : TYPE(qs_environment_type), POINTER :: qs_env
693 : TYPE(ls_scf_env_type) :: ls_scf_env
694 :
695 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_nonscf_energy'
696 :
697 : INTEGER :: handle, ispin, nspin
698 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao
699 : TYPE(mp_para_env_type), POINTER :: para_env
700 : TYPE(qs_energy_type), POINTER :: energy
701 : TYPE(qs_rho_type), POINTER :: rho
702 :
703 72 : NULLIFY (energy, rho, rho_ao)
704 72 : CALL timeset(routineN, handle)
705 72 : IF (qs_env%qmmm) THEN
706 0 : CPABORT("NYA")
707 : END IF
708 :
709 72 : nspin = ls_scf_env%nspins
710 72 : CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
711 72 : CALL qs_rho_get(rho, rho_ao=rho_ao)
712 :
713 : ! set the new density matrix
714 144 : DO ispin = 1, nspin
715 : CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
716 144 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
717 : END DO
718 :
719 72 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
720 :
721 : ! band energy : Tr(PH)
722 72 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
723 72 : CALL calculate_ptrace(matrix_ks, rho_ao, energy%band, nspin, .TRUE.)
724 : ! core energy : Tr(Ph)
725 72 : energy%total = energy%total - energy%core
726 72 : CALL get_qs_env(qs_env, matrix_h=matrix_h)
727 72 : CALL calculate_ptrace(matrix_h, rho_ao, energy%core, nspin)
728 :
729 72 : CALL timestop(handle)
730 :
731 72 : END SUBROUTINE ls_nonscf_energy
732 :
733 : ! **************************************************************************************************
734 : !> \brief update CP2K/tblite total energy after an LS_SCF KS rebuild.
735 : !> \param qs_env ...
736 : !> \param energy ...
737 : ! **************************************************************************************************
738 3786 : SUBROUTINE ls_scf_tblite_energy(qs_env, energy)
739 : TYPE(qs_environment_type), POINTER :: qs_env
740 : TYPE(qs_energy_type), POINTER :: energy
741 :
742 : TYPE(dft_control_type), POINTER :: dft_control
743 : TYPE(tblite_type), POINTER :: tb
744 :
745 3786 : NULLIFY (dft_control, tb)
746 3786 : CALL get_qs_env(qs_env, dft_control=dft_control, tb_tblite=tb)
747 :
748 3786 : IF (dft_control%qs_control%xtb .AND. dft_control%qs_control%xtb_control%do_tblite) THEN
749 80 : CPASSERT(ASSOCIATED(tb))
750 80 : CALL tb_get_energy(qs_env, tb, energy)
751 : END IF
752 :
753 3786 : END SUBROUTINE ls_scf_tblite_energy
754 :
755 : ! **************************************************************************************************
756 : !> \brief ...
757 : !> \param qs_env ...
758 : !> \param ls_scf_env ...
759 : !> \param matrix_p_ls ...
760 : !> \param unit_nr ...
761 : !> \param title ...
762 : !> \param stride ...
763 : ! **************************************************************************************************
764 6 : SUBROUTINE write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
765 : TYPE(qs_environment_type), POINTER :: qs_env
766 : TYPE(ls_scf_env_type) :: ls_scf_env
767 : TYPE(dbcsr_type), INTENT(IN) :: matrix_p_ls
768 : INTEGER, INTENT(IN) :: unit_nr
769 : CHARACTER(LEN=*), INTENT(IN) :: title
770 : INTEGER, DIMENSION(:), POINTER :: stride
771 :
772 : CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_to_cube'
773 :
774 : INTEGER :: handle
775 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
776 : TYPE(dbcsr_type), TARGET :: matrix_p_qs
777 : TYPE(particle_list_type), POINTER :: particles
778 : TYPE(pw_c1d_gs_type) :: wf_g
779 : TYPE(pw_env_type), POINTER :: pw_env
780 6 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
781 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
782 : TYPE(pw_r3d_rs_type) :: wf_r
783 : TYPE(qs_ks_env_type), POINTER :: ks_env
784 : TYPE(qs_subsys_type), POINTER :: subsys
785 :
786 6 : CALL timeset(routineN, handle)
787 :
788 6 : NULLIFY (ks_env, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, matrix_ks)
789 :
790 : CALL get_qs_env(qs_env, &
791 : ks_env=ks_env, &
792 : subsys=subsys, &
793 : pw_env=pw_env, &
794 6 : matrix_ks=matrix_ks)
795 :
796 6 : CALL qs_subsys_get(subsys, particles=particles)
797 :
798 : ! convert the density matrix (ls style) to QS style
799 6 : CALL dbcsr_copy(matrix_p_qs, matrix_ks(1)%matrix)
800 6 : CALL dbcsr_set(matrix_p_qs, 0.0_dp) !zero matrix creation
801 6 : CALL matrix_ls_to_qs(matrix_p_qs, matrix_p_ls, ls_scf_env%ls_mstruct, covariant=.FALSE.)
802 :
803 : ! Print total electronic density
804 : CALL pw_env_get(pw_env=pw_env, &
805 : auxbas_pw_pool=auxbas_pw_pool, &
806 6 : pw_pools=pw_pools)
807 6 : CALL auxbas_pw_pool%create_pw(pw=wf_r)
808 6 : CALL pw_zero(wf_r)
809 6 : CALL auxbas_pw_pool%create_pw(pw=wf_g)
810 6 : CALL pw_zero(wf_g)
811 : CALL calculate_rho_elec(matrix_p=matrix_p_qs, &
812 : rho=wf_r, &
813 : rho_gspace=wf_g, &
814 6 : ks_env=ks_env)
815 :
816 : ! write this to a cube
817 : CALL cp_pw_to_cube(wf_r, unit_nr=unit_nr, title=title, &
818 6 : particles=particles, stride=stride)
819 :
820 : !free memory
821 6 : CALL auxbas_pw_pool%give_back_pw(wf_r)
822 6 : CALL auxbas_pw_pool%give_back_pw(wf_g)
823 6 : CALL dbcsr_release(matrix_p_qs)
824 :
825 6 : CALL timestop(handle)
826 :
827 6 : END SUBROUTINE write_matrix_to_cube
828 :
829 : ! **************************************************************************************************
830 : !> \brief Initialize g-space density mixing
831 : !> \param qs_env ...
832 : !> \param ls_scf_env ...
833 : ! **************************************************************************************************
834 0 : SUBROUTINE rho_mixing_ls_init(qs_env, ls_scf_env)
835 : TYPE(qs_environment_type), POINTER :: qs_env
836 : TYPE(ls_scf_env_type) :: ls_scf_env
837 :
838 : CHARACTER(len=*), PARAMETER :: routineN = 'rho_mixing_ls_init'
839 :
840 : INTEGER :: handle
841 : TYPE(dft_control_type), POINTER :: dft_control
842 : TYPE(qs_rho_type), POINTER :: rho
843 0 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
844 :
845 0 : CALL timeset(routineN, handle)
846 :
847 0 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
848 :
849 : CALL mixing_allocate(qs_env, ls_scf_env%density_mixing_method, nspins=ls_scf_env%nspins, &
850 0 : mixing_store=ls_scf_env%mixing_store)
851 0 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
852 0 : IF (dft_control%qs_control%gapw) THEN
853 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom)
854 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
855 0 : ls_scf_env%para_env, rho_atom=rho_atom)
856 0 : ELSE IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
857 0 : CALL charge_mixing_init(ls_scf_env%mixing_store)
858 0 : ELSE IF (dft_control%qs_control%semi_empirical) THEN
859 0 : CPABORT('SE Code not possible')
860 : ELSE
861 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
862 0 : ls_scf_env%para_env)
863 : END IF
864 : END IF
865 0 : CALL timestop(handle)
866 0 : END SUBROUTINE rho_mixing_ls_init
867 :
868 : END MODULE dm_ls_scf_qs
|
