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