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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
19 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_get_default_unit_nr,&
22 : cp_logger_type
23 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
24 : USE dbcsr_api, ONLY: &
25 : dbcsr_complete_redistribute, dbcsr_copy, dbcsr_copy_into_existing, dbcsr_create, &
26 : dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_hold, &
27 : dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
28 : dbcsr_finalize, dbcsr_get_info, dbcsr_multiply, dbcsr_nblkrows_total, dbcsr_p_type, &
29 : dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_real_8
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 174 : 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 174 : 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 174 : 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 174 : 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 If set dbcsr_copy_into_existing will be used, 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 : IF (my_keep_sparsity) THEN
331 2232 : CALL dbcsr_copy_into_existing(matrix_qs, matrix_declustered) ! preserve sparsity of matrix_qs
332 : ELSE
333 0 : CALL dbcsr_copy(matrix_qs, matrix_declustered) ! overwrite sparsity of matrix_qs
334 : END IF
335 2232 : CALL dbcsr_release(matrix_declustered)
336 :
337 : ELSE ! using pao
338 1576 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
339 : CALL dbcsr_create(matrix_declustered, &
340 : template=matrix_qs, &
341 : row_blk_size=pao_blk_sizes, &
342 1576 : col_blk_size=pao_blk_sizes)
343 :
344 1576 : CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
345 :
346 1576 : matrix_trafo => ls_mstruct%matrix_B ! contra-variant
347 1576 : IF (covariant) matrix_trafo => ls_mstruct%matrix_A ! co-variant
348 1576 : CALL dbcsr_create(matrix_tmp1, template=matrix_trafo)
349 1576 : CALL dbcsr_create(matrix_tmp2, template=matrix_qs)
350 1576 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_declustered, 0.0_dp, matrix_tmp1)
351 1576 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp1, matrix_trafo, 0.0_dp, matrix_tmp2)
352 1576 : IF (my_keep_sparsity) THEN
353 1298 : CALL dbcsr_copy_into_existing(matrix_qs, matrix_tmp2) ! preserve sparsity of matrix_qs
354 : ELSE
355 278 : CALL dbcsr_copy(matrix_qs, matrix_tmp2) ! overwrite sparsity of matrix_qs
356 : END IF
357 1576 : CALL dbcsr_release(matrix_declustered)
358 1576 : CALL dbcsr_release(matrix_tmp1)
359 1576 : CALL dbcsr_release(matrix_tmp2)
360 : END IF
361 :
362 3808 : CALL timestop(handle)
363 :
364 3808 : END SUBROUTINE matrix_ls_to_qs
365 :
366 : ! **************************************************************************************************
367 : !> \brief Reverses molecular blocking and reduction to single precision if enabled
368 : !> \param matrix_out ...
369 : !> \param matrix_in ...
370 : !> \param ls_mstruct ...
371 : !> \author Ole Schuett
372 : ! **************************************************************************************************
373 11800 : SUBROUTINE matrix_decluster(matrix_out, matrix_in, ls_mstruct)
374 : TYPE(dbcsr_type) :: matrix_out, matrix_in
375 : TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
376 :
377 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_decluster'
378 :
379 : INTEGER :: handle
380 :
381 11800 : CALL timeset(routineN, handle)
382 :
383 22714 : SELECT CASE (ls_mstruct%cluster_type)
384 : CASE (ls_cluster_atomic)
385 10914 : CALL dbcsr_copy(matrix_out, matrix_in) ! takes care of an eventual data_type conversion
386 :
387 : CASE (ls_cluster_molecular)
388 : ! perform the magic complete redistribute copy
389 886 : CALL dbcsr_complete_redistribute(matrix_in, matrix_out)
390 :
391 : CASE DEFAULT
392 11800 : CPABORT("Unknown LS cluster type")
393 : END SELECT
394 :
395 11800 : CALL timestop(handle)
396 :
397 11800 : END SUBROUTINE matrix_decluster
398 :
399 : ! **************************************************************************************************
400 : !> \brief further required initialization of QS.
401 : !> Might be factored-out since this seems common code with the other SCF.
402 : !> \param qs_env ...
403 : !> \par History
404 : !> 2010.10 created [Joost VandeVondele]
405 : !> \author Joost VandeVondele
406 : ! **************************************************************************************************
407 814 : SUBROUTINE ls_scf_init_qs(qs_env)
408 : TYPE(qs_environment_type), POINTER :: qs_env
409 :
410 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_qs'
411 :
412 : INTEGER :: handle, ispin, nspin, unit_nr
413 : TYPE(cp_logger_type), POINTER :: logger
414 814 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
415 : TYPE(dft_control_type), POINTER :: dft_control
416 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
417 814 : POINTER :: sab_orb
418 : TYPE(qs_ks_env_type), POINTER :: ks_env
419 :
420 814 : NULLIFY (sab_orb)
421 814 : CALL timeset(routineN, handle)
422 :
423 : ! get a useful output_unit
424 814 : logger => cp_get_default_logger()
425 814 : IF (logger%para_env%is_source()) THEN
426 407 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
427 : ELSE
428 : unit_nr = -1
429 : END IF
430 :
431 : ! get basic quantities from the qs_env
432 : CALL get_qs_env(qs_env, dft_control=dft_control, &
433 : matrix_s=matrix_s, &
434 : matrix_ks=matrix_ks, &
435 : ks_env=ks_env, &
436 814 : sab_orb=sab_orb)
437 :
438 814 : nspin = dft_control%nspins
439 :
440 : ! we might have to create matrix_ks
441 814 : IF (.NOT. ASSOCIATED(matrix_ks)) THEN
442 0 : CALL dbcsr_allocate_matrix_set(matrix_ks, nspin)
443 0 : DO ispin = 1, nspin
444 0 : ALLOCATE (matrix_ks(ispin)%matrix)
445 0 : CALL dbcsr_create(matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
446 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix, sab_orb)
447 0 : CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
448 : END DO
449 0 : CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
450 : END IF
451 :
452 814 : CALL timestop(handle)
453 :
454 814 : END SUBROUTINE ls_scf_init_qs
455 :
456 : ! **************************************************************************************************
457 : !> \brief get an atomic initial guess
458 : !> \param qs_env ...
459 : !> \param energy ...
460 : !> \par History
461 : !> 2012.11 created [Joost VandeVondele]
462 : !> \author Joost VandeVondele
463 : ! **************************************************************************************************
464 308 : SUBROUTINE ls_scf_qs_atomic_guess(qs_env, energy)
465 : TYPE(qs_environment_type), POINTER :: qs_env
466 : REAL(KIND=dp) :: energy
467 :
468 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_qs_atomic_guess'
469 :
470 : INTEGER :: handle, nspin, unit_nr
471 : INTEGER, DIMENSION(2) :: nelectron_spin
472 : LOGICAL :: has_unit_metric
473 308 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
474 : TYPE(cp_logger_type), POINTER :: logger
475 308 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
476 : TYPE(dft_control_type), POINTER :: dft_control
477 : TYPE(mp_para_env_type), POINTER :: para_env
478 308 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
479 : TYPE(qs_energy_type), POINTER :: qs_energy
480 308 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
481 : TYPE(qs_ks_env_type), POINTER :: ks_env
482 : TYPE(qs_rho_type), POINTER :: rho
483 :
484 308 : CALL timeset(routineN, handle)
485 308 : NULLIFY (rho, rho_ao)
486 :
487 : ! get a useful output_unit
488 308 : logger => cp_get_default_logger()
489 308 : IF (logger%para_env%is_source()) THEN
490 154 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
491 : ELSE
492 154 : unit_nr = -1
493 : END IF
494 :
495 : ! get basic quantities from the qs_env
496 : CALL get_qs_env(qs_env, dft_control=dft_control, &
497 : matrix_s=matrix_s, &
498 : matrix_ks=matrix_ks, &
499 : ks_env=ks_env, &
500 : energy=qs_energy, &
501 : atomic_kind_set=atomic_kind_set, &
502 : qs_kind_set=qs_kind_set, &
503 : particle_set=particle_set, &
504 : has_unit_metric=has_unit_metric, &
505 : para_env=para_env, &
506 : nelectron_spin=nelectron_spin, &
507 308 : rho=rho)
508 :
509 308 : CALL qs_rho_get(rho, rho_ao=rho_ao)
510 :
511 308 : nspin = dft_control%nspins
512 :
513 : ! create an initial atomic guess
514 308 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
515 : dft_control%qs_control%xtb) THEN
516 : CALL calculate_mopac_dm(rho_ao, matrix_s(1)%matrix, has_unit_metric, &
517 : dft_control, particle_set, atomic_kind_set, qs_kind_set, &
518 72 : nspin, nelectron_spin, para_env)
519 : ELSE
520 : CALL calculate_atomic_block_dm(rho_ao, matrix_s(1)%matrix, atomic_kind_set, qs_kind_set, &
521 236 : nspin, nelectron_spin, unit_nr, para_env)
522 : END IF
523 :
524 308 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
525 308 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
526 308 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
527 :
528 308 : energy = qs_energy%total
529 :
530 308 : CALL timestop(handle)
531 :
532 308 : 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 3058 : 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 3058 : 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 3058 : NULLIFY (energy, rho, rho_ao)
561 3058 : CALL timeset(routineN, handle)
562 :
563 3058 : logger => cp_get_default_logger()
564 3058 : IF (logger%para_env%is_source()) THEN
565 1529 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
566 : ELSE
567 : unit_nr = -1
568 : END IF
569 :
570 3058 : nspin = ls_scf_env%nspins
571 3058 : CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
572 3058 : CALL qs_rho_get(rho, rho_ao=rho_ao)
573 :
574 : ! set the new density matrix
575 6342 : DO ispin = 1, nspin
576 : CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
577 6342 : 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 3058 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
582 3058 : IF (ls_scf_env%do_rho_mixing) THEN
583 36 : IF (ls_scf_env%density_mixing_method .EQ. direct_mixing_nr) &
584 0 : CPABORT("Direct P mixing not implemented in linear scaling SCF. ")
585 36 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
586 36 : 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 34 : iscf - 1)
590 34 : IF (unit_nr > 0) THEN
591 : WRITE (unit_nr, '(A57)') &
592 17 : "*********************************************************"
593 : WRITE (unit_nr, '(A13,F5.3,A20,A6,A7,I3)') &
594 17 : " Using ALPHA=", ls_scf_env%mixing_store%alpha, &
595 34 : " 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 17 : " rho_nw=", ls_scf_env%mixing_store%alpha, "*rho + ", &
598 34 : 1.0_dp - ls_scf_env%mixing_store%alpha, "*rho_old"
599 : WRITE (unit_nr, '(A57)') &
600 17 : "*********************************************************"
601 : END IF
602 : END IF
603 : END IF
604 : END IF
605 :
606 3058 : 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 3058 : just_energy=.FALSE., print_active=.TRUE.)
609 3058 : energy_new = energy%total
610 :
611 3058 : CALL timestop(handle)
612 :
613 3058 : END SUBROUTINE ls_scf_dm_to_ks
614 :
615 : ! **************************************************************************************************
616 : !> \brief ...
617 : !> \param qs_env ...
618 : !> \param ls_scf_env ...
619 : !> \param matrix_p_ls ...
620 : !> \param unit_nr ...
621 : !> \param title ...
622 : !> \param stride ...
623 : ! **************************************************************************************************
624 6 : SUBROUTINE write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
625 : TYPE(qs_environment_type), POINTER :: qs_env
626 : TYPE(ls_scf_env_type) :: ls_scf_env
627 : TYPE(dbcsr_type), INTENT(IN) :: matrix_p_ls
628 : INTEGER, INTENT(IN) :: unit_nr
629 : CHARACTER(LEN=*), INTENT(IN) :: title
630 : INTEGER, DIMENSION(:), POINTER :: stride
631 :
632 : CHARACTER(len=*), PARAMETER :: routineN = 'write_matrix_to_cube'
633 :
634 : INTEGER :: handle
635 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
636 : TYPE(dbcsr_type), TARGET :: matrix_p_qs
637 : TYPE(particle_list_type), POINTER :: particles
638 : TYPE(pw_c1d_gs_type) :: wf_g
639 : TYPE(pw_env_type), POINTER :: pw_env
640 6 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
641 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
642 : TYPE(pw_r3d_rs_type) :: wf_r
643 : TYPE(qs_ks_env_type), POINTER :: ks_env
644 : TYPE(qs_subsys_type), POINTER :: subsys
645 :
646 6 : CALL timeset(routineN, handle)
647 :
648 6 : NULLIFY (ks_env, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, matrix_ks)
649 :
650 : CALL get_qs_env(qs_env, &
651 : ks_env=ks_env, &
652 : subsys=subsys, &
653 : pw_env=pw_env, &
654 6 : matrix_ks=matrix_ks)
655 :
656 6 : CALL qs_subsys_get(subsys, particles=particles)
657 :
658 : ! convert the density matrix (ls style) to QS style
659 6 : CALL dbcsr_copy(matrix_p_qs, matrix_ks(1)%matrix)
660 6 : CALL dbcsr_set(matrix_p_qs, 0.0_dp) !zero matrix creation
661 6 : CALL matrix_ls_to_qs(matrix_p_qs, matrix_p_ls, ls_scf_env%ls_mstruct, covariant=.FALSE.)
662 :
663 : ! Print total electronic density
664 : CALL pw_env_get(pw_env=pw_env, &
665 : auxbas_pw_pool=auxbas_pw_pool, &
666 6 : pw_pools=pw_pools)
667 6 : CALL auxbas_pw_pool%create_pw(pw=wf_r)
668 6 : CALL pw_zero(wf_r)
669 6 : CALL auxbas_pw_pool%create_pw(pw=wf_g)
670 6 : CALL pw_zero(wf_g)
671 : CALL calculate_rho_elec(matrix_p=matrix_p_qs, &
672 : rho=wf_r, &
673 : rho_gspace=wf_g, &
674 6 : ks_env=ks_env)
675 :
676 : ! write this to a cube
677 : CALL cp_pw_to_cube(wf_r, unit_nr=unit_nr, title=title, &
678 6 : particles=particles, stride=stride)
679 :
680 : !free memory
681 6 : CALL auxbas_pw_pool%give_back_pw(wf_r)
682 6 : CALL auxbas_pw_pool%give_back_pw(wf_g)
683 6 : CALL dbcsr_release(matrix_p_qs)
684 :
685 6 : CALL timestop(handle)
686 :
687 6 : END SUBROUTINE write_matrix_to_cube
688 :
689 : ! **************************************************************************************************
690 : !> \brief Initialize g-space density mixing
691 : !> \param qs_env ...
692 : !> \param ls_scf_env ...
693 : ! **************************************************************************************************
694 2 : SUBROUTINE rho_mixing_ls_init(qs_env, ls_scf_env)
695 : TYPE(qs_environment_type), POINTER :: qs_env
696 : TYPE(ls_scf_env_type) :: ls_scf_env
697 :
698 : CHARACTER(len=*), PARAMETER :: routineN = 'rho_mixing_ls_init'
699 :
700 : INTEGER :: handle
701 : TYPE(dft_control_type), POINTER :: dft_control
702 : TYPE(qs_rho_type), POINTER :: rho
703 2 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
704 :
705 2 : CALL timeset(routineN, handle)
706 :
707 2 : CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
708 :
709 : CALL mixing_allocate(qs_env, ls_scf_env%density_mixing_method, nspins=ls_scf_env%nspins, &
710 2 : mixing_store=ls_scf_env%mixing_store)
711 2 : IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
712 2 : IF (dft_control%qs_control%gapw) THEN
713 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom)
714 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
715 0 : ls_scf_env%para_env, rho_atom=rho_atom)
716 2 : ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
717 0 : CALL charge_mixing_init(ls_scf_env%mixing_store)
718 2 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
719 0 : CPABORT('SE Code not possible')
720 : ELSE
721 : CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
722 2 : ls_scf_env%para_env)
723 : END IF
724 : END IF
725 2 : CALL timestop(handle)
726 2 : END SUBROUTINE rho_mixing_ls_init
727 :
728 : END MODULE dm_ls_scf_qs
|
