Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_nblkrows_total, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, &
23 : dbcsr_type_real_8
24 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_get_default_unit_nr,&
28 : cp_logger_type
29 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
30 : USE dm_ls_scf_types, ONLY: ls_cluster_atomic,&
31 : ls_cluster_molecular,&
32 : ls_mstruct_type,&
33 : ls_scf_env_type
34 : USE input_constants, ONLY: 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_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_initial_guess, ONLY: calculate_mopac_dm
57 : USE qs_kind_types, ONLY: qs_kind_type
58 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
59 : USE qs_ks_types, ONLY: qs_ks_did_change,&
60 : qs_ks_env_type,&
61 : set_ks_env
62 : USE qs_mixing_utils, ONLY: charge_mixing_init,&
63 : mixing_allocate,&
64 : mixing_init
65 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
66 : USE qs_rho_atom_types, ONLY: rho_atom_type
67 : USE qs_rho_methods, ONLY: qs_rho_update_rho
68 : USE qs_rho_types, ONLY: qs_rho_get,&
69 : qs_rho_type
70 : USE qs_subsys_types, ONLY: qs_subsys_get,&
71 : qs_subsys_type
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : PRIVATE
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_qs'
79 :
80 : PUBLIC :: matrix_ls_create, matrix_qs_to_ls, matrix_ls_to_qs, ls_scf_init_qs, &
81 : ls_scf_dm_to_ks, ls_scf_qs_atomic_guess, write_matrix_to_cube, rho_mixing_ls_init, &
82 : matrix_decluster
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief create a matrix for use (and as a template) in ls based on a qs template
88 : !> \param matrix_ls ...
89 : !> \param matrix_qs ...
90 : !> \param ls_mstruct ...
91 : !> \par History
92 : !> 2011.03 created [Joost VandeVondele]
93 : !> 2015.09 add support for PAO [Ole Schuett]
94 : !> \author Joost VandeVondele
95 : ! **************************************************************************************************
96 844 : SUBROUTINE matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
97 : TYPE(dbcsr_type) :: matrix_ls, matrix_qs
98 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
99 :
100 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_ls_create'
101 :
102 : CHARACTER(len=default_string_length) :: name
103 : INTEGER :: handle, iatom, imol, jatom, &
104 : ls_data_type, natom, nmol
105 844 : INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: atom_to_cluster, atom_to_cluster_primus, &
106 844 : clustered_blk_sizes, primus_of_mol
107 844 : INTEGER, DIMENSION(:), POINTER :: clustered_col_dist, clustered_row_dist, &
108 844 : ls_blk_sizes, ls_col_dist, ls_row_dist
109 : TYPE(dbcsr_distribution_type) :: ls_dist, ls_dist_clustered
110 :
111 844 : CALL timeset(routineN, handle)
112 :
113 : ! Defaults -----------------------------------------------------------------------------------
114 844 : CALL dbcsr_get_info(matrix_qs, col_blk_size=ls_blk_sizes, distribution=ls_dist)
115 844 : CALL dbcsr_distribution_hold(ls_dist)
116 844 : CALL dbcsr_distribution_get(ls_dist, row_dist=ls_row_dist, col_dist=ls_col_dist)
117 844 : ls_data_type = dbcsr_type_real_8
118 :
119 : ! PAO ----------------------------------------------------------------------------------------
120 844 : IF (ls_mstruct%do_pao) THEN
121 496 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=ls_blk_sizes)
122 : END IF
123 :
124 : ! Clustering ---------------------------------------------------------------------------------
125 902 : 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 : natom = dbcsr_nblkrows_total(matrix_qs)
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 58 : ls_dist = ls_dist_clustered
178 :
179 : CASE DEFAULT
180 844 : CPABORT("Unknown LS cluster type")
181 : END SELECT
182 :
183 : ! Create actual matrix -----------------------------------------------------------------------
184 844 : 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 : data_type=ls_data_type, &
190 : row_blk_size=ls_blk_sizes, &
191 844 : col_blk_size=ls_blk_sizes)
192 844 : CALL dbcsr_distribution_release(ls_dist)
193 844 : CALL dbcsr_finalize(matrix_ls)
194 :
195 844 : CALL timestop(handle)
196 :
197 1688 : END SUBROUTINE matrix_ls_create
198 :
199 : ! **************************************************************************************************
200 : !> \brief first link to QS, copy a QS matrix to LS matrix
201 : !> used to isolate QS style matrices from LS style
202 : !> will be useful for future features (e.g. precision, symmetry, blocking, ...)
203 : !> \param matrix_ls ...
204 : !> \param matrix_qs ...
205 : !> \param ls_mstruct ...
206 : !> \param covariant ...
207 : !> \par History
208 : !> 2010.10 created [Joost VandeVondele]
209 : !> 2015.09 add support for PAO [Ole Schuett]
210 : !> \author Joost VandeVondele
211 : ! **************************************************************************************************
212 27698 : SUBROUTINE matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
213 : TYPE(dbcsr_type) :: matrix_ls, matrix_qs
214 : TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
215 : LOGICAL, INTENT(IN) :: covariant
216 :
217 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_qs_to_ls'
218 :
219 : INTEGER :: handle
220 27698 : INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
221 : TYPE(dbcsr_type) :: matrix_pao, matrix_tmp
222 : TYPE(dbcsr_type), POINTER :: matrix_trafo
223 :
224 27698 : CALL timeset(routineN, handle)
225 :
226 27698 : IF (.NOT. ls_mstruct%do_pao) THEN
227 2328 : CALL matrix_cluster(matrix_ls, matrix_qs, ls_mstruct)
228 :
229 : ELSE ! using pao
230 25370 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
231 : CALL dbcsr_create(matrix_pao, &
232 : matrix_type="N", &
233 : template=matrix_qs, &
234 : row_blk_size=pao_blk_sizes, &
235 25370 : col_blk_size=pao_blk_sizes)
236 :
237 25370 : matrix_trafo => ls_mstruct%matrix_A ! contra-variant
238 25370 : IF (covariant) matrix_trafo => ls_mstruct%matrix_B ! co-variant
239 25370 : CALL dbcsr_create(matrix_tmp, template=matrix_trafo)
240 :
241 25370 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qs, matrix_trafo, 0.0_dp, matrix_tmp)
242 25370 : CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_trafo, matrix_tmp, 0.0_dp, matrix_pao)
243 25370 : CALL dbcsr_release(matrix_tmp)
244 :
245 25370 : CALL matrix_cluster(matrix_ls, matrix_pao, ls_mstruct)
246 25370 : CALL dbcsr_release(matrix_pao)
247 : END IF
248 :
249 27698 : CALL timestop(handle)
250 :
251 27698 : END SUBROUTINE matrix_qs_to_ls
252 :
253 : ! **************************************************************************************************
254 : !> \brief Performs molecular blocking and reduction to single precision if enabled
255 : !> \param matrix_out ...
256 : !> \param matrix_in ...
257 : !> \param ls_mstruct ...
258 : !> \author Ole Schuett
259 : ! **************************************************************************************************
260 27698 : SUBROUTINE matrix_cluster(matrix_out, matrix_in, ls_mstruct)
261 : TYPE(dbcsr_type) :: matrix_out, matrix_in
262 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
263 :
264 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_cluster'
265 :
266 : INTEGER :: handle
267 : TYPE(dbcsr_type) :: matrix_in_nosym
268 :
269 27698 : CALL timeset(routineN, handle)
270 :
271 53836 : SELECT CASE (ls_mstruct%cluster_type)
272 : CASE (ls_cluster_atomic)
273 26138 : CALL dbcsr_copy(matrix_out, matrix_in) ! takes care of an eventual data_type conversion
274 :
275 : CASE (ls_cluster_molecular)
276 : ! desymmetrize the qs matrix
277 1560 : CALL dbcsr_create(matrix_in_nosym, template=matrix_in, matrix_type="N")
278 1560 : CALL dbcsr_desymmetrize(matrix_in, matrix_in_nosym)
279 :
280 : ! perform the magic complete redistribute copy
281 1560 : CALL dbcsr_complete_redistribute(matrix_in_nosym, matrix_out);
282 1560 : CALL dbcsr_release(matrix_in_nosym)
283 :
284 : CASE DEFAULT
285 27698 : CPABORT("Unknown LS cluster type")
286 : END SELECT
287 :
288 27698 : CALL timestop(handle)
289 :
290 27698 : END SUBROUTINE matrix_cluster
291 :
292 : ! **************************************************************************************************
293 : !> \brief second link to QS, copy a LS matrix to QS matrix
294 : !> used to isolate QS style matrices from LS style
295 : !> will be useful for future features (e.g. precision, symmetry, blocking, ...)
296 : !> \param matrix_qs ...
297 : !> \param matrix_ls ...
298 : !> \param ls_mstruct ...
299 : !> \param covariant ...
300 : !> \param keep_sparsity will be passed on to dbcsr_copy, by default set to .TRUE.
301 : !> \par History
302 : !> 2010.10 created [Joost VandeVondele]
303 : !> 2015.09 add support for PAO [Ole Schuett]
304 : !> \author Joost VandeVondele
305 : ! **************************************************************************************************
306 3808 : SUBROUTINE matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
307 : TYPE(dbcsr_type) :: matrix_qs, matrix_ls
308 : TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
309 : LOGICAL :: covariant
310 : LOGICAL, OPTIONAL :: keep_sparsity
311 :
312 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_ls_to_qs'
313 :
314 : INTEGER :: handle
315 3808 : INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
316 : LOGICAL :: my_keep_sparsity
317 : TYPE(dbcsr_type) :: matrix_declustered, matrix_tmp1, &
318 : matrix_tmp2
319 : TYPE(dbcsr_type), POINTER :: matrix_trafo
320 :
321 3808 : CALL timeset(routineN, handle)
322 :
323 3808 : my_keep_sparsity = .TRUE.
324 3808 : IF (PRESENT(keep_sparsity)) &
325 278 : my_keep_sparsity = keep_sparsity
326 :
327 3808 : IF (.NOT. ls_mstruct%do_pao) THEN
328 2232 : CALL dbcsr_create(matrix_declustered, template=matrix_qs)
329 2232 : CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
330 2232 : CALL dbcsr_copy(matrix_qs, matrix_declustered, keep_sparsity=my_keep_sparsity)
331 2232 : CALL dbcsr_release(matrix_declustered)
332 :
333 : ELSE ! using pao
334 1576 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
335 : CALL dbcsr_create(matrix_declustered, &
336 : template=matrix_qs, &
337 : row_blk_size=pao_blk_sizes, &
338 1576 : col_blk_size=pao_blk_sizes)
339 :
340 1576 : CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
341 :
342 1576 : matrix_trafo => ls_mstruct%matrix_B ! contra-variant
343 1576 : IF (covariant) matrix_trafo => ls_mstruct%matrix_A ! co-variant
344 1576 : CALL dbcsr_create(matrix_tmp1, template=matrix_trafo)
345 1576 : CALL dbcsr_create(matrix_tmp2, template=matrix_qs)
346 1576 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_declustered, 0.0_dp, matrix_tmp1)
347 1576 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp1, matrix_trafo, 0.0_dp, matrix_tmp2)
348 1576 : CALL dbcsr_copy(matrix_qs, matrix_tmp2, keep_sparsity=my_keep_sparsity)
349 1576 : CALL dbcsr_release(matrix_declustered)
350 1576 : CALL dbcsr_release(matrix_tmp1)
351 1576 : CALL dbcsr_release(matrix_tmp2)
352 : END IF
353 :
354 3808 : CALL timestop(handle)
355 :
356 3808 : END SUBROUTINE matrix_ls_to_qs
357 :
358 : ! **************************************************************************************************
359 : !> \brief Reverses molecular blocking and reduction to single precision if enabled
360 : !> \param matrix_out ...
361 : !> \param matrix_in ...
362 : !> \param ls_mstruct ...
363 : !> \author Ole Schuett
364 : ! **************************************************************************************************
365 11794 : SUBROUTINE matrix_decluster(matrix_out, matrix_in, ls_mstruct)
366 : TYPE(dbcsr_type) :: matrix_out, matrix_in
367 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
368 :
369 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_decluster'
370 :
371 : INTEGER :: handle
372 :
373 11794 : CALL timeset(routineN, handle)
374 :
375 22702 : SELECT CASE (ls_mstruct%cluster_type)
376 : CASE (ls_cluster_atomic)
377 10908 : CALL dbcsr_copy(matrix_out, matrix_in) ! takes care of an eventual data_type conversion
378 :
379 : CASE (ls_cluster_molecular)
380 : ! perform the magic complete redistribute copy
381 886 : CALL dbcsr_complete_redistribute(matrix_in, matrix_out)
382 :
383 : CASE DEFAULT
384 11794 : CPABORT("Unknown LS cluster type")
385 : END SELECT
386 :
387 11794 : CALL timestop(handle)
388 :
389 11794 : END SUBROUTINE matrix_decluster
390 :
391 : ! **************************************************************************************************
392 : !> \brief further required initialization of QS.
393 : !> Might be factored-out since this seems common code with the other SCF.
394 : !> \param qs_env ...
395 : !> \par History
396 : !> 2010.10 created [Joost VandeVondele]
397 : !> \author Joost VandeVondele
398 : ! **************************************************************************************************
399 814 : SUBROUTINE ls_scf_init_qs(qs_env)
400 : TYPE(qs_environment_type), POINTER :: qs_env
401 :
402 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_qs'
403 :
404 : INTEGER :: handle, ispin, nspin, unit_nr
405 : TYPE(cp_logger_type), POINTER :: logger
406 814 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
407 : TYPE(dft_control_type), POINTER :: dft_control
408 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
409 814 : POINTER :: sab_orb
410 : TYPE(qs_ks_env_type), POINTER :: ks_env
411 :
412 814 : NULLIFY (sab_orb)
413 814 : CALL timeset(routineN, handle)
414 :
415 : ! get a useful output_unit
416 814 : logger => cp_get_default_logger()
417 814 : IF (logger%para_env%is_source()) THEN
418 407 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
419 : ELSE
420 : unit_nr = -1
421 : END IF
422 :
423 : ! get basic quantities from the qs_env
424 : CALL get_qs_env(qs_env, dft_control=dft_control, &
425 : matrix_s=matrix_s, &
426 : matrix_ks=matrix_ks, &
427 : ks_env=ks_env, &
428 814 : sab_orb=sab_orb)
429 :
430 814 : nspin = dft_control%nspins
431 :
432 : ! we might have to create matrix_ks
433 814 : IF (.NOT. ASSOCIATED(matrix_ks)) THEN
434 0 : CALL dbcsr_allocate_matrix_set(matrix_ks, nspin)
435 0 : DO ispin = 1, nspin
436 0 : ALLOCATE (matrix_ks(ispin)%matrix)
437 0 : CALL dbcsr_create(matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
438 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix, sab_orb)
439 0 : CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
440 : END DO
441 0 : CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
442 : END IF
443 :
444 814 : CALL timestop(handle)
445 :
446 814 : END SUBROUTINE ls_scf_init_qs
447 :
448 : ! **************************************************************************************************
449 : !> \brief get an atomic initial guess
450 : !> \param qs_env ...
451 : !> \param energy ...
452 : !> \par History
453 : !> 2012.11 created [Joost VandeVondele]
454 : !> \author Joost VandeVondele
455 : ! **************************************************************************************************
456 308 : SUBROUTINE ls_scf_qs_atomic_guess(qs_env, energy)
457 : TYPE(qs_environment_type), POINTER :: qs_env
458 : REAL(KIND=dp) :: energy
459 :
460 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_qs_atomic_guess'
461 :
462 : INTEGER :: handle, nspin, unit_nr
463 : INTEGER, DIMENSION(2) :: nelectron_spin
464 : LOGICAL :: has_unit_metric
465 308 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
466 : TYPE(cp_logger_type), POINTER :: logger
467 308 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
468 : TYPE(dft_control_type), POINTER :: dft_control
469 : TYPE(mp_para_env_type), POINTER :: para_env
470 308 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
471 : TYPE(qs_energy_type), POINTER :: qs_energy
472 308 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
473 : TYPE(qs_ks_env_type), POINTER :: ks_env
474 : TYPE(qs_rho_type), POINTER :: rho
475 :
476 308 : CALL timeset(routineN, handle)
477 308 : NULLIFY (rho, rho_ao)
478 :
479 : ! get a useful output_unit
480 308 : logger => cp_get_default_logger()
481 308 : IF (logger%para_env%is_source()) THEN
482 154 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
483 : ELSE
484 154 : unit_nr = -1
485 : END IF
486 :
487 : ! get basic quantities from the qs_env
488 : CALL get_qs_env(qs_env, dft_control=dft_control, &
489 : matrix_s=matrix_s, &
490 : matrix_ks=matrix_ks, &
491 : ks_env=ks_env, &
492 : energy=qs_energy, &
493 : atomic_kind_set=atomic_kind_set, &
494 : qs_kind_set=qs_kind_set, &
495 : particle_set=particle_set, &
496 : has_unit_metric=has_unit_metric, &
497 : para_env=para_env, &
498 : nelectron_spin=nelectron_spin, &
499 308 : rho=rho)
500 :
501 308 : CALL qs_rho_get(rho, rho_ao=rho_ao)
502 :
503 308 : nspin = dft_control%nspins
504 :
505 : ! create an initial atomic guess
506 308 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
507 : dft_control%qs_control%xtb) THEN
508 : CALL calculate_mopac_dm(rho_ao, matrix_s(1)%matrix, has_unit_metric, &
509 : dft_control, particle_set, atomic_kind_set, qs_kind_set, &
510 72 : nspin, nelectron_spin, para_env)
511 : ELSE
512 : CALL calculate_atomic_block_dm(rho_ao, matrix_s(1)%matrix, atomic_kind_set, qs_kind_set, &
513 236 : nspin, nelectron_spin, unit_nr, para_env)
514 : END IF
515 :
516 308 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
517 308 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
518 308 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
519 :
520 308 : energy = qs_energy%total
521 :
522 308 : CALL timestop(handle)
523 :
524 308 : END SUBROUTINE ls_scf_qs_atomic_guess
525 :
526 : ! **************************************************************************************************
527 : !> \brief use the density matrix in ls_scf_env to compute the new energy and KS matrix
528 : !> \param qs_env ...
529 : !> \param ls_scf_env ...
530 : !> \param energy_new ...
531 : !> \param iscf ...
532 : !> \par History
533 : !> 2011.04 created [Joost VandeVondele]
534 : !> 2015.02 added gspace density mixing [Patrick Seewald]
535 : !> \author Joost VandeVondele
536 : ! **************************************************************************************************
537 3058 : SUBROUTINE ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
538 : TYPE(qs_environment_type), POINTER :: qs_env
539 : TYPE(ls_scf_env_type) :: ls_scf_env
540 : REAL(KIND=dp) :: energy_new
541 : INTEGER, INTENT(IN) :: iscf
542 :
543 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_dm_to_ks'
544 :
545 : INTEGER :: handle, ispin, nspin, unit_nr
546 : TYPE(cp_logger_type), POINTER :: logger
547 3058 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
548 : TYPE(mp_para_env_type), POINTER :: para_env
549 : TYPE(qs_energy_type), POINTER :: energy
550 : TYPE(qs_rho_type), POINTER :: rho
551 :
552 3058 : NULLIFY (energy, rho, rho_ao)
553 3058 : CALL timeset(routineN, handle)
554 :
555 3058 : logger => cp_get_default_logger()
556 3058 : IF (logger%para_env%is_source()) THEN
557 1529 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
558 : ELSE
559 : unit_nr = -1
560 : END IF
561 :
562 3058 : nspin = ls_scf_env%nspins
563 3058 : CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
564 3058 : CALL qs_rho_get(rho, rho_ao=rho_ao)
565 :
566 : ! set the new density matrix
567 6342 : DO ispin = 1, nspin
568 : CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
569 6342 : ls_scf_env%ls_mstruct, covariant=.FALSE.)
570 : END DO
571 :
572 : ! compute the corresponding KS matrix and new energy, mix density if requested
573 3058 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
574 3058 : IF (ls_scf_env%do_rho_mixing) THEN
575 36 : IF (ls_scf_env%density_mixing_method .EQ. direct_mixing_nr) &
576 0 : CPABORT("Direct P mixing not implemented in linear scaling SCF. ")
577 36 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
578 36 : IF (iscf .GT. MAX(ls_scf_env%mixing_store%nskip_mixing, 1)) THEN
579 : CALL gspace_mixing(qs_env, ls_scf_env%density_mixing_method, &
580 : ls_scf_env%mixing_store, rho, para_env, &
581 34 : iscf - 1)
582 34 : IF (unit_nr > 0) THEN
583 : WRITE (unit_nr, '(A57)') &
584 17 : "*********************************************************"
585 : WRITE (unit_nr, '(A13,F5.3,A20,A6,A7,I3)') &
586 17 : " Using ALPHA=", ls_scf_env%mixing_store%alpha, &
587 34 : " to mix rho: method=", ls_scf_env%mixing_store%iter_method, ", iscf=", iscf
588 : WRITE (unit_nr, '(A8,F5.3,A6,F5.3,A8)') &
589 17 : " rho_nw=", ls_scf_env%mixing_store%alpha, "*rho + ", &
590 34 : 1.0_dp - ls_scf_env%mixing_store%alpha, "*rho_old"
591 : WRITE (unit_nr, '(A57)') &
592 17 : "*********************************************************"
593 : END IF
594 : END IF
595 : END IF
596 : END IF
597 :
598 3058 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
599 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
600 3058 : just_energy=.FALSE., print_active=.TRUE.)
601 3058 : energy_new = energy%total
602 :
603 3058 : CALL timestop(handle)
604 :
605 3058 : END SUBROUTINE ls_scf_dm_to_ks
606 :
607 : ! **************************************************************************************************
608 : !> \brief ...
609 : !> \param qs_env ...
610 : !> \param ls_scf_env ...
611 : !> \param matrix_p_ls ...
612 : !> \param unit_nr ...
613 : !> \param title ...
614 : !> \param stride ...
615 : ! **************************************************************************************************
616 6 : SUBROUTINE write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
617 : TYPE(qs_environment_type), POINTER :: qs_env
618 : TYPE(ls_scf_env_type) :: ls_scf_env
619 : TYPE(dbcsr_type), INTENT(IN) :: matrix_p_ls
620 : INTEGER, INTENT(IN) :: unit_nr
621 : CHARACTER(LEN=*), INTENT(IN) :: title
622 : INTEGER, DIMENSION(:), POINTER :: stride
623 :
624 : CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_to_cube'
625 :
626 : INTEGER :: handle
627 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
628 : TYPE(dbcsr_type), TARGET :: matrix_p_qs
629 : TYPE(particle_list_type), POINTER :: particles
630 : TYPE(pw_c1d_gs_type) :: wf_g
631 : TYPE(pw_env_type), POINTER :: pw_env
632 6 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
633 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
634 : TYPE(pw_r3d_rs_type) :: wf_r
635 : TYPE(qs_ks_env_type), POINTER :: ks_env
636 : TYPE(qs_subsys_type), POINTER :: subsys
637 :
638 6 : CALL timeset(routineN, handle)
639 :
640 6 : NULLIFY (ks_env, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, matrix_ks)
641 :
642 : CALL get_qs_env(qs_env, &
643 : ks_env=ks_env, &
644 : subsys=subsys, &
645 : pw_env=pw_env, &
646 6 : matrix_ks=matrix_ks)
647 :
648 6 : CALL qs_subsys_get(subsys, particles=particles)
649 :
650 : ! convert the density matrix (ls style) to QS style
651 6 : CALL dbcsr_copy(matrix_p_qs, matrix_ks(1)%matrix)
652 6 : CALL dbcsr_set(matrix_p_qs, 0.0_dp) !zero matrix creation
653 6 : CALL matrix_ls_to_qs(matrix_p_qs, matrix_p_ls, ls_scf_env%ls_mstruct, covariant=.FALSE.)
654 :
655 : ! Print total electronic density
656 : CALL pw_env_get(pw_env=pw_env, &
657 : auxbas_pw_pool=auxbas_pw_pool, &
658 6 : pw_pools=pw_pools)
659 6 : CALL auxbas_pw_pool%create_pw(pw=wf_r)
660 6 : CALL pw_zero(wf_r)
661 6 : CALL auxbas_pw_pool%create_pw(pw=wf_g)
662 6 : CALL pw_zero(wf_g)
663 : CALL calculate_rho_elec(matrix_p=matrix_p_qs, &
664 : rho=wf_r, &
665 : rho_gspace=wf_g, &
666 6 : ks_env=ks_env)
667 :
668 : ! write this to a cube
669 : CALL cp_pw_to_cube(wf_r, unit_nr=unit_nr, title=title, &
670 6 : particles=particles, stride=stride)
671 :
672 : !free memory
673 6 : CALL auxbas_pw_pool%give_back_pw(wf_r)
674 6 : CALL auxbas_pw_pool%give_back_pw(wf_g)
675 6 : CALL dbcsr_release(matrix_p_qs)
676 :
677 6 : CALL timestop(handle)
678 :
679 6 : END SUBROUTINE write_matrix_to_cube
680 :
681 : ! **************************************************************************************************
682 : !> \brief Initialize g-space density mixing
683 : !> \param qs_env ...
684 : !> \param ls_scf_env ...
685 : ! **************************************************************************************************
686 2 : SUBROUTINE rho_mixing_ls_init(qs_env, ls_scf_env)
687 : TYPE(qs_environment_type), POINTER :: qs_env
688 : TYPE(ls_scf_env_type) :: ls_scf_env
689 :
690 : CHARACTER(len=*), PARAMETER :: routineN = 'rho_mixing_ls_init'
691 :
692 : INTEGER :: handle
693 : TYPE(dft_control_type), POINTER :: dft_control
694 : TYPE(qs_rho_type), POINTER :: rho
695 2 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
696 :
697 2 : CALL timeset(routineN, handle)
698 :
699 2 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
700 :
701 : CALL mixing_allocate(qs_env, ls_scf_env%density_mixing_method, nspins=ls_scf_env%nspins, &
702 2 : mixing_store=ls_scf_env%mixing_store)
703 2 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
704 2 : IF (dft_control%qs_control%gapw) THEN
705 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom)
706 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
707 0 : ls_scf_env%para_env, rho_atom=rho_atom)
708 2 : ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
709 0 : CALL charge_mixing_init(ls_scf_env%mixing_store)
710 2 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
711 0 : CPABORT('SE Code not possible')
712 : ELSE
713 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
714 2 : ls_scf_env%para_env)
715 : END IF
716 : END IF
717 2 : CALL timestop(handle)
718 2 : END SUBROUTINE rho_mixing_ls_init
719 :
720 : END MODULE dm_ls_scf_qs
|
