(git:ccc2433)
dm_ls_scf_qs.F
Go to the documentation of this file.
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 ! **************************************************************************************************
16  USE atomic_kind_types, ONLY: atomic_kind_type
17  USE cp_control_types, ONLY: dft_control_type
22  cp_logger_type
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
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
52  USE qs_energy_types, ONLY: qs_energy_type
53  USE qs_environment_types, ONLY: get_qs_env,&
54  qs_environment_type
57  USE qs_kind_types, ONLY: qs_kind_type
59  USE qs_ks_types, ONLY: qs_ks_did_change,&
60  qs_ks_env_type,&
65  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
66  USE qs_rho_atom_types, ONLY: rho_atom_type
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 
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  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  INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: atom_to_cluster, atom_to_cluster_primus, &
106  clustered_blk_sizes, primus_of_mol
107  INTEGER, DIMENSION(:), POINTER :: clustered_col_dist, clustered_row_dist, &
108  ls_blk_sizes, ls_col_dist, ls_row_dist
109  TYPE(dbcsr_distribution_type) :: ls_dist, ls_dist_clustered
110 
111  CALL timeset(routinen, handle)
112 
113  ! Defaults -----------------------------------------------------------------------------------
114  CALL dbcsr_get_info(matrix_qs, col_blk_size=ls_blk_sizes, distribution=ls_dist)
115  CALL dbcsr_distribution_hold(ls_dist)
116  CALL dbcsr_distribution_get(ls_dist, row_dist=ls_row_dist, col_dist=ls_col_dist)
117  ls_data_type = dbcsr_type_real_8
118 
119  ! PAO ----------------------------------------------------------------------------------------
120  IF (ls_mstruct%do_pao) THEN
121  CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=ls_blk_sizes)
122  END IF
123 
124  ! Clustering ---------------------------------------------------------------------------------
125  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  natom = dbcsr_nblkrows_total(matrix_qs)
131  nmol = maxval(ls_mstruct%atom_to_molecule)
132  ALLOCATE (atom_to_cluster_primus(natom))
133  ALLOCATE (atom_to_cluster(natom))
134  ALLOCATE (primus_of_mol(nmol))
135  DO iatom = 1, natom
136  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  DO jatom = iatom, 1, -1
141  IF (ls_mstruct%atom_to_molecule(jatom) == atom_to_cluster(iatom)) THEN
142  atom_to_cluster_primus(iatom) = jatom
143  ELSE
144  EXIT
145  END IF
146  END DO
147  primus_of_mol(atom_to_cluster(iatom)) = atom_to_cluster_primus(iatom)
148  END DO
149 
150  ! row
151  ALLOCATE (clustered_row_dist(nmol))
152  DO imol = 1, nmol
153  clustered_row_dist(imol) = ls_row_dist(primus_of_mol(imol))
154  END DO
155 
156  ! col
157  ALLOCATE (clustered_col_dist(nmol))
158  DO imol = 1, nmol
159  clustered_col_dist(imol) = ls_col_dist(primus_of_mol(imol))
160  END DO
161 
162  ALLOCATE (clustered_blk_sizes(nmol))
163  clustered_blk_sizes = 0
164  DO iatom = 1, natom
165  clustered_blk_sizes(atom_to_cluster(iatom)) = clustered_blk_sizes(atom_to_cluster(iatom)) + &
166  ls_blk_sizes(iatom)
167  END DO
168  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  reuse_arrays=.true.)
176  CALL dbcsr_distribution_release(ls_dist)
177  ls_dist = ls_dist_clustered
178 
179  CASE DEFAULT
180  cpabort("Unknown LS cluster type")
181  END SELECT
182 
183  ! Create actual matrix -----------------------------------------------------------------------
184  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  col_blk_size=ls_blk_sizes)
192  CALL dbcsr_distribution_release(ls_dist)
193  CALL dbcsr_finalize(matrix_ls)
194 
195  CALL timestop(handle)
196 
197  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  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  INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
221  TYPE(dbcsr_type) :: matrix_pao, matrix_tmp
222  TYPE(dbcsr_type), POINTER :: matrix_trafo
223 
224  CALL timeset(routinen, handle)
225 
226  IF (.NOT. ls_mstruct%do_pao) THEN
227  CALL matrix_cluster(matrix_ls, matrix_qs, ls_mstruct)
228 
229  ELSE ! using pao
230  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  col_blk_size=pao_blk_sizes)
236 
237  matrix_trafo => ls_mstruct%matrix_A ! contra-variant
238  IF (covariant) matrix_trafo => ls_mstruct%matrix_B ! co-variant
239  CALL dbcsr_create(matrix_tmp, template=matrix_trafo)
240 
241  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qs, matrix_trafo, 0.0_dp, matrix_tmp)
242  CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_trafo, matrix_tmp, 0.0_dp, matrix_pao)
243  CALL dbcsr_release(matrix_tmp)
244 
245  CALL matrix_cluster(matrix_ls, matrix_pao, ls_mstruct)
246  CALL dbcsr_release(matrix_pao)
247  END IF
248 
249  CALL timestop(handle)
250 
251  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  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  CALL timeset(routinen, handle)
270 
271  SELECT CASE (ls_mstruct%cluster_type)
272  CASE (ls_cluster_atomic)
273  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  CALL dbcsr_create(matrix_in_nosym, template=matrix_in, matrix_type="N")
278  CALL dbcsr_desymmetrize(matrix_in, matrix_in_nosym)
279 
280  ! perform the magic complete redistribute copy
281  CALL dbcsr_complete_redistribute(matrix_in_nosym, matrix_out);
282  CALL dbcsr_release(matrix_in_nosym)
283 
284  CASE DEFAULT
285  cpabort("Unknown LS cluster type")
286  END SELECT
287 
288  CALL timestop(handle)
289 
290  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  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  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  CALL timeset(routinen, handle)
322 
323  my_keep_sparsity = .true.
324  IF (PRESENT(keep_sparsity)) &
325  my_keep_sparsity = keep_sparsity
326 
327  IF (.NOT. ls_mstruct%do_pao) THEN
328  CALL dbcsr_create(matrix_declustered, template=matrix_qs)
329  CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
330  IF (my_keep_sparsity) THEN
331  CALL dbcsr_copy_into_existing(matrix_qs, matrix_declustered) ! preserve sparsity of matrix_qs
332  ELSE
333  CALL dbcsr_copy(matrix_qs, matrix_declustered) ! overwrite sparsity of matrix_qs
334  END IF
335  CALL dbcsr_release(matrix_declustered)
336 
337  ELSE ! using pao
338  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  col_blk_size=pao_blk_sizes)
343 
344  CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
345 
346  matrix_trafo => ls_mstruct%matrix_B ! contra-variant
347  IF (covariant) matrix_trafo => ls_mstruct%matrix_A ! co-variant
348  CALL dbcsr_create(matrix_tmp1, template=matrix_trafo)
349  CALL dbcsr_create(matrix_tmp2, template=matrix_qs)
350  CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_declustered, 0.0_dp, matrix_tmp1)
351  CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp1, matrix_trafo, 0.0_dp, matrix_tmp2)
352  IF (my_keep_sparsity) THEN
353  CALL dbcsr_copy_into_existing(matrix_qs, matrix_tmp2) ! preserve sparsity of matrix_qs
354  ELSE
355  CALL dbcsr_copy(matrix_qs, matrix_tmp2) ! overwrite sparsity of matrix_qs
356  END IF
357  CALL dbcsr_release(matrix_declustered)
358  CALL dbcsr_release(matrix_tmp1)
359  CALL dbcsr_release(matrix_tmp2)
360  END IF
361 
362  CALL timestop(handle)
363 
364  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  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  CALL timeset(routinen, handle)
382 
383  SELECT CASE (ls_mstruct%cluster_type)
384  CASE (ls_cluster_atomic)
385  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  CALL dbcsr_complete_redistribute(matrix_in, matrix_out)
390 
391  CASE DEFAULT
392  cpabort("Unknown LS cluster type")
393  END SELECT
394 
395  CALL timestop(handle)
396 
397  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  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  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  POINTER :: sab_orb
418  TYPE(qs_ks_env_type), POINTER :: ks_env
419 
420  NULLIFY (sab_orb)
421  CALL timeset(routinen, handle)
422 
423  ! get a useful output_unit
424  logger => cp_get_default_logger()
425  IF (logger%para_env%is_source()) THEN
426  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  sab_orb=sab_orb)
437 
438  nspin = dft_control%nspins
439 
440  ! we might have to create matrix_ks
441  IF (.NOT. ASSOCIATED(matrix_ks)) THEN
442  CALL dbcsr_allocate_matrix_set(matrix_ks, nspin)
443  DO ispin = 1, nspin
444  ALLOCATE (matrix_ks(ispin)%matrix)
445  CALL dbcsr_create(matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
446  CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix, sab_orb)
447  CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
448  END DO
449  CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
450  END IF
451 
452  CALL timestop(handle)
453 
454  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  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  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
474  TYPE(cp_logger_type), POINTER :: logger
475  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
479  TYPE(qs_energy_type), POINTER :: qs_energy
480  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  CALL timeset(routinen, handle)
485  NULLIFY (rho, rho_ao)
486 
487  ! get a useful output_unit
488  logger => cp_get_default_logger()
489  IF (logger%para_env%is_source()) THEN
490  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
491  ELSE
492  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  rho=rho)
508 
509  CALL qs_rho_get(rho, rho_ao=rho_ao)
510 
511  nspin = dft_control%nspins
512 
513  ! create an initial atomic guess
514  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  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  nspin, nelectron_spin, unit_nr, para_env)
522  END IF
523 
524  CALL qs_rho_update_rho(rho, qs_env=qs_env)
525  CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
526  CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., just_energy=.false.)
527 
528  energy = qs_energy%total
529 
530  CALL timestop(handle)
531 
532  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  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  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  NULLIFY (energy, rho, rho_ao)
561  CALL timeset(routinen, handle)
562 
563  logger => cp_get_default_logger()
564  IF (logger%para_env%is_source()) THEN
565  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
566  ELSE
567  unit_nr = -1
568  END IF
569 
570  nspin = ls_scf_env%nspins
571  CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
572  CALL qs_rho_get(rho, rho_ao=rho_ao)
573 
574  ! set the new density matrix
575  DO ispin = 1, nspin
576  CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
577  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  CALL qs_rho_update_rho(rho, qs_env=qs_env)
582  IF (ls_scf_env%do_rho_mixing) THEN
583  IF (ls_scf_env%density_mixing_method .EQ. direct_mixing_nr) &
584  cpabort("Direct P mixing not implemented in linear scaling SCF. ")
585  IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
586  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  iscf - 1)
590  IF (unit_nr > 0) THEN
591  WRITE (unit_nr, '(A57)') &
592  "*********************************************************"
593  WRITE (unit_nr, '(A13,F5.3,A20,A6,A7,I3)') &
594  " Using ALPHA=", ls_scf_env%mixing_store%alpha, &
595  " 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  " rho_nw=", ls_scf_env%mixing_store%alpha, "*rho + ", &
598  1.0_dp - ls_scf_env%mixing_store%alpha, "*rho_old"
599  WRITE (unit_nr, '(A57)') &
600  "*********************************************************"
601  END IF
602  END IF
603  END IF
604  END IF
605 
606  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  just_energy=.false., print_active=.true.)
609  energy_new = energy%total
610 
611  CALL timestop(handle)
612 
613  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  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  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  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  CALL timeset(routinen, handle)
647 
648  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  matrix_ks=matrix_ks)
655 
656  CALL qs_subsys_get(subsys, particles=particles)
657 
658  ! convert the density matrix (ls style) to QS style
659  CALL dbcsr_copy(matrix_p_qs, matrix_ks(1)%matrix)
660  CALL dbcsr_set(matrix_p_qs, 0.0_dp) !zero matrix creation
661  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  pw_pools=pw_pools)
667  CALL auxbas_pw_pool%create_pw(pw=wf_r)
668  CALL pw_zero(wf_r)
669  CALL auxbas_pw_pool%create_pw(pw=wf_g)
670  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  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  particles=particles, stride=stride)
679 
680  !free memory
681  CALL auxbas_pw_pool%give_back_pw(wf_r)
682  CALL auxbas_pw_pool%give_back_pw(wf_g)
683  CALL dbcsr_release(matrix_p_qs)
684 
685  CALL timestop(handle)
686 
687  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  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  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
704 
705  CALL timeset(routinen, handle)
706 
707  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  mixing_store=ls_scf_env%mixing_store)
711  IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
712  IF (dft_control%qs_control%gapw) THEN
713  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  ls_scf_env%para_env, rho_atom=rho_atom)
716  ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
717  CALL charge_mixing_init(ls_scf_env%mixing_store)
718  ELSEIF (dft_control%qs_control%semi_empirical) THEN
719  cpabort('SE Code not possible')
720  ELSE
721  CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
722  ls_scf_env%para_env)
723  END IF
724  END IF
725  CALL timestop(handle)
726  END SUBROUTINE rho_mixing_ls_init
727 
728 END MODULE dm_ls_scf_qs
Define the atomic kind types and their sub types.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
Definition: dm_ls_scf_qs.F:15
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
Definition: dm_ls_scf_qs.F:307
subroutine, public matrix_decluster(matrix_out, matrix_in, ls_mstruct)
Reverses molecular blocking and reduction to single precision if enabled.
Definition: dm_ls_scf_qs.F:374
subroutine, public ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
use the density matrix in ls_scf_env to compute the new energy and KS matrix
Definition: dm_ls_scf_qs.F:546
subroutine, public write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
...
Definition: dm_ls_scf_qs.F:625
subroutine, public rho_mixing_ls_init(qs_env, ls_scf_env)
Initialize g-space density mixing.
Definition: dm_ls_scf_qs.F:695
subroutine, public matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
create a matrix for use (and as a template) in ls based on a qs template
Definition: dm_ls_scf_qs.F:97
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
Definition: dm_ls_scf_qs.F:213
subroutine, public ls_scf_init_qs(qs_env)
further required initialization of QS. Might be factored-out since this seems common code with the ot...
Definition: dm_ls_scf_qs.F:408
subroutine, public ls_scf_qs_atomic_guess(qs_env, energy)
get an atomic initial guess
Definition: dm_ls_scf_qs.F:465
Types needed for a linear scaling quickstep SCF run based on the density matrix.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ls_cluster_molecular
integer, parameter, public ls_cluster_atomic
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Interface to the message passing library MPI.
represent a simple array based list of the given type
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
Routine to return block diagonal density matrix. Blocks correspond to the atomic densities.
subroutine, public calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, ounit, para_env)
returns a block diagonal density matrix. Blocks correspond to the atomic densities.
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
module that contains the definitions of the scf types
integer, parameter, public direct_mixing_nr
integer, parameter, public gspace_mixing_nr
Perform a QUICKSTEP wavefunction optimization (single point)
Definition: qs_energy.F:14
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
Routines to somehow generate an initial guess.
subroutine, public calculate_mopac_dm(pmat, matrix_s, has_unit_metric, dft_control, particle_set, atomic_kind_set, qs_kind_set, nspin, nelectron_spin, para_env)
returns a block diagonal density matrix. Blocks correspond to the mopac initial guess.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition: qs_ks_types.F:919
elemental subroutine, public charge_mixing_init(mixing_store)
initialiation needed when charge mixing is used
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
initialiation needed when gspace mixing is used
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
Define the neighbor list data types and the corresponding functionality.
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...