(git:ed6f26b)
Loading...
Searching...
No Matches
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-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines for a linear scaling quickstep SCF run based on the density
10!> matrix, with a focus on the interface between dm_ls_scf and qs
11!> \par History
12!> 2011.04 created [Joost VandeVondele]
13!> \author Joost VandeVondele
14! **************************************************************************************************
18 USE cp_dbcsr_api, ONLY: &
29 USE dm_ls_scf_types, ONLY: ls_cluster_atomic,&
30 ls_cluster_molecular,&
35 USE kinds, ONLY: default_string_length,&
36 dp
40 USE pw_env_types, ONLY: pw_env_get,&
42 USE pw_methods, ONLY: pw_zero
43 USE pw_pool_types, ONLY: pw_pool_p_type,&
45 USE pw_types, ONLY: pw_c1d_gs_type,&
61 USE qs_ks_types, ONLY: qs_ks_did_change,&
70 USE qs_rho_types, ONLY: qs_rho_get,&
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77
78 PRIVATE
79
80 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_qs'
81
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief create a matrix for use (and as a template) in ls based on a qs template
90!> \param matrix_ls ...
91!> \param matrix_qs ...
92!> \param ls_mstruct ...
93!> \par History
94!> 2011.03 created [Joost VandeVondele]
95!> 2015.09 add support for PAO [Ole Schuett]
96!> \author Joost VandeVondele
97! **************************************************************************************************
98 SUBROUTINE matrix_ls_create(matrix_ls, matrix_qs, ls_mstruct)
99 TYPE(dbcsr_type) :: matrix_ls, matrix_qs
100 TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
101
102 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_ls_create'
103
104 CHARACTER(len=default_string_length) :: name
105 INTEGER :: handle, iatom, imol, jatom, natom, nmol
106 INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: atom_to_cluster, atom_to_cluster_primus, &
107 clustered_blk_sizes, primus_of_mol
108 INTEGER, DIMENSION(:), POINTER :: clustered_col_dist, clustered_row_dist, &
109 ls_blk_sizes, ls_col_dist, ls_row_dist
110 TYPE(dbcsr_distribution_type) :: ls_dist, ls_dist_clustered
111
112 CALL timeset(routinen, handle)
113
114 ! Defaults -----------------------------------------------------------------------------------
115 CALL dbcsr_get_info(matrix_qs, col_blk_size=ls_blk_sizes, distribution=ls_dist)
116 CALL dbcsr_distribution_hold(ls_dist)
117 CALL dbcsr_distribution_get(ls_dist, row_dist=ls_row_dist, col_dist=ls_col_dist)
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
129 ! create format of the clustered matrix
130 CALL dbcsr_get_info(matrix_qs, nblkrows_total=natom)
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 row_blk_size=ls_blk_sizes, &
190 col_blk_size=ls_blk_sizes)
191 CALL dbcsr_distribution_release(ls_dist)
192 CALL dbcsr_finalize(matrix_ls)
193
194 CALL timestop(handle)
195
196 END SUBROUTINE matrix_ls_create
197
198! **************************************************************************************************
199!> \brief first link to QS, copy a QS matrix to LS matrix
200!> used to isolate QS style matrices from LS style
201!> will be useful for future features (e.g. precision, symmetry, blocking, ...)
202!> \param matrix_ls ...
203!> \param matrix_qs ...
204!> \param ls_mstruct ...
205!> \param covariant ...
206!> \par History
207!> 2010.10 created [Joost VandeVondele]
208!> 2015.09 add support for PAO [Ole Schuett]
209!> \author Joost VandeVondele
210! **************************************************************************************************
211 SUBROUTINE matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
212 TYPE(dbcsr_type) :: matrix_ls, matrix_qs
213 TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
214 LOGICAL, INTENT(IN) :: covariant
215
216 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_qs_to_ls'
217
218 INTEGER :: handle
219 INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
220 TYPE(dbcsr_type) :: matrix_pao, matrix_tmp
221 TYPE(dbcsr_type), POINTER :: matrix_trafo
222
223 CALL timeset(routinen, handle)
224
225 IF (.NOT. ls_mstruct%do_pao) THEN
226 CALL matrix_cluster(matrix_ls, matrix_qs, ls_mstruct)
227
228 ELSE ! using pao
229 CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
230 CALL dbcsr_create(matrix_pao, &
231 matrix_type="N", &
232 template=matrix_qs, &
233 row_blk_size=pao_blk_sizes, &
234 col_blk_size=pao_blk_sizes)
235
236 matrix_trafo => ls_mstruct%matrix_A ! contra-variant
237 IF (covariant) matrix_trafo => ls_mstruct%matrix_B ! co-variant
238 CALL dbcsr_create(matrix_tmp, template=matrix_trafo)
239
240 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qs, matrix_trafo, 0.0_dp, matrix_tmp)
241 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_trafo, matrix_tmp, 0.0_dp, matrix_pao)
242 CALL dbcsr_release(matrix_tmp)
243
244 CALL matrix_cluster(matrix_ls, matrix_pao, ls_mstruct)
245 CALL dbcsr_release(matrix_pao)
246 END IF
247
248 CALL timestop(handle)
249
250 END SUBROUTINE matrix_qs_to_ls
251
252! **************************************************************************************************
253!> \brief Performs molecular blocking and reduction to single precision if enabled
254!> \param matrix_out ...
255!> \param matrix_in ...
256!> \param ls_mstruct ...
257!> \author Ole Schuett
258! **************************************************************************************************
259 SUBROUTINE matrix_cluster(matrix_out, matrix_in, ls_mstruct)
260 TYPE(dbcsr_type) :: matrix_out, matrix_in
261 TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
262
263 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_cluster'
264
265 INTEGER :: handle
266 TYPE(dbcsr_type) :: matrix_in_nosym
267
268 CALL timeset(routinen, handle)
269
270 SELECT CASE (ls_mstruct%cluster_type)
271 CASE (ls_cluster_atomic)
272 CALL dbcsr_copy(matrix_out, matrix_in)
273
275 ! desymmetrize the qs matrix
276 CALL dbcsr_create(matrix_in_nosym, template=matrix_in, matrix_type="N")
277 CALL dbcsr_desymmetrize(matrix_in, matrix_in_nosym)
278
279 ! perform the magic complete redistribute copy
280 CALL dbcsr_complete_redistribute(matrix_in_nosym, matrix_out)
281 CALL dbcsr_release(matrix_in_nosym)
282
283 CASE DEFAULT
284 cpabort("Unknown LS cluster type")
285 END SELECT
286
287 CALL timestop(handle)
288
289 END SUBROUTINE matrix_cluster
290
291! **************************************************************************************************
292!> \brief second link to QS, copy a LS matrix to QS matrix
293!> used to isolate QS style matrices from LS style
294!> will be useful for future features (e.g. precision, symmetry, blocking, ...)
295!> \param matrix_qs ...
296!> \param matrix_ls ...
297!> \param ls_mstruct ...
298!> \param covariant ...
299!> \param keep_sparsity will be passed on to dbcsr_copy, by default set to .TRUE.
300!> \par History
301!> 2010.10 created [Joost VandeVondele]
302!> 2015.09 add support for PAO [Ole Schuett]
303!> \author Joost VandeVondele
304! **************************************************************************************************
305 SUBROUTINE matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
306 TYPE(dbcsr_type) :: matrix_qs, matrix_ls
307 TYPE(ls_mstruct_type), INTENT(IN), TARGET :: ls_mstruct
308 LOGICAL :: covariant
309 LOGICAL, OPTIONAL :: keep_sparsity
310
311 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_ls_to_qs'
312
313 INTEGER :: handle
314 INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
315 LOGICAL :: my_keep_sparsity
316 TYPE(dbcsr_type) :: matrix_declustered, matrix_tmp1, &
317 matrix_tmp2
318 TYPE(dbcsr_type), POINTER :: matrix_trafo
319
320 CALL timeset(routinen, handle)
321
322 my_keep_sparsity = .true.
323 IF (PRESENT(keep_sparsity)) &
324 my_keep_sparsity = keep_sparsity
325
326 IF (.NOT. ls_mstruct%do_pao) THEN
327 CALL dbcsr_create(matrix_declustered, template=matrix_qs)
328 CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
329 CALL dbcsr_copy(matrix_qs, matrix_declustered, keep_sparsity=my_keep_sparsity)
330 CALL dbcsr_release(matrix_declustered)
331
332 ELSE ! using pao
333 CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
334 CALL dbcsr_create(matrix_declustered, &
335 template=matrix_qs, &
336 row_blk_size=pao_blk_sizes, &
337 col_blk_size=pao_blk_sizes)
338
339 CALL matrix_decluster(matrix_declustered, matrix_ls, ls_mstruct)
340
341 matrix_trafo => ls_mstruct%matrix_B ! contra-variant
342 IF (covariant) matrix_trafo => ls_mstruct%matrix_A ! co-variant
343 CALL dbcsr_create(matrix_tmp1, template=matrix_trafo)
344 CALL dbcsr_create(matrix_tmp2, template=matrix_qs)
345 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_declustered, 0.0_dp, matrix_tmp1)
346 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp1, matrix_trafo, 0.0_dp, matrix_tmp2)
347 CALL dbcsr_copy(matrix_qs, matrix_tmp2, keep_sparsity=my_keep_sparsity)
348 CALL dbcsr_release(matrix_declustered)
349 CALL dbcsr_release(matrix_tmp1)
350 CALL dbcsr_release(matrix_tmp2)
351 END IF
352
353 CALL timestop(handle)
354
355 END SUBROUTINE matrix_ls_to_qs
356
357! **************************************************************************************************
358!> \brief Reverses molecular blocking and reduction to single precision if enabled
359!> \param matrix_out ...
360!> \param matrix_in ...
361!> \param ls_mstruct ...
362!> \author Ole Schuett
363! **************************************************************************************************
364 SUBROUTINE matrix_decluster(matrix_out, matrix_in, ls_mstruct)
365 TYPE(dbcsr_type) :: matrix_out, matrix_in
366 TYPE(ls_mstruct_type), INTENT(IN) :: ls_mstruct
367
368 CHARACTER(len=*), PARAMETER :: routinen = 'matrix_decluster'
369
370 INTEGER :: handle
371
372 CALL timeset(routinen, handle)
373
374 SELECT CASE (ls_mstruct%cluster_type)
375 CASE (ls_cluster_atomic)
376 CALL dbcsr_copy(matrix_out, matrix_in)
377
379 ! perform the magic complete redistribute copy
380 CALL dbcsr_complete_redistribute(matrix_in, matrix_out)
381
382 CASE DEFAULT
383 cpabort("Unknown LS cluster type")
384 END SELECT
385
386 CALL timestop(handle)
387
388 END SUBROUTINE matrix_decluster
389
390! **************************************************************************************************
391!> \brief further required initialization of QS.
392!> Might be factored-out since this seems common code with the other SCF.
393!> \param qs_env ...
394!> \par History
395!> 2010.10 created [Joost VandeVondele]
396!> \author Joost VandeVondele
397! **************************************************************************************************
398 SUBROUTINE ls_scf_init_qs(qs_env)
399 TYPE(qs_environment_type), POINTER :: qs_env
400
401 CHARACTER(len=*), PARAMETER :: routinen = 'ls_scf_init_qs'
402
403 INTEGER :: handle, ispin, nspin, unit_nr
404 TYPE(cp_logger_type), POINTER :: logger
405 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
406 TYPE(dft_control_type), POINTER :: dft_control
407 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
408 POINTER :: sab_orb
409 TYPE(qs_ks_env_type), POINTER :: ks_env
410
411 NULLIFY (sab_orb)
412 CALL timeset(routinen, handle)
413
414 ! get a useful output_unit
415 logger => cp_get_default_logger()
416 IF (logger%para_env%is_source()) THEN
417 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
418 ELSE
419 unit_nr = -1
420 END IF
421
422 ! get basic quantities from the qs_env
423 CALL get_qs_env(qs_env, dft_control=dft_control, &
424 matrix_s=matrix_s, &
425 matrix_ks=matrix_ks, &
426 ks_env=ks_env, &
427 sab_orb=sab_orb)
428
429 nspin = dft_control%nspins
430
431 ! we might have to create matrix_ks
432 IF (.NOT. ASSOCIATED(matrix_ks)) THEN
433 CALL dbcsr_allocate_matrix_set(matrix_ks, nspin)
434 DO ispin = 1, nspin
435 ALLOCATE (matrix_ks(ispin)%matrix)
436 CALL dbcsr_create(matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
437 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix, sab_orb)
438 CALL dbcsr_set(matrix_ks(ispin)%matrix, 0.0_dp)
439 END DO
440 CALL set_ks_env(ks_env, matrix_ks=matrix_ks)
441 END IF
442
443 CALL timestop(handle)
444
445 END SUBROUTINE ls_scf_init_qs
446
447! **************************************************************************************************
448!> \brief get an atomic initial guess
449!> \param qs_env ...
450!> \param ls_scf_env ...
451!> \param energy ...
452!> \param nonscf ...
453!> \par History
454!> 2012.11 created [Joost VandeVondele]
455!> \author Joost VandeVondele
456! **************************************************************************************************
457 SUBROUTINE ls_scf_qs_atomic_guess(qs_env, ls_scf_env, energy, nonscf)
458 TYPE(qs_environment_type), POINTER :: qs_env
459 TYPE(ls_scf_env_type) :: ls_scf_env
460 REAL(kind=dp) :: energy
461 LOGICAL, INTENT(IN), OPTIONAL :: nonscf
462
463 CHARACTER(len=*), PARAMETER :: routinen = 'ls_scf_qs_atomic_guess'
464
465 INTEGER :: handle, nspin, unit_nr
466 INTEGER, DIMENSION(2) :: nelectron_spin
467 LOGICAL :: do_scf, has_unit_metric
468 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
469 TYPE(cp_logger_type), POINTER :: logger
470 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
471 TYPE(dft_control_type), POINTER :: dft_control
472 TYPE(mp_para_env_type), POINTER :: para_env
473 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
474 TYPE(qs_energy_type), POINTER :: qs_energy
475 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
476 TYPE(qs_ks_env_type), POINTER :: ks_env
477 TYPE(qs_rho_type), POINTER :: rho
478
479 CALL timeset(routinen, handle)
480 NULLIFY (rho, rho_ao)
481
482 ! get a useful output_unit
483 logger => cp_get_default_logger()
484 IF (logger%para_env%is_source()) THEN
485 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
486 ELSE
487 unit_nr = -1
488 END IF
489
490 ! get basic quantities from the qs_env
491 CALL get_qs_env(qs_env, dft_control=dft_control, &
492 matrix_s=matrix_s, &
493 matrix_ks=matrix_ks, &
494 ks_env=ks_env, &
495 energy=qs_energy, &
496 atomic_kind_set=atomic_kind_set, &
497 qs_kind_set=qs_kind_set, &
498 particle_set=particle_set, &
499 has_unit_metric=has_unit_metric, &
500 para_env=para_env, &
501 nelectron_spin=nelectron_spin, &
502 rho=rho)
503
504 CALL qs_rho_get(rho, rho_ao=rho_ao)
505
506 nspin = dft_control%nspins
507
508 ! create an initial atomic guess
509 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical .OR. &
510 dft_control%qs_control%xtb) THEN
511 CALL calculate_mopac_dm(rho_ao, matrix_s(1)%matrix, has_unit_metric, &
512 dft_control, particle_set, atomic_kind_set, qs_kind_set, &
513 nspin, nelectron_spin, para_env)
514 ELSE
515 CALL calculate_atomic_block_dm(rho_ao, matrix_s(1)%matrix, atomic_kind_set, qs_kind_set, &
516 nspin, nelectron_spin, unit_nr, para_env)
517 END IF
518
519 do_scf = .true.
520 IF (PRESENT(nonscf)) do_scf = .NOT. nonscf
521 IF (do_scf) THEN
522 CALL qs_rho_update_rho(rho, qs_env=qs_env)
523 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
524 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., just_energy=.false.)
525 energy = qs_energy%total
526 ELSE
527 CALL ls_nonscf_ks(qs_env, ls_scf_env, energy)
528 END IF
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 use the external density in ls_scf_env to compute the new KS matrix
617!> \param qs_env ...
618!> \param ls_scf_env ...
619!> \param energy_new ...
620! **************************************************************************************************
621 SUBROUTINE ls_nonscf_ks(qs_env, ls_scf_env, energy_new)
622 TYPE(qs_environment_type), POINTER :: qs_env
623 TYPE(ls_scf_env_type) :: ls_scf_env
624 REAL(kind=dp) :: energy_new
625
626 CHARACTER(len=*), PARAMETER :: routinen = 'ls_nonscf_ks'
627
628 INTEGER :: handle, ispin, nspin
629 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
630 TYPE(harris_type), POINTER :: harris_env
631 TYPE(mp_para_env_type), POINTER :: para_env
632 TYPE(qs_energy_type), POINTER :: energy
633 TYPE(qs_rho_type), POINTER :: rho
634
635 NULLIFY (energy, rho, rho_ao)
636 CALL timeset(routinen, handle)
637
638 nspin = ls_scf_env%nspins
639 CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
640 CALL qs_rho_get(rho, rho_ao=rho_ao)
641
642 ! set the new density matrix
643 DO ispin = 1, nspin
644 CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
645 ls_scf_env%ls_mstruct, covariant=.false.)
646 END DO
647
648 IF (qs_env%harris_method) THEN
649 CALL get_qs_env(qs_env, harris_env=harris_env)
650 CALL harris_density_update(qs_env, harris_env)
651 END IF
652 ! compute the corresponding KS matrix and new energy
653 CALL qs_rho_update_rho(rho, qs_env=qs_env)
654 IF (ls_scf_env%do_rho_mixing) THEN
655 cpabort("P mixing not implemented in linear scaling NONSCF. ")
656 END IF
657
658 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
659 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., &
660 just_energy=.false., print_active=.true.)
661 energy_new = energy%total
662
663 CALL timestop(handle)
664
665 END SUBROUTINE ls_nonscf_ks
666
667! **************************************************************************************************
668!> \brief use the new density matrix in ls_scf_env to compute the new energy
669!> \param qs_env ...
670!> \param ls_scf_env ...
671! **************************************************************************************************
672 SUBROUTINE ls_nonscf_energy(qs_env, ls_scf_env)
673 TYPE(qs_environment_type), POINTER :: qs_env
674 TYPE(ls_scf_env_type) :: ls_scf_env
675
676 CHARACTER(len=*), PARAMETER :: routinen = 'ls_nonscf_energy'
677
678 INTEGER :: handle, ispin, nspin
679 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao
680 TYPE(mp_para_env_type), POINTER :: para_env
681 TYPE(qs_energy_type), POINTER :: energy
682 TYPE(qs_rho_type), POINTER :: rho
683
684 NULLIFY (energy, rho, rho_ao)
685 CALL timeset(routinen, handle)
686 IF (qs_env%qmmm) THEN
687 cpabort("NYA")
688 END IF
689
690 nspin = ls_scf_env%nspins
691 CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho)
692 CALL qs_rho_get(rho, rho_ao=rho_ao)
693
694 ! set the new density matrix
695 DO ispin = 1, nspin
696 CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), &
697 ls_scf_env%ls_mstruct, covariant=.false.)
698 END DO
699
700 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
701
702 ! band energy : Tr(PH)
703 CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
704 CALL calculate_ptrace(matrix_ks, rho_ao, energy%band, nspin)
705 ! core energy : Tr(Ph)
706 energy%total = energy%total - energy%core
707 CALL get_qs_env(qs_env, matrix_h=matrix_h)
708 CALL calculate_ptrace(matrix_h, rho_ao, energy%core, nspin)
709
710 CALL timestop(handle)
711
712 END SUBROUTINE ls_nonscf_energy
713
714! **************************************************************************************************
715!> \brief ...
716!> \param qs_env ...
717!> \param ls_scf_env ...
718!> \param matrix_p_ls ...
719!> \param unit_nr ...
720!> \param title ...
721!> \param stride ...
722! **************************************************************************************************
723 SUBROUTINE write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
724 TYPE(qs_environment_type), POINTER :: qs_env
725 TYPE(ls_scf_env_type) :: ls_scf_env
726 TYPE(dbcsr_type), INTENT(IN) :: matrix_p_ls
727 INTEGER, INTENT(IN) :: unit_nr
728 CHARACTER(LEN=*), INTENT(IN) :: title
729 INTEGER, DIMENSION(:), POINTER :: stride
730
731 CHARACTER(len=*), PARAMETER :: routinen = 'write_matrix_to_cube'
732
733 INTEGER :: handle
734 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
735 TYPE(dbcsr_type), TARGET :: matrix_p_qs
736 TYPE(particle_list_type), POINTER :: particles
737 TYPE(pw_c1d_gs_type) :: wf_g
738 TYPE(pw_env_type), POINTER :: pw_env
739 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
740 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
741 TYPE(pw_r3d_rs_type) :: wf_r
742 TYPE(qs_ks_env_type), POINTER :: ks_env
743 TYPE(qs_subsys_type), POINTER :: subsys
744
745 CALL timeset(routinen, handle)
746
747 NULLIFY (ks_env, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, matrix_ks)
748
749 CALL get_qs_env(qs_env, &
750 ks_env=ks_env, &
751 subsys=subsys, &
752 pw_env=pw_env, &
753 matrix_ks=matrix_ks)
754
755 CALL qs_subsys_get(subsys, particles=particles)
756
757 ! convert the density matrix (ls style) to QS style
758 CALL dbcsr_copy(matrix_p_qs, matrix_ks(1)%matrix)
759 CALL dbcsr_set(matrix_p_qs, 0.0_dp) !zero matrix creation
760 CALL matrix_ls_to_qs(matrix_p_qs, matrix_p_ls, ls_scf_env%ls_mstruct, covariant=.false.)
761
762 ! Print total electronic density
763 CALL pw_env_get(pw_env=pw_env, &
764 auxbas_pw_pool=auxbas_pw_pool, &
765 pw_pools=pw_pools)
766 CALL auxbas_pw_pool%create_pw(pw=wf_r)
767 CALL pw_zero(wf_r)
768 CALL auxbas_pw_pool%create_pw(pw=wf_g)
769 CALL pw_zero(wf_g)
770 CALL calculate_rho_elec(matrix_p=matrix_p_qs, &
771 rho=wf_r, &
772 rho_gspace=wf_g, &
773 ks_env=ks_env)
774
775 ! write this to a cube
776 CALL cp_pw_to_cube(wf_r, unit_nr=unit_nr, title=title, &
777 particles=particles, stride=stride)
778
779 !free memory
780 CALL auxbas_pw_pool%give_back_pw(wf_r)
781 CALL auxbas_pw_pool%give_back_pw(wf_g)
782 CALL dbcsr_release(matrix_p_qs)
783
784 CALL timestop(handle)
785
786 END SUBROUTINE write_matrix_to_cube
787
788! **************************************************************************************************
789!> \brief Initialize g-space density mixing
790!> \param qs_env ...
791!> \param ls_scf_env ...
792! **************************************************************************************************
793 SUBROUTINE rho_mixing_ls_init(qs_env, ls_scf_env)
794 TYPE(qs_environment_type), POINTER :: qs_env
795 TYPE(ls_scf_env_type) :: ls_scf_env
796
797 CHARACTER(len=*), PARAMETER :: routinen = 'rho_mixing_ls_init'
798
799 INTEGER :: handle
800 TYPE(dft_control_type), POINTER :: dft_control
801 TYPE(qs_rho_type), POINTER :: rho
802 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
803
804 CALL timeset(routinen, handle)
805
806 CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho)
807
808 CALL mixing_allocate(qs_env, ls_scf_env%density_mixing_method, nspins=ls_scf_env%nspins, &
809 mixing_store=ls_scf_env%mixing_store)
810 IF (ls_scf_env%density_mixing_method >= gspace_mixing_nr) THEN
811 IF (dft_control%qs_control%gapw) THEN
812 CALL get_qs_env(qs_env, rho_atom_set=rho_atom)
813 CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
814 ls_scf_env%para_env, rho_atom=rho_atom)
815 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
816 CALL charge_mixing_init(ls_scf_env%mixing_store)
817 ELSEIF (dft_control%qs_control%semi_empirical) THEN
818 cpabort('SE Code not possible')
819 ELSE
820 CALL mixing_init(ls_scf_env%density_mixing_method, rho, ls_scf_env%mixing_store, &
821 ls_scf_env%para_env)
822 END IF
823 END IF
824 CALL timestop(handle)
825 END SUBROUTINE rho_mixing_ls_init
826
827END 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...
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_distribution_hold(dist)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
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...
subroutine, public ls_scf_qs_atomic_guess(qs_env, ls_scf_env, energy, nonscf)
get an atomic initial guess
subroutine, public ls_nonscf_ks(qs_env, ls_scf_env, energy_new)
use the external density in ls_scf_env to compute the new KS matrix
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...
subroutine, public matrix_decluster(matrix_out, matrix_in, ls_mstruct)
Reverses molecular blocking and reduction to single precision if enabled.
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
subroutine, public write_matrix_to_cube(qs_env, ls_scf_env, matrix_p_ls, unit_nr, title, stride)
...
subroutine, public rho_mixing_ls_init(qs_env, ls_scf_env)
Initialize g-space density mixing.
subroutine, public ls_nonscf_energy(qs_env, ls_scf_env)
use the new density matrix in ls_scf_env to compute the new energy
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
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 ...
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...
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
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
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
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
Calculation of the energies concerning the core charge distribution.
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_pp, 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, harris_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, eeq, 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.
Types needed for a for a Harris model calculation.
Harris method environment setup and handling.
subroutine, public harris_density_update(qs_env, harris_env)
...
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.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
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_pp, 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)
...
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...
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...
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...
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)
...
Provides all information about an atomic kind.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
contained for different pw related things
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Contains information on the Harris method.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.