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