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