(git:15c1bfc)
Loading...
Searching...
No Matches
energy_corrections.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 an energy correction on top of a Kohn-Sham calculation
10!> \par History
11!> 03.2014 created
12!> 09.2019 Moved from KG to Kohn-Sham
13!> 08.2022 Add Density-Corrected DFT methods
14!> 04.2023 Add hybrid functionals for DC-DFT
15!> 10.2024 Add external energy method
16!> \author JGH
17! **************************************************************************************************
25 USE bibliography, ONLY: belleflamme2023,&
26 cite_reference
27 USE cell_types, ONLY: cell_type,&
28 pbc
31 USE cp_dbcsr_api, ONLY: &
34 dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
42 USE cp_files, ONLY: close_file,&
51 USE cp_fm_types, ONLY: cp_fm_create,&
60 USE cp_output_handling, ONLY: cp_p_file,&
84 USE ec_external, ONLY: ec_ext_energy
85 USE ec_methods, ONLY: create_kernel,&
91 USE hfx_exx, ONLY: add_exx_to_rhs,&
93 USE input_constants, ONLY: &
106 USE kinds, ONLY: default_path_length,&
108 dp
110 USE mathlib, ONLY: det_3x3
115 USE periodic_table, ONLY: ptable
116 USE physcon, ONLY: bohr,&
117 debye,&
118 pascal
123 USE pw_env_types, ONLY: pw_env_get,&
125 USE pw_grid_types, ONLY: pw_grid_type
126 USE pw_methods, ONLY: pw_axpy,&
127 pw_copy,&
129 pw_scale,&
131 pw_zero
134 USE pw_pool_types, ONLY: pw_pool_p_type,&
136 USE pw_types, ONLY: pw_c1d_gs_type,&
152 USE qs_force_types, ONLY: qs_force_type,&
155 USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
156 integrate_v_rspace
157 USE qs_kind_types, ONLY: get_qs_kind,&
163 USE qs_ks_types, ONLY: qs_ks_env_type
166 USE qs_mo_types, ONLY: deallocate_mo_set,&
167 get_mo_set,&
177 USE qs_oce_types, ONLY: allocate_oce_set,&
182 USE qs_rho_types, ONLY: qs_rho_get,&
184 USE qs_vxc, ONLY: qs_vxc_create
188 USE string_utilities, ONLY: uppercase
192 USE trexio_utils, ONLY: write_trexio
200#include "./base/base_uses.f90"
201
202 IMPLICIT NONE
203
204 PRIVATE
205
206 ! Global parameters
207
208 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'energy_corrections'
209
210 PUBLIC :: energy_correction
211
212CONTAINS
213
214! **************************************************************************************************
215!> \brief Energy Correction to a Kohn-Sham simulation
216!> Available energy corrections: (1) Harris energy functional
217!> (2) Density-corrected DFT
218!> (3) Energy from external source
219!>
220!> \param qs_env ...
221!> \param ec_init ...
222!> \param calculate_forces ...
223!> \par History
224!> 03.2014 created
225!> \author JGH
226! **************************************************************************************************
227 SUBROUTINE energy_correction(qs_env, ec_init, calculate_forces)
228 TYPE(qs_environment_type), POINTER :: qs_env
229 LOGICAL, INTENT(IN), OPTIONAL :: ec_init, calculate_forces
230
231 CHARACTER(len=*), PARAMETER :: routinen = 'energy_correction'
232
233 INTEGER :: handle, unit_nr
234 LOGICAL :: my_calc_forces
235 TYPE(cp_logger_type), POINTER :: logger
236 TYPE(energy_correction_type), POINTER :: ec_env
237 TYPE(qs_energy_type), POINTER :: energy
238 TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force
239 TYPE(virial_type), POINTER :: virial
240
241 CALL timeset(routinen, handle)
242
243 logger => cp_get_default_logger()
244 IF (logger%para_env%is_source()) THEN
245 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
246 ELSE
247 unit_nr = -1
248 END IF
249
250 CALL cite_reference(belleflamme2023)
251
252 NULLIFY (ec_env)
253 CALL get_qs_env(qs_env, ec_env=ec_env)
254
255 ! Skip energy correction if ground-state is NOT converged
256 IF (.NOT. ec_env%do_skip) THEN
257
258 ec_env%should_update = .true.
259 IF (PRESENT(ec_init)) ec_env%should_update = ec_init
260
261 my_calc_forces = .false.
262 IF (PRESENT(calculate_forces)) my_calc_forces = calculate_forces
263
264 IF (ec_env%should_update) THEN
265 ec_env%old_etotal = 0.0_dp
266 ec_env%etotal = 0.0_dp
267 ec_env%eband = 0.0_dp
268 ec_env%ehartree = 0.0_dp
269 ec_env%ex = 0.0_dp
270 ec_env%exc = 0.0_dp
271 ec_env%vhxc = 0.0_dp
272 ec_env%edispersion = 0.0_dp
273 ec_env%exc_aux_fit = 0.0_dp
274
275 ! Save total energy of reference calculation
276 CALL get_qs_env(qs_env, energy=energy)
277 ec_env%old_etotal = energy%total
278
279 END IF
280
281 IF (my_calc_forces) THEN
282 IF (unit_nr > 0) THEN
283 WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 25), &
284 " Energy Correction Forces ", repeat("-", 26), "!"
285 END IF
286 CALL get_qs_env(qs_env, force=ks_force, virial=virial)
287 CALL zero_qs_force(ks_force)
288 CALL zero_virial(virial, reset=.false.)
289 ELSE
290 IF (unit_nr > 0) THEN
291 WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 29), &
292 " Energy Correction ", repeat("-", 29), "!"
293 END IF
294 END IF
295
296 ! Perform the energy correction
297 CALL energy_correction_low(qs_env, ec_env, my_calc_forces, unit_nr)
298
299 ! Update total energy in qs environment and amount fo correction
300 IF (ec_env%should_update) THEN
301 energy%nonscf_correction = ec_env%etotal - ec_env%old_etotal
302 energy%total = ec_env%etotal
303 END IF
304
305 IF (.NOT. my_calc_forces .AND. unit_nr > 0) THEN
306 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Energy Correction ", energy%nonscf_correction
307 END IF
308 IF (unit_nr > 0) THEN
309 WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
310 END IF
311
312 ELSE
313
314 ! Ground-state energy calculation did not converge,
315 ! do not calculate energy correction
316 IF (unit_nr > 0) THEN
317 WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
318 WRITE (unit_nr, '(T2,A,A,A,A,A)') "!", repeat("-", 26), &
319 " Skip Energy Correction ", repeat("-", 27), "!"
320 WRITE (unit_nr, '(T2,A,A,A)') "!", repeat("-", 77), "!"
321 END IF
322
323 END IF
324
325 CALL timestop(handle)
326
327 END SUBROUTINE energy_correction
328
329! **************************************************************************************************
330!> \brief Energy Correction to a Kohn-Sham simulation
331!>
332!> \param qs_env ...
333!> \param ec_env ...
334!> \param calculate_forces ...
335!> \param unit_nr ...
336!> \par History
337!> 03.2014 created
338!> \author JGH
339! **************************************************************************************************
340 SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr)
341 TYPE(qs_environment_type), POINTER :: qs_env
342 TYPE(energy_correction_type), POINTER :: ec_env
343 LOGICAL, INTENT(IN) :: calculate_forces
344 INTEGER, INTENT(IN) :: unit_nr
345
346 INTEGER :: ispin, nkind, nspins
347 LOGICAL :: debug_f, gapw, gapw_xc
348 REAL(kind=dp) :: eps_fit, exc
349 TYPE(dft_control_type), POINTER :: dft_control
350 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
351 POINTER :: sap_oce
352 TYPE(oce_matrix_type), POINTER :: oce
353 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
354 TYPE(pw_env_type), POINTER :: pw_env
355 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
356 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
357
358 IF (ec_env%should_update) THEN
359 CALL ec_build_neighborlist(qs_env, ec_env)
360 CALL ks_ref_potential(qs_env, &
361 ec_env%vh_rspace, &
362 ec_env%vxc_rspace, &
363 ec_env%vtau_rspace, &
364 ec_env%vadmm_rspace, &
365 ec_env%ehartree, exc)
366
367 SELECT CASE (ec_env%energy_functional)
369
370 CALL ec_build_core_hamiltonian(qs_env, ec_env)
371 CALL ec_build_ks_matrix(qs_env, ec_env)
372
373 IF (ec_env%mao) THEN
374 ! MAO basis
375 IF (ASSOCIATED(ec_env%mao_coef)) CALL dbcsr_deallocate_matrix_set(ec_env%mao_coef)
376 NULLIFY (ec_env%mao_coef)
377 CALL mao_generate_basis(qs_env, ec_env%mao_coef, ref_basis_set="HARRIS", molecular=.true., &
378 max_iter=ec_env%mao_max_iter, eps_grad=ec_env%mao_eps_grad, &
379 eps1_mao=ec_env%mao_eps1, iolevel=ec_env%mao_iolevel, unit_nr=unit_nr)
380 END IF
381
382 CALL ec_ks_solver(qs_env, ec_env)
383
384 CALL evaluate_ec_core_matrix_traces(qs_env, ec_env)
385
386 CASE (ec_functional_dc)
387
388 ! Prepare Density-corrected DFT (DC-DFT) calculation
389 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.false.)
390
391 ! Rebuild KS matrix with DC-DFT XC functional evaluated in ground-state density.
392 ! KS matrix might contain unwanted contributions
393 ! Calculate Hartree and XC related energies here
394 CALL ec_build_ks_matrix(qs_env, ec_env)
395
396 CASE (ec_functional_ext)
397
398 CALL ec_ext_energy(qs_env, ec_env, calculate_forces=.false.)
399
400 CASE DEFAULT
401 cpabort("unknown energy correction")
402 END SELECT
403
404 ! dispersion through pairpotentials
405 CALL ec_disp(qs_env, ec_env, calculate_forces=.false.)
406
407 ! Calculate total energy
408 CALL ec_energy(ec_env, unit_nr)
409
410 END IF
411
412 IF (calculate_forces) THEN
413
414 debug_f = ec_env%debug_forces .OR. ec_env%debug_stress
415
416 CALL ec_disp(qs_env, ec_env, calculate_forces=.true.)
417
418 SELECT CASE (ec_env%energy_functional)
420
421 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
422 ec_env%matrix_p, &
423 ec_env%matrix_s, &
424 ec_env%matrix_w)
425 CALL ec_build_ks_matrix_force(qs_env, ec_env)
426 IF (ec_env%debug_external) THEN
427 CALL write_response_interface(qs_env, ec_env)
428 CALL init_response_deriv(qs_env, ec_env)
429 END IF
430
431 CASE (ec_functional_dc)
432
433 ! Prepare Density-corrected DFT (DC-DFT) calculation
434 ! by getting ground-state matrices
435 CALL ec_dc_energy(qs_env, ec_env, calculate_forces=.true.)
436
437 CALL ec_build_core_hamiltonian_force(qs_env, ec_env, &
438 ec_env%matrix_p, &
439 ec_env%matrix_s, &
440 ec_env%matrix_w)
441 CALL ec_dc_build_ks_matrix_force(qs_env, ec_env)
442 IF (ec_env%debug_external) THEN
443 CALL write_response_interface(qs_env, ec_env)
444 CALL init_response_deriv(qs_env, ec_env)
445 END IF
446
447 CASE (ec_functional_ext)
448
449 CALL init_response_deriv(qs_env, ec_env)
450 CALL ec_ext_energy(qs_env, ec_env, calculate_forces=.true.)
451
452 CASE DEFAULT
453 cpabort("unknown energy correction")
454 END SELECT
455
456 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
457 nspins = dft_control%nspins
458 gapw = dft_control%qs_control%gapw
459 gapw_xc = dft_control%qs_control%gapw_xc
460 IF (gapw .OR. gapw_xc) THEN
461 CALL get_qs_env(qs_env=qs_env, nkind=nkind, &
462 qs_kind_set=qs_kind_set, particle_set=particle_set)
463 CALL get_qs_env(qs_env=qs_env, oce=oce, sap_oce=sap_oce)
464 CALL create_oce_set(oce)
465 CALL allocate_oce_set(oce, nkind)
466 eps_fit = dft_control%qs_control%gapw_control%eps_fit
467 CALL build_oce_matrices(oce%intac, .true., 1, qs_kind_set, particle_set, &
468 sap_oce, eps_fit)
469 CALL set_qs_env(qs_env, oce=oce)
470 END IF
471
472 CALL response_calculation(qs_env, ec_env)
473
474 ! Allocate response density on real space grid for use in properties
475 ! Calculated in response_force
476 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
477
478 cpassert(ASSOCIATED(pw_env))
479 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
480 ALLOCATE (ec_env%rhoz_r(nspins))
481 DO ispin = 1, nspins
482 CALL auxbas_pw_pool%create_pw(ec_env%rhoz_r(ispin))
483 END DO
484
485 CALL response_force(qs_env, &
486 vh_rspace=ec_env%vh_rspace, &
487 vxc_rspace=ec_env%vxc_rspace, &
488 vtau_rspace=ec_env%vtau_rspace, &
489 vadmm_rspace=ec_env%vadmm_rspace, &
490 matrix_hz=ec_env%matrix_hz, &
491 matrix_pz=ec_env%matrix_z, &
492 matrix_pz_admm=ec_env%z_admm, &
493 matrix_wz=ec_env%matrix_wz, &
494 rhopz_r=ec_env%rhoz_r, &
495 zehartree=ec_env%ehartree, &
496 zexc=ec_env%exc, &
497 zexc_aux_fit=ec_env%exc_aux_fit, &
498 p_env=ec_env%p_env, &
499 debug=debug_f)
500
501 CALL output_response_deriv(qs_env, ec_env, unit_nr)
502
503 CALL ec_properties(qs_env, ec_env)
504
505 ! Deallocate Harris density and response density on grid
506 IF (ASSOCIATED(ec_env%rhoout_r)) THEN
507 DO ispin = 1, nspins
508 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin))
509 END DO
510 DEALLOCATE (ec_env%rhoout_r)
511 END IF
512 IF (ASSOCIATED(ec_env%rhoz_r)) THEN
513 DO ispin = 1, nspins
514 CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin))
515 END DO
516 DEALLOCATE (ec_env%rhoz_r)
517 END IF
518
519 ! Deallocate matrices
520 IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
521 IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
522 IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
523 IF (ASSOCIATED(ec_env%matrix_t)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_t)
524 IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
525 IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
526 IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
527 IF (ASSOCIATED(ec_env%matrix_wz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_wz)
528 IF (ASSOCIATED(ec_env%matrix_z)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_z)
529
530 END IF
531
532 END SUBROUTINE energy_correction_low
533
534! **************************************************************************************************
535!> \brief Output response information to TREXIO file (for testing external method read in)
536!> \param qs_env ...
537!> \param ec_env ...
538!> \author JHU
539! **************************************************************************************************
540 SUBROUTINE write_response_interface(qs_env, ec_env)
541 TYPE(qs_environment_type), POINTER :: qs_env
542 TYPE(energy_correction_type), POINTER :: ec_env
543
544 TYPE(section_vals_type), POINTER :: section, trexio_section
545
546 section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%TREXIO")
547 NULLIFY (trexio_section)
548 CALL section_vals_duplicate(section, trexio_section)
549 CALL section_vals_val_set(trexio_section, "FILENAME", c_val=ec_env%exresp_fn)
550 CALL section_vals_val_set(trexio_section, "CARTESIAN", l_val=.false.)
551 CALL write_trexio(qs_env, trexio_section, ec_env%matrix_hz)
552
553 END SUBROUTINE write_response_interface
554
555! **************************************************************************************************
556!> \brief Initialize arrays for response derivatives
557!> \param qs_env ...
558!> \param ec_env ...
559!> \author JHU
560! **************************************************************************************************
561 SUBROUTINE init_response_deriv(qs_env, ec_env)
562 TYPE(qs_environment_type), POINTER :: qs_env
563 TYPE(energy_correction_type), POINTER :: ec_env
564
565 INTEGER :: natom
566 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
567 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
568 TYPE(virial_type), POINTER :: virial
569
570 CALL get_qs_env(qs_env, natom=natom)
571 ALLOCATE (ec_env%rf(3, natom))
572 ec_env%rf = 0.0_dp
573 ec_env%rpv = 0.0_dp
574 CALL get_qs_env(qs_env, force=force, virial=virial)
575
576 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
577 CALL total_qs_force(ec_env%rf, force, atomic_kind_set)
578
579 IF (virial%pv_availability .AND. (.NOT. virial%pv_numer)) THEN
580 ec_env%rpv = virial%pv_virial
581 END IF
582
583 END SUBROUTINE init_response_deriv
584
585! **************************************************************************************************
586!> \brief Write the reponse forces to file
587!> \param qs_env ...
588!> \param ec_env ...
589!> \param unit_nr ...
590!> \author JHU
591! **************************************************************************************************
592 SUBROUTINE output_response_deriv(qs_env, ec_env, unit_nr)
593 TYPE(qs_environment_type), POINTER :: qs_env
594 TYPE(energy_correction_type), POINTER :: ec_env
595 INTEGER, INTENT(IN) :: unit_nr
596
597 CHARACTER(LEN=default_string_length) :: unit_string
598 INTEGER :: funit, ia, natom
599 REAL(kind=dp) :: evol, fconv
600 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
601 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
602 TYPE(cell_type), POINTER :: cell
603 TYPE(mp_para_env_type), POINTER :: para_env
604 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
605 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
606 TYPE(virial_type), POINTER :: virial
607
608 IF (ASSOCIATED(ec_env%rf)) THEN
609 CALL get_qs_env(qs_env, natom=natom)
610 ALLOCATE (ftot(3, natom))
611 ftot = 0.0_dp
612 CALL get_qs_env(qs_env, force=force, virial=virial, para_env=para_env)
613
614 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
615 CALL total_qs_force(ftot, force, atomic_kind_set)
616 ec_env%rf(1:3, 1:natom) = ftot(1:3, 1:natom) - ec_env%rf(1:3, 1:natom)
617 CALL para_env%sum(ec_env%rf)
618 DEALLOCATE (ftot)
619
620 IF (virial%pv_availability .AND. (.NOT. virial%pv_numer)) THEN
621 ec_env%rpv = virial%pv_virial - ec_env%rpv
622 CALL para_env%sum(ec_env%rpv)
623 ! Volume terms
624 evol = ec_env%exc + ec_env%exc_aux_fit + 2.0_dp*ec_env%ehartree
625 ec_env%rpv(1, 1) = ec_env%rpv(1, 1) - evol
626 ec_env%rpv(2, 2) = ec_env%rpv(2, 2) - evol
627 ec_env%rpv(3, 3) = ec_env%rpv(3, 3) - evol
628 END IF
629
630 CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
631 ! Conversion factor a.u. -> GPa
632 unit_string = "GPa"
633 fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, trim(unit_string))
634 IF (unit_nr > 0) THEN
635 WRITE (unit_nr, '(T2,A)') "Write EXTERNAL Response Derivative: "//trim(ec_env%exresult_fn)
636
637 CALL open_file(ec_env%exresult_fn, file_status="REPLACE", file_form="FORMATTED", &
638 file_action="WRITE", unit_number=funit)
639 WRITE (funit, "(T8,A,T58,A)") "COORDINATES [Bohr]", "RESPONSE FORCES [Hartree/Bohr]"
640 DO ia = 1, natom
641 WRITE (funit, "(2(3F15.8,5x))") particle_set(ia)%r(1:3), ec_env%rf(1:3, ia)
642 END DO
643 WRITE (funit, *)
644 WRITE (funit, "(T8,A,T58,A)") "CELL [Bohr]", "RESPONSE PRESSURE [GPa]"
645 DO ia = 1, 3
646 WRITE (funit, "(3F15.8,5x,3F15.8)") cell%hmat(ia, 1:3), -fconv*ec_env%rpv(ia, 1:3)
647 END DO
648
649 CALL close_file(funit)
650 END IF
651 END IF
652
653 END SUBROUTINE output_response_deriv
654
655! **************************************************************************************************
656!> \brief Calculates the traces of the core matrices and the density matrix.
657!> \param qs_env ...
658!> \param ec_env ...
659!> \author Ole Schuett
660!> adapted for energy correction fbelle
661! **************************************************************************************************
662 SUBROUTINE evaluate_ec_core_matrix_traces(qs_env, ec_env)
663 TYPE(qs_environment_type), POINTER :: qs_env
664 TYPE(energy_correction_type), POINTER :: ec_env
665
666 CHARACTER(LEN=*), PARAMETER :: routinen = 'evaluate_ec_core_matrix_traces'
667
668 INTEGER :: handle
669 TYPE(dft_control_type), POINTER :: dft_control
670 TYPE(qs_energy_type), POINTER :: energy
671
672 CALL timeset(routinen, handle)
673 NULLIFY (energy)
674
675 CALL get_qs_env(qs_env, dft_control=dft_control, energy=energy)
676
677 ! Core hamiltonian energy
678 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, energy%core, dft_control%nspins)
679
680 ! kinetic energy
681 CALL calculate_ptrace(ec_env%matrix_t, ec_env%matrix_p, energy%kinetic, dft_control%nspins)
682
683 CALL timestop(handle)
684
685 END SUBROUTINE evaluate_ec_core_matrix_traces
686
687! **************************************************************************************************
688!> \brief Prepare DC-DFT calculation by copying unaffected ground-state matrices (core Hamiltonian,
689!> density matrix) into energy correction environment and rebuild the overlap matrix
690!>
691!> \param qs_env ...
692!> \param ec_env ...
693!> \param calculate_forces ...
694!> \par History
695!> 07.2022 created
696!> \author fbelle
697! **************************************************************************************************
698 SUBROUTINE ec_dc_energy(qs_env, ec_env, calculate_forces)
699 TYPE(qs_environment_type), POINTER :: qs_env
700 TYPE(energy_correction_type), POINTER :: ec_env
701 LOGICAL, INTENT(IN) :: calculate_forces
702
703 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_dc_energy'
704
705 CHARACTER(LEN=default_string_length) :: headline
706 INTEGER :: handle, ispin, nspins
707 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w
708 TYPE(dft_control_type), POINTER :: dft_control
709 TYPE(qs_ks_env_type), POINTER :: ks_env
710 TYPE(qs_rho_type), POINTER :: rho
711
712 CALL timeset(routinen, handle)
713
714 NULLIFY (dft_control, ks_env, matrix_h, matrix_p, matrix_s, matrix_w, rho)
715 CALL get_qs_env(qs_env=qs_env, &
716 dft_control=dft_control, &
717 ks_env=ks_env, &
718 matrix_h_kp=matrix_h, &
719 matrix_s_kp=matrix_s, &
720 matrix_w_kp=matrix_w, &
721 rho=rho)
722 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
723 nspins = dft_control%nspins
724
725 ! For density-corrected DFT only the ground-state matrices are required
726 ! Comply with ec_env environment for property calculations later
727 CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
728 matrix_name="OVERLAP MATRIX", &
729 basis_type_a="HARRIS", &
730 basis_type_b="HARRIS", &
731 sab_nl=ec_env%sab_orb)
732
733 ! Core Hamiltonian matrix
734 IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
735 CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
736 headline = "CORE HAMILTONIAN MATRIX"
737 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
738 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=trim(headline), &
739 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
740 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, ec_env%sab_orb)
741 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
742
743 ! Density matrix
744 IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
745 CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
746 headline = "DENSITY MATRIX"
747 DO ispin = 1, nspins
748 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
749 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
750 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
751 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
752 CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
753 END DO
754
755 IF (calculate_forces) THEN
756
757 ! Energy-weighted density matrix
758 IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
759 CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
760 headline = "ENERGY-WEIGHTED DENSITY MATRIX"
761 DO ispin = 1, nspins
762 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
763 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
764 template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
765 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
766 CALL dbcsr_copy(ec_env%matrix_w(ispin, 1)%matrix, matrix_w(ispin, 1)%matrix)
767 END DO
768
769 END IF
770
771 ! External field (nonperiodic case)
772 ec_env%efield_nuclear = 0.0_dp
773 ec_env%efield_elec = 0.0_dp
774 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces=.false.)
775
776 CALL timestop(handle)
777
778 END SUBROUTINE ec_dc_energy
779
780! **************************************************************************************************
781!> \brief Kohn-Sham matrix contributions to force in DC-DFT
782!> also calculate right-hand-side matrix B for response equations AX=B
783!> \param qs_env ...
784!> \param ec_env ...
785!> \par History
786!> 08.2022 adapted from qs_ks_build_kohn_sham_matrix
787!> \author fbelle
788! **************************************************************************************************
789 SUBROUTINE ec_dc_build_ks_matrix_force(qs_env, ec_env)
790 TYPE(qs_environment_type), POINTER :: qs_env
791 TYPE(energy_correction_type), POINTER :: ec_env
792
793 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_dc_build_ks_matrix_force'
794
795 CHARACTER(LEN=default_string_length) :: unit_string
796 INTEGER :: handle, i, iounit, ispin, natom, nspins
797 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
798 use_virial
799 REAL(dp) :: dummy_real, dummy_real2(2), ehartree, &
800 eovrl, exc, fconv
801 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
802 REAL(dp), DIMENSION(3) :: fodeb, fodeb2
803 REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
804 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
805 TYPE(cell_type), POINTER :: cell
806 TYPE(cp_logger_type), POINTER :: logger
807 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, scrm
808 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
809 TYPE(dft_control_type), POINTER :: dft_control
810 TYPE(mp_para_env_type), POINTER :: para_env
811 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
812 POINTER :: sab_orb
813 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
814 TYPE(pw_env_type), POINTER :: pw_env
815 TYPE(pw_grid_type), POINTER :: pw_grid
816 TYPE(pw_poisson_type), POINTER :: poisson_env
817 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
818 TYPE(pw_r3d_rs_type) :: v_hartree_rspace
819 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace, v_rspace_in, &
820 v_tau_rspace
821 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
822 TYPE(qs_ks_env_type), POINTER :: ks_env
823 TYPE(qs_rho_type), POINTER :: rho
824 TYPE(section_vals_type), POINTER :: ec_hfx_sections
825 TYPE(virial_type), POINTER :: virial
826
827 CALL timeset(routinen, handle)
828
829 debug_forces = ec_env%debug_forces
830 debug_stress = ec_env%debug_stress
831
832 logger => cp_get_default_logger()
833 IF (logger%para_env%is_source()) THEN
834 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
835 ELSE
836 iounit = -1
837 END IF
838
839 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, matrix_ks, &
840 matrix_p, matrix_s, para_env, pw_env, rho, sab_orb, virial)
841 CALL get_qs_env(qs_env=qs_env, &
842 cell=cell, &
843 dft_control=dft_control, &
844 force=force, &
845 ks_env=ks_env, &
846 matrix_ks=matrix_ks, &
847 matrix_s=matrix_s, &
848 para_env=para_env, &
849 pw_env=pw_env, &
850 rho=rho, &
851 sab_orb=sab_orb, &
852 virial=virial)
853 cpassert(ASSOCIATED(pw_env))
854
855 nspins = dft_control%nspins
856 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
857
858 fconv = 1.0e-9_dp*pascal/cell%deth
859 IF (debug_stress .AND. use_virial) THEN
860 sttot = virial%pv_virial
861 END IF
862
863 ! Get density matrix of reference calculation
864 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
865
866 NULLIFY (auxbas_pw_pool, poisson_env)
867 ! gets the tmp grids
868 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
869 poisson_env=poisson_env)
870
871 ! Calculate the Hartree potential
872 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
873 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
874 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
875
876 ! Get the total input density in g-space [ions + electrons]
877 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
878
879 ! v_H[n_in]
880 IF (use_virial) THEN
881
882 ! Stress tensor - Volume and Green function contribution
883 h_stress(:, :) = 0.0_dp
884 CALL pw_poisson_solve(poisson_env, &
885 density=rho_tot_gspace, &
886 ehartree=ehartree, &
887 vhartree=v_hartree_gspace, &
888 h_stress=h_stress)
889
890 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
891 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
892
893 IF (debug_stress) THEN
894 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
895 CALL para_env%sum(stdeb)
896 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
897 'STRESS| GREEN 1st V_H[n_in]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
898 END IF
899
900 ELSE
901 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
902 v_hartree_gspace)
903 END IF
904
905 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
906 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
907
908 ! Save density on real space grid for use in properties
909 CALL qs_rho_get(rho, rho_r=rho_r)
910 ALLOCATE (ec_env%rhoout_r(nspins))
911 DO ispin = 1, nspins
912 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
913 CALL pw_copy(rho_r(ispin), ec_env%rhoout_r(ispin))
914 END DO
915
916 ! Getting nuclear force contribution from the core charge density
917 ! Vh(rho_c + rho_in)
918 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
919 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
920 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
921 IF (debug_forces) THEN
922 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
923 CALL para_env%sum(fodeb)
924 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb
925 END IF
926 IF (debug_stress .AND. use_virial) THEN
927 stdeb = fconv*(virial%pv_ehartree - stdeb)
928 CALL para_env%sum(stdeb)
929 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
930 'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb)
931 END IF
932
933 ! v_XC[n_in]_DC
934 ! v_rspace and v_tau_rspace are generated from the auxbas pool
935 NULLIFY (v_rspace, v_tau_rspace)
936
937 ! only activate stress calculation if
938 IF (use_virial) virial%pv_calculate = .true.
939
940 ! Exchange-correlation potential
941 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
942 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
943
944 IF (.NOT. ASSOCIATED(v_rspace)) THEN
945 ALLOCATE (v_rspace(nspins))
946 DO ispin = 1, nspins
947 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
948 CALL pw_zero(v_rspace(ispin))
949 END DO
950 END IF
951
952 IF (use_virial) THEN
953 virial%pv_exc = virial%pv_exc - virial%pv_xc
954 virial%pv_virial = virial%pv_virial - virial%pv_xc
955 ! virial%pv_xc will be zeroed in the xc routines
956 END IF
957
958 ! initialize srcm matrix
959 NULLIFY (scrm)
960 CALL dbcsr_allocate_matrix_set(scrm, nspins)
961 DO ispin = 1, nspins
962 ALLOCATE (scrm(ispin)%matrix)
963 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
964 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
965 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
966 END DO
967
968 pw_grid => v_hartree_rspace%pw_grid
969 ALLOCATE (v_rspace_in(nspins))
970 DO ispin = 1, nspins
971 CALL v_rspace_in(ispin)%create(pw_grid)
972 END DO
973
974 ! v_rspace_in = v_H[n_in] + v_xc[n_in] calculated in ks_ref_potential
975 DO ispin = 1, nspins
976 ! v_xc[n_in]_GS
977 CALL pw_transfer(ec_env%vxc_rspace(ispin), v_rspace_in(ispin))
978 ! add v_H[n_in]
979 CALL pw_axpy(ec_env%vh_rspace, v_rspace_in(ispin))
980 END DO
981
982!------------------------------------------------
983
984 ! If hybrid functional in DC-DFT
985 ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF")
986 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
987
988 IF (do_ec_hfx) THEN
989
990 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
991 IF (debug_forces) fodeb2(1:3) = force(1)%overlap_admm(1:3, 1)
992
993 ! Calculate direct HFX forces here
994 ! Virial contribution (fock_4c) done inside calculate_exx
995 dummy_real = 0.0_dp
996 CALL calculate_exx(qs_env=qs_env, &
997 unit_nr=iounit, &
998 hfx_sections=ec_hfx_sections, &
999 x_data=ec_env%x_data, &
1000 do_gw=.false., &
1001 do_admm=ec_env%do_ec_admm, &
1002 calc_forces=.true., &
1003 reuse_hfx=ec_env%reuse_hfx, &
1004 do_im_time=.false., &
1005 e_ex_from_gw=dummy_real, &
1006 e_admm_from_gw=dummy_real2, &
1007 t3=dummy_real)
1008
1009 IF (debug_forces) THEN
1010 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
1011 CALL para_env%sum(fodeb)
1012 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC ", fodeb
1013
1014 fodeb2(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb2(1:3)
1015 CALL para_env%sum(fodeb2)
1016 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx_DC*S ", fodeb2
1017 END IF
1018 IF (debug_stress .AND. use_virial) THEN
1019 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
1020 CALL para_env%sum(stdeb)
1021 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1022 'STRESS| P*hfx_DC ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1023 END IF
1024
1025 END IF
1026
1027!------------------------------------------------
1028
1029 ! Stress-tensor contribution derivative of integrand
1030 ! int v_Hxc[n^în]*n^out
1031 IF (use_virial) THEN
1032 pv_loc = virial%pv_virial
1033 END IF
1034
1035 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1036 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1037
1038 DO ispin = 1, nspins
1039 ! Add v_H[n_in] + v_xc[n_in] = v_rspace
1040 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1041 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
1042 ! integrate over potential <a|V|b>
1043 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1044 hmat=scrm(ispin), &
1045 pmat=matrix_p(ispin, 1), &
1046 qs_env=qs_env, &
1047 calculate_forces=.true., &
1048 basis_type="HARRIS", &
1049 task_list_external=ec_env%task_list)
1050 END DO
1051
1052 IF (debug_forces) THEN
1053 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1054 CALL para_env%sum(fodeb)
1055 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb
1056 END IF
1057 IF (debug_stress .AND. use_virial) THEN
1058 stdeb = fconv*(virial%pv_virial - stdeb)
1059 CALL para_env%sum(stdeb)
1060 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1061 'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1062 END IF
1063
1064 IF (ASSOCIATED(v_tau_rspace)) THEN
1065 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1066 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1067 DO ispin = 1, nspins
1068 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1069 ! integrate over Tau-potential <nabla.a|V|nabla.b>
1070 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1071 hmat=scrm(ispin), &
1072 pmat=matrix_p(ispin, 1), &
1073 qs_env=qs_env, &
1074 calculate_forces=.true., &
1075 compute_tau=.true., &
1076 basis_type="HARRIS", &
1077 task_list_external=ec_env%task_list)
1078 END DO
1079
1080 IF (debug_forces) THEN
1081 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1082 CALL para_env%sum(fodeb)
1083 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb
1084 END IF
1085 IF (debug_stress .AND. use_virial) THEN
1086 stdeb = fconv*(virial%pv_virial - stdeb)
1087 CALL para_env%sum(stdeb)
1088 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1089 'STRESS| INT Pout*dVhxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1090 END IF
1091 END IF
1092
1093 ! Stress-tensor
1094 IF (use_virial) THEN
1095 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1096 END IF
1097
1098 ! delete scrm matrix
1100
1101 !----------------------------------------------------
1102 ! Right-hand-side matrix B for linear response equations AX = B
1103 !----------------------------------------------------
1104
1105 ! RHS = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC * E_X[n] - alpha_gs * E_X[n]
1106 ! = int v_Hxc[n]_DC - v_Hxc[n]_GS dr + alpha_DC / alpha_GS * E_X[n]_GS - E_X[n]_GS
1107 !
1108 ! with v_Hxc[n] = v_H[n] + v_xc[n]
1109 !
1110 ! Actually v_H[n_in] same for DC and GS, just there for convenience
1111 ! v_xc[n_in]_GS = 0 if GS is HF BUT =/0 if hybrid
1112 ! so, we keep this general form
1113
1114 NULLIFY (ec_env%matrix_hz)
1115 CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
1116 DO ispin = 1, nspins
1117 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
1118 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
1119 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
1120 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
1121 END DO
1122
1123 DO ispin = 1, nspins
1124 ! v_rspace = v_rspace - v_rspace_in
1125 ! = v_Hxc[n_in]_DC - v_Hxc[n_in]_GS
1126 CALL pw_axpy(v_rspace_in(ispin), v_rspace(ispin), -1.0_dp)
1127 END DO
1128
1129 DO ispin = 1, nspins
1130 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1131 hmat=ec_env%matrix_hz(ispin), &
1132 pmat=matrix_p(ispin, 1), &
1133 qs_env=qs_env, &
1134 calculate_forces=.false., &
1135 basis_type="HARRIS", &
1136 task_list_external=ec_env%task_list)
1137 END DO
1138
1139 ! Check if mGGA functionals are used
1140 IF (dft_control%use_kinetic_energy_density) THEN
1141
1142 ! If DC-DFT without mGGA functional, this needs to be allocated now.
1143 IF (.NOT. ASSOCIATED(v_tau_rspace)) THEN
1144 ALLOCATE (v_tau_rspace(nspins))
1145 DO ispin = 1, nspins
1146 CALL auxbas_pw_pool%create_pw(v_tau_rspace(ispin))
1147 CALL pw_zero(v_tau_rspace(ispin))
1148 END DO
1149 END IF
1150
1151 DO ispin = 1, nspins
1152 ! v_tau_rspace = v_Hxc_tau[n_in]_DC - v_Hxc_tau[n_in]_GS
1153 IF (ASSOCIATED(ec_env%vtau_rspace)) THEN
1154 CALL pw_axpy(ec_env%vtau_rspace(ispin), v_tau_rspace(ispin), -1.0_dp)
1155 END IF
1156 ! integrate over Tau-potential <nabla.a|V|nabla.b>
1157 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1158 hmat=ec_env%matrix_hz(ispin), &
1159 pmat=matrix_p(ispin, 1), &
1160 qs_env=qs_env, &
1161 calculate_forces=.false., compute_tau=.true., &
1162 basis_type="HARRIS", &
1163 task_list_external=ec_env%task_list)
1164 END DO
1165 END IF
1166
1167 ! Need to also subtract HFX contribution of reference calculation from ec_env%matrix_hz
1168 ! and/or add HFX contribution if DC-DFT ueses hybrid XC-functional
1169 CALL add_exx_to_rhs(rhs=ec_env%matrix_hz, &
1170 qs_env=qs_env, &
1171 ext_hfx_section=ec_hfx_sections, &
1172 x_data=ec_env%x_data, &
1173 recalc_integrals=.false., &
1174 do_admm=ec_env%do_ec_admm, &
1175 do_ec=.true., &
1176 do_exx=.false., &
1177 reuse_hfx=ec_env%reuse_hfx)
1178
1179 ! Core overlap
1180 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
1181 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
1182 CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
1183 IF (debug_forces) THEN
1184 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
1185 CALL para_env%sum(fodeb)
1186 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb
1187 END IF
1188 IF (debug_stress .AND. use_virial) THEN
1189 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
1190 CALL para_env%sum(stdeb)
1191 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1192 'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1193 END IF
1194
1195 IF (debug_forces) THEN
1196 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
1197 ALLOCATE (ftot(3, natom))
1198 CALL total_qs_force(ftot, force, atomic_kind_set)
1199 fodeb(1:3) = ftot(1:3, 1)
1200 DEALLOCATE (ftot)
1201 CALL para_env%sum(fodeb)
1202 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Force Explicit", fodeb
1203 END IF
1204
1205 ! return pw grids
1206 DO ispin = 1, nspins
1207 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1208 CALL auxbas_pw_pool%give_back_pw(v_rspace_in(ispin))
1209 IF (ASSOCIATED(v_tau_rspace)) THEN
1210 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1211 END IF
1212 END DO
1213
1214 DEALLOCATE (v_rspace, v_rspace_in)
1215 IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
1216 !
1217 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1218 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1219 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1220
1221 ! Stress tensor - volume terms need to be stored,
1222 ! for a sign correction in QS at the end of qs_force
1223 IF (use_virial) THEN
1224 IF (qs_env%energy_correction) THEN
1225 ec_env%ehartree = ehartree
1226 ec_env%exc = exc
1227 END IF
1228 END IF
1229
1230 IF (debug_stress .AND. use_virial) THEN
1231 ! In total: -1.0*E_H
1232 stdeb = -1.0_dp*fconv*ehartree
1233 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1234 'STRESS| VOL 1st v_H[n_in]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb)
1235
1236 stdeb = -1.0_dp*fconv*exc
1237 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1238 'STRESS| VOL 1st E_XC_DC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb)
1239
1240 ! For debugging, create a second virial environment,
1241 ! apply volume terms immediately
1242 block
1243 TYPE(virial_type) :: virdeb
1244 virdeb = virial
1245
1246 CALL para_env%sum(virdeb%pv_overlap)
1247 CALL para_env%sum(virdeb%pv_ekinetic)
1248 CALL para_env%sum(virdeb%pv_ppl)
1249 CALL para_env%sum(virdeb%pv_ppnl)
1250 CALL para_env%sum(virdeb%pv_ecore_overlap)
1251 CALL para_env%sum(virdeb%pv_ehartree)
1252 CALL para_env%sum(virdeb%pv_exc)
1253 CALL para_env%sum(virdeb%pv_exx)
1254 CALL para_env%sum(virdeb%pv_vdw)
1255 CALL para_env%sum(virdeb%pv_mp2)
1256 CALL para_env%sum(virdeb%pv_nlcc)
1257 CALL para_env%sum(virdeb%pv_gapw)
1258 CALL para_env%sum(virdeb%pv_lrigpw)
1259 CALL para_env%sum(virdeb%pv_virial)
1260 CALL symmetrize_virial(virdeb)
1261
1262 ! apply stress-tensor 1st terms
1263 DO i = 1, 3
1264 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*ehartree
1265 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc &
1266 - 2.0_dp*ehartree
1267 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc
1268 ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
1269 ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
1270 ! There should be a more elegant solution to that ...
1271 END DO
1272
1273 CALL para_env%sum(sttot)
1274 stdeb = fconv*(virdeb%pv_virial - sttot)
1275 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1276 'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1277
1278 stdeb = fconv*(virdeb%pv_virial)
1279 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1280 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1281
1282 unit_string = "GPa" ! old default
1283 CALL write_stress_tensor_components(virdeb, iounit, cell, unit_string)
1284 CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, unit_string, .false.)
1285
1286 END block
1287 END IF
1288
1289 CALL timestop(handle)
1290
1291 END SUBROUTINE ec_dc_build_ks_matrix_force
1292
1293! **************************************************************************************************
1294!> \brief ...
1295!> \param qs_env ...
1296!> \param ec_env ...
1297!> \param calculate_forces ...
1298! **************************************************************************************************
1299 SUBROUTINE ec_disp(qs_env, ec_env, calculate_forces)
1300 TYPE(qs_environment_type), POINTER :: qs_env
1301 TYPE(energy_correction_type), POINTER :: ec_env
1302 LOGICAL, INTENT(IN) :: calculate_forces
1303
1304 REAL(kind=dp) :: edisp
1305
1306 CALL calculate_dispersion_pairpot(qs_env, ec_env%dispersion_env, edisp, calculate_forces)
1307 IF (.NOT. calculate_forces) ec_env%edispersion = ec_env%edispersion + edisp
1308
1309 END SUBROUTINE ec_disp
1310
1311! **************************************************************************************************
1312!> \brief Construction of the Core Hamiltonian Matrix
1313!> Short version of qs_core_hamiltonian
1314!> \param qs_env ...
1315!> \param ec_env ...
1316!> \author Creation (03.2014,JGH)
1317! **************************************************************************************************
1318 SUBROUTINE ec_build_core_hamiltonian(qs_env, ec_env)
1319 TYPE(qs_environment_type), POINTER :: qs_env
1320 TYPE(energy_correction_type), POINTER :: ec_env
1321
1322 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_core_hamiltonian'
1323
1324 CHARACTER(LEN=default_string_length) :: basis_type
1325 INTEGER :: handle, nder, nimages
1326 LOGICAL :: calculate_forces, use_virial
1327 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1328 TYPE(dft_control_type), POINTER :: dft_control
1329 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1330 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1331 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1332 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1333 TYPE(qs_ks_env_type), POINTER :: ks_env
1334
1335 CALL timeset(routinen, handle)
1336
1337 NULLIFY (atomic_kind_set, dft_control, ks_env, particle_set, &
1338 qs_kind_set)
1339
1340 CALL get_qs_env(qs_env=qs_env, &
1341 atomic_kind_set=atomic_kind_set, &
1342 dft_control=dft_control, &
1343 particle_set=particle_set, &
1344 qs_kind_set=qs_kind_set, &
1345 ks_env=ks_env)
1346
1347 ! no k-points possible
1348 nimages = dft_control%nimages
1349 IF (nimages /= 1) THEN
1350 cpabort("K-points for Harris functional not implemented")
1351 END IF
1352
1353 ! check for GAPW/GAPW_XC
1354 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1355 cpabort("Harris functional for GAPW not implemented")
1356 END IF
1357
1358 ! Do not calculate forces or stress tensor here
1359 use_virial = .false.
1360 calculate_forces = .false.
1361
1362 ! get neighbor lists, we need the full sab_orb list from the ec_env
1363 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
1364 sab_orb => ec_env%sab_orb
1365 sac_ae => ec_env%sac_ae
1366 sac_ppl => ec_env%sac_ppl
1367 sap_ppnl => ec_env%sap_ppnl
1368
1369 basis_type = "HARRIS"
1370
1371 nder = 0
1372 ! Overlap and kinetic energy matrices
1373 CALL build_overlap_matrix(ks_env, matrixkp_s=ec_env%matrix_s, &
1374 matrix_name="OVERLAP MATRIX", &
1375 basis_type_a=basis_type, &
1376 basis_type_b=basis_type, &
1377 sab_nl=sab_orb)
1378 CALL build_kinetic_matrix(ks_env, matrixkp_t=ec_env%matrix_t, &
1379 matrix_name="KINETIC ENERGY MATRIX", &
1380 basis_type=basis_type, &
1381 sab_nl=sab_orb)
1382
1383 ! initialize H matrix
1384 CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
1385 ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
1386 CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, template=ec_env%matrix_s(1, 1)%matrix)
1387 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
1388
1389 ! add kinetic energy
1390 CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, ec_env%matrix_t(1, 1)%matrix, &
1391 keep_sparsity=.true., name="CORE HAMILTONIAN MATRIX")
1392
1393 CALL core_matrices(qs_env, ec_env%matrix_h, ec_env%matrix_p, calculate_forces, nder, &
1394 ec_env=ec_env, ec_env_matrices=.true., basis_type=basis_type)
1395
1396 ! External field (nonperiodic case)
1397 ec_env%efield_nuclear = 0.0_dp
1398 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1399
1400 CALL timestop(handle)
1401
1402 END SUBROUTINE ec_build_core_hamiltonian
1403
1404! **************************************************************************************************
1405!> \brief Solve KS equation for a given matrix
1406!> calculate the complete KS matrix
1407!> \param qs_env ...
1408!> \param ec_env ...
1409!> \par History
1410!> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH]
1411!> \author JGH
1412! **************************************************************************************************
1413 SUBROUTINE ec_build_ks_matrix(qs_env, ec_env)
1414 TYPE(qs_environment_type), POINTER :: qs_env
1415 TYPE(energy_correction_type), POINTER :: ec_env
1416
1417 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_ks_matrix'
1418
1419 CHARACTER(LEN=default_string_length) :: headline
1420 INTEGER :: handle, iounit, ispin, nspins
1421 LOGICAL :: calculate_forces, &
1422 do_adiabatic_rescaling, do_ec_hfx, &
1423 hfx_treat_lsd_in_core, use_virial
1424 REAL(dp) :: dummy_real, dummy_real2(2), eexc, evhxc, &
1425 t3
1426 TYPE(cp_logger_type), POINTER :: logger
1427 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_mat
1428 TYPE(dft_control_type), POINTER :: dft_control
1429 TYPE(pw_env_type), POINTER :: pw_env
1430 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1431 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r, v_rspace, v_tau_rspace
1432 TYPE(qs_energy_type), POINTER :: energy
1433 TYPE(qs_ks_env_type), POINTER :: ks_env
1434 TYPE(qs_rho_type), POINTER :: rho
1435 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
1436 ec_hfx_sections, ec_section
1437
1438 CALL timeset(routinen, handle)
1439
1440 logger => cp_get_default_logger()
1441 IF (logger%para_env%is_source()) THEN
1442 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1443 ELSE
1444 iounit = -1
1445 END IF
1446
1447 ! get all information on the electronic density
1448 NULLIFY (auxbas_pw_pool, dft_control, energy, ks_env, rho, rho_r, tau_r)
1449 CALL get_qs_env(qs_env=qs_env, &
1450 dft_control=dft_control, &
1451 ks_env=ks_env, &
1452 rho=rho)
1453 nspins = dft_control%nspins
1454 calculate_forces = .false.
1455 use_virial = .false.
1456
1457 ! Kohn-Sham matrix
1458 IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
1459 CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
1460 DO ispin = 1, nspins
1461 headline = "KOHN-SHAM MATRIX"
1462 ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
1463 CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=trim(headline), &
1464 template=ec_env%matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
1465 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%sab_orb)
1466 CALL dbcsr_set(ec_env%matrix_ks(ispin, 1)%matrix, 0.0_dp)
1467 END DO
1468
1469 NULLIFY (pw_env)
1470 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1471 cpassert(ASSOCIATED(pw_env))
1472
1473 ! Exact exchange contribution (hybrid functionals)
1474 ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
1475 ec_hfx_sections => section_vals_get_subs_vals(ec_section, "XC%HF")
1476 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
1477
1478 IF (do_ec_hfx) THEN
1479
1480 ! Check what works
1481 adiabatic_rescaling_section => section_vals_get_subs_vals(ec_section, "XC%ADIABATIC_RESCALING")
1482 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1483 IF (do_adiabatic_rescaling) THEN
1484 CALL cp_abort(__location__, "Adiabatic rescaling NYI for energy correction")
1485 END IF
1486 CALL section_vals_val_get(ec_hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core)
1487 IF (hfx_treat_lsd_in_core) THEN
1488 CALL cp_abort(__location__, "HFX_TREAT_LSD_IN_CORE NYI for energy correction")
1489 END IF
1490
1491 ! calculate the density matrix for the fitted mo_coeffs
1492 IF (dft_control%do_admm) THEN
1493
1494 IF (dft_control%do_admm_mo) THEN
1495 IF (qs_env%run_rtp) THEN
1496 CALL rtp_admm_calc_rho_aux(qs_env)
1497 ELSE
1498 CALL admm_mo_calc_rho_aux(qs_env)
1499 END IF
1500 ELSEIF (dft_control%do_admm_dm) THEN
1501 CALL admm_dm_calc_rho_aux(qs_env)
1502 END IF
1503 END IF
1504
1505 ! Get exact exchange energy
1506 dummy_real = 0.0_dp
1507 t3 = 0.0_dp
1508 CALL get_qs_env(qs_env, energy=energy)
1509 CALL calculate_exx(qs_env=qs_env, &
1510 unit_nr=iounit, &
1511 hfx_sections=ec_hfx_sections, &
1512 x_data=ec_env%x_data, &
1513 do_gw=.false., &
1514 do_admm=ec_env%do_ec_admm, &
1515 calc_forces=.false., &
1516 reuse_hfx=ec_env%reuse_hfx, &
1517 do_im_time=.false., &
1518 e_ex_from_gw=dummy_real, &
1519 e_admm_from_gw=dummy_real2, &
1520 t3=dummy_real)
1521
1522 ! Save exchange energy
1523 ec_env%ex = energy%ex
1524 ! Save EXX ADMM XC correction
1525 IF (ec_env%do_ec_admm) THEN
1526 ec_env%exc_aux_fit = energy%exc_aux_fit + energy%exc
1527 END IF
1528
1529 ! Add exact echange contribution of EC to EC Hamiltonian
1530 ! do_ec = .FALSE prevents subtraction of HFX contribution of reference calculation
1531 ! do_exx = .FALSE. prevents subtraction of reference XC contribution
1532 ks_mat => ec_env%matrix_ks(:, 1)
1533 CALL add_exx_to_rhs(rhs=ks_mat, &
1534 qs_env=qs_env, &
1535 ext_hfx_section=ec_hfx_sections, &
1536 x_data=ec_env%x_data, &
1537 recalc_integrals=.false., &
1538 do_admm=ec_env%do_ec_admm, &
1539 do_ec=.false., &
1540 do_exx=.false., &
1541 reuse_hfx=ec_env%reuse_hfx)
1542
1543 END IF
1544
1545 ! v_rspace and v_tau_rspace are generated from the auxbas pool
1546 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1547 NULLIFY (v_rspace, v_tau_rspace)
1548 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1549 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
1550
1551 IF (.NOT. ASSOCIATED(v_rspace)) THEN
1552 ALLOCATE (v_rspace(nspins))
1553 DO ispin = 1, nspins
1554 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
1555 CALL pw_zero(v_rspace(ispin))
1556 END DO
1557 END IF
1558
1559 evhxc = 0.0_dp
1560 CALL qs_rho_get(rho, rho_r=rho_r)
1561 IF (ASSOCIATED(v_tau_rspace)) THEN
1562 CALL qs_rho_get(rho, tau_r=tau_r)
1563 END IF
1564 DO ispin = 1, nspins
1565 ! Add v_hartree + v_xc = v_rspace
1566 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
1567 CALL pw_axpy(ec_env%vh_rspace, v_rspace(ispin))
1568 ! integrate over potential <a|V|b>
1569 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1570 hmat=ec_env%matrix_ks(ispin, 1), &
1571 qs_env=qs_env, &
1572 calculate_forces=.false., &
1573 basis_type="HARRIS", &
1574 task_list_external=ec_env%task_list)
1575
1576 IF (ASSOCIATED(v_tau_rspace)) THEN
1577 ! integrate over Tau-potential <nabla.a|V|nabla.b>
1578 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1579 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1580 hmat=ec_env%matrix_ks(ispin, 1), &
1581 qs_env=qs_env, &
1582 calculate_forces=.false., &
1583 compute_tau=.true., &
1584 basis_type="HARRIS", &
1585 task_list_external=ec_env%task_list)
1586 END IF
1587
1588 ! calclulate Int(vhxc*rho)dr and Int(vtau*tau)dr
1589 evhxc = evhxc + pw_integral_ab(rho_r(ispin), v_rspace(ispin))/ &
1590 v_rspace(1)%pw_grid%dvol
1591 IF (ASSOCIATED(v_tau_rspace)) THEN
1592 evhxc = evhxc + pw_integral_ab(tau_r(ispin), v_tau_rspace(ispin))/ &
1593 v_tau_rspace(ispin)%pw_grid%dvol
1594 END IF
1595
1596 END DO
1597
1598 ! return pw grids
1599 DO ispin = 1, nspins
1600 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1601 IF (ASSOCIATED(v_tau_rspace)) THEN
1602 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1603 END IF
1604 END DO
1605 DEALLOCATE (v_rspace)
1606 IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
1607
1608 ! energies
1609 ec_env%exc = eexc
1610 ec_env%vhxc = evhxc
1611
1612 ! add the core matrix
1613 DO ispin = 1, nspins
1614 CALL dbcsr_add(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_h(1, 1)%matrix, &
1615 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1616 CALL dbcsr_filter(ec_env%matrix_ks(ispin, 1)%matrix, &
1617 dft_control%qs_control%eps_filter_matrix)
1618 END DO
1619
1620 CALL timestop(handle)
1621
1622 END SUBROUTINE ec_build_ks_matrix
1623
1624! **************************************************************************************************
1625!> \brief Construction of the Core Hamiltonian Matrix
1626!> Short version of qs_core_hamiltonian
1627!> \param qs_env ...
1628!> \param ec_env ...
1629!> \param matrix_p ...
1630!> \param matrix_s ...
1631!> \param matrix_w ...
1632!> \author Creation (03.2014,JGH)
1633! **************************************************************************************************
1634 SUBROUTINE ec_build_core_hamiltonian_force(qs_env, ec_env, matrix_p, matrix_s, matrix_w)
1635 TYPE(qs_environment_type), POINTER :: qs_env
1636 TYPE(energy_correction_type), POINTER :: ec_env
1637 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s, matrix_w
1638
1639 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_core_hamiltonian_force'
1640
1641 CHARACTER(LEN=default_string_length) :: basis_type
1642 INTEGER :: handle, iounit, nder, nimages
1643 LOGICAL :: calculate_forces, debug_forces, &
1644 debug_stress, use_virial
1645 REAL(kind=dp) :: fconv
1646 REAL(kind=dp), DIMENSION(3) :: fodeb
1647 REAL(kind=dp), DIMENSION(3, 3) :: stdeb, sttot
1648 TYPE(cell_type), POINTER :: cell
1649 TYPE(cp_logger_type), POINTER :: logger
1650 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
1651 TYPE(dft_control_type), POINTER :: dft_control
1652 TYPE(mp_para_env_type), POINTER :: para_env
1653 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1654 POINTER :: sab_orb
1655 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1656 TYPE(qs_ks_env_type), POINTER :: ks_env
1657 TYPE(virial_type), POINTER :: virial
1658
1659 CALL timeset(routinen, handle)
1660
1661 debug_forces = ec_env%debug_forces
1662 debug_stress = ec_env%debug_stress
1663
1664 logger => cp_get_default_logger()
1665 IF (logger%para_env%is_source()) THEN
1666 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1667 ELSE
1668 iounit = -1
1669 END IF
1670
1671 calculate_forces = .true.
1672
1673 basis_type = "HARRIS"
1674
1675 ! no k-points possible
1676 NULLIFY (cell, dft_control, force, ks_env, para_env, virial)
1677 CALL get_qs_env(qs_env=qs_env, &
1678 cell=cell, &
1679 dft_control=dft_control, &
1680 force=force, &
1681 ks_env=ks_env, &
1682 para_env=para_env, &
1683 virial=virial)
1684 nimages = dft_control%nimages
1685 IF (nimages /= 1) THEN
1686 cpabort("K-points for Harris functional not implemented")
1687 END IF
1688
1689 ! check for virial
1690 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1691
1692 fconv = 1.0e-9_dp*pascal/cell%deth
1693 IF (debug_stress .AND. use_virial) THEN
1694 sttot = virial%pv_virial
1695 END IF
1696
1697 ! check for GAPW/GAPW_XC
1698 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1699 IF (ec_env%energy_functional == ec_functional_harris) THEN
1700 cpabort("Harris functional for GAPW not implemented")
1701 END IF
1702 END IF
1703
1704 ! get neighbor lists, we need the full sab_orb list from the ec_env
1705 sab_orb => ec_env%sab_orb
1706
1707 ! initialize src matrix
1708 NULLIFY (scrm)
1709 CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
1710 ALLOCATE (scrm(1, 1)%matrix)
1711 CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_s(1, 1)%matrix)
1712 CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
1713
1714 nder = 1
1715 IF (SIZE(matrix_p, 1) == 2) THEN
1716 CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
1717 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1718 END IF
1719
1720 ! Overlap and kinetic energy matrices
1721 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
1722 IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
1723 CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
1724 matrix_name="OVERLAP MATRIX", &
1725 basis_type_a=basis_type, &
1726 basis_type_b=basis_type, &
1727 sab_nl=sab_orb, calculate_forces=.true., &
1728 matrixkp_p=matrix_w)
1729
1730 IF (debug_forces) THEN
1731 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
1732 CALL para_env%sum(fodeb)
1733 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
1734 END IF
1735 IF (debug_stress .AND. use_virial) THEN
1736 stdeb = fconv*(virial%pv_overlap - stdeb)
1737 CALL para_env%sum(stdeb)
1738 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1739 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
1740 END IF
1741
1742 CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, &
1743 calculate_forces=.true., &
1744 sab_orb=sab_orb, &
1745 basis_type=basis_type, &
1746 debug_forces=debug_forces, debug_stress=debug_stress)
1747
1748 CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, &
1749 ec_env=ec_env, ec_env_matrices=.false., basis_type=basis_type, &
1750 debug_forces=debug_forces, debug_stress=debug_stress)
1751
1752 ! External field (nonperiodic case)
1753 ec_env%efield_nuclear = 0.0_dp
1754 IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%efield(1:3, 1)
1755 CALL ec_efield_local_operator(qs_env, ec_env, calculate_forces)
1756 IF (calculate_forces .AND. debug_forces) THEN
1757 fodeb(1:3) = force(1)%efield(1:3, 1) - fodeb(1:3)
1758 CALL para_env%sum(fodeb)
1759 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dEfield", fodeb
1760 END IF
1761 IF (debug_stress .AND. use_virial) THEN
1762 stdeb = fconv*(virial%pv_virial - sttot)
1763 CALL para_env%sum(stdeb)
1764 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1765 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1766 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") ' '
1767 END IF
1768
1769 ! delete scr matrix
1770 CALL dbcsr_deallocate_matrix_set(scrm)
1771
1772 CALL timestop(handle)
1773
1774 END SUBROUTINE ec_build_core_hamiltonian_force
1775
1776! **************************************************************************************************
1777!> \brief Solve KS equation for a given matrix
1778!> \brief calculate the complete KS matrix
1779!> \param qs_env ...
1780!> \param ec_env ...
1781!> \par History
1782!> 03.2014 adapted from qs_ks_build_kohn_sham_matrix [JGH]
1783!> \author JGH
1784! **************************************************************************************************
1785 SUBROUTINE ec_build_ks_matrix_force(qs_env, ec_env)
1786 TYPE(qs_environment_type), POINTER :: qs_env
1787 TYPE(energy_correction_type), POINTER :: ec_env
1788
1789 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_ks_matrix_force'
1790
1791 CHARACTER(LEN=default_string_length) :: unit_string
1792 INTEGER :: handle, i, iounit, ispin, natom, nspins
1793 LOGICAL :: debug_forces, debug_stress, do_ec_hfx, &
1794 use_virial
1795 REAL(dp) :: dehartree, dummy_real, dummy_real2(2), &
1796 eexc, ehartree, eovrl, exc, fconv
1797 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ftot
1798 REAL(dp), DIMENSION(3) :: fodeb
1799 REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot
1800 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1801 TYPE(cell_type), POINTER :: cell
1802 TYPE(cp_logger_type), POINTER :: logger
1803 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, scrm
1804 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1805 TYPE(dft_control_type), POINTER :: dft_control
1806 TYPE(mp_para_env_type), POINTER :: para_env
1807 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1808 POINTER :: sab_orb
1809 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhodn_tot_gspace, &
1810 v_hartree_gspace
1811 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rhoout_g
1812 TYPE(pw_c1d_gs_type), POINTER :: rho_core
1813 TYPE(pw_env_type), POINTER :: pw_env
1814 TYPE(pw_poisson_type), POINTER :: poisson_env
1815 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1816 TYPE(pw_r3d_rs_type) :: dv_hartree_rspace, v_hartree_rspace, &
1817 vtot_rspace
1818 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rhoout_r, tau_r, tauout_r, &
1819 v_rspace, v_tau_rspace, v_xc, v_xc_tau
1820 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1821 TYPE(qs_ks_env_type), POINTER :: ks_env
1822 TYPE(qs_rho_type), POINTER :: rho
1823 TYPE(section_vals_type), POINTER :: ec_hfx_sections, xc_section
1824 TYPE(virial_type), POINTER :: virial
1825
1826 CALL timeset(routinen, handle)
1827
1828 debug_forces = ec_env%debug_forces
1829 debug_stress = ec_env%debug_stress
1830
1831 logger => cp_get_default_logger()
1832 IF (logger%para_env%is_source()) THEN
1833 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1834 ELSE
1835 iounit = -1
1836 END IF
1837
1838 ! get all information on the electronic density
1839 NULLIFY (atomic_kind_set, cell, dft_control, force, ks_env, &
1840 matrix_ks, matrix_p, matrix_s, para_env, rho, rho_core, &
1841 rho_g, rho_r, sab_orb, tau_r, virial)
1842 CALL get_qs_env(qs_env=qs_env, &
1843 cell=cell, &
1844 dft_control=dft_control, &
1845 force=force, &
1846 ks_env=ks_env, &
1847 matrix_ks=matrix_ks, &
1848 para_env=para_env, &
1849 rho=rho, &
1850 sab_orb=sab_orb, &
1851 virial=virial)
1852
1853 nspins = dft_control%nspins
1854 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1855
1856 ! Conversion factor a.u. -> GPa
1857 unit_string = "GPa"
1858 fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, trim(unit_string))
1859
1860 IF (debug_stress .AND. use_virial) THEN
1861 sttot = virial%pv_virial
1862 END IF
1863
1864 NULLIFY (pw_env)
1865 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1866 cpassert(ASSOCIATED(pw_env))
1867
1868 NULLIFY (auxbas_pw_pool, poisson_env)
1869 ! gets the tmp grids
1870 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1871 poisson_env=poisson_env)
1872
1873 ! Calculate the Hartree potential
1874 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
1875 CALL auxbas_pw_pool%create_pw(rhodn_tot_gspace)
1876 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
1877
1878 CALL pw_transfer(ec_env%vh_rspace, v_hartree_rspace)
1879
1880 ! calculate output density on grid
1881 ! rho_in(R): CALL qs_rho_get(rho, rho_r=rho_r)
1882 ! rho_in(G): CALL qs_rho_get(rho, rho_g=rho_g)
1883 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
1884 NULLIFY (rhoout_r, rhoout_g)
1885 ALLOCATE (rhoout_r(nspins), rhoout_g(nspins))
1886 DO ispin = 1, nspins
1887 CALL auxbas_pw_pool%create_pw(rhoout_r(ispin))
1888 CALL auxbas_pw_pool%create_pw(rhoout_g(ispin))
1889 END DO
1890 CALL auxbas_pw_pool%create_pw(dv_hartree_rspace)
1891 CALL auxbas_pw_pool%create_pw(vtot_rspace)
1892
1893 CALL pw_zero(rhodn_tot_gspace)
1894 DO ispin = 1, nspins
1895 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1896 rho=rhoout_r(ispin), &
1897 rho_gspace=rhoout_g(ispin), &
1898 basis_type="HARRIS", &
1899 task_list_external=ec_env%task_list)
1900 END DO
1901
1902 ! Save Harris on real space grid for use in properties
1903 ALLOCATE (ec_env%rhoout_r(nspins))
1904 DO ispin = 1, nspins
1905 CALL auxbas_pw_pool%create_pw(ec_env%rhoout_r(ispin))
1906 CALL pw_copy(rhoout_r(ispin), ec_env%rhoout_r(ispin))
1907 END DO
1908
1909 NULLIFY (tauout_r)
1910 IF (dft_control%use_kinetic_energy_density) THEN
1911 block
1912 TYPE(pw_c1d_gs_type) :: tauout_g
1913 ALLOCATE (tauout_r(nspins))
1914 DO ispin = 1, nspins
1915 CALL auxbas_pw_pool%create_pw(tauout_r(ispin))
1916 END DO
1917 CALL auxbas_pw_pool%create_pw(tauout_g)
1918
1919 DO ispin = 1, nspins
1920 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=ec_env%matrix_p(ispin, 1)%matrix, &
1921 rho=tauout_r(ispin), &
1922 rho_gspace=tauout_g, &
1923 compute_tau=.true., &
1924 basis_type="HARRIS", &
1925 task_list_external=ec_env%task_list)
1926 END DO
1927
1928 CALL auxbas_pw_pool%give_back_pw(tauout_g)
1929 END block
1930 END IF
1931
1932 IF (use_virial) THEN
1933
1934 ! Calculate the Hartree potential
1935 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1936
1937 ! Get the total input density in g-space [ions + electrons]
1938 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1939
1940 ! make rho_tot_gspace with output density
1941 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1942 CALL pw_copy(rho_core, rhodn_tot_gspace)
1943 DO ispin = 1, dft_control%nspins
1944 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
1945 END DO
1946
1947 ! Volume and Green function terms
1948 h_stress(:, :) = 0.0_dp
1949 CALL pw_poisson_solve(poisson_env, &
1950 density=rho_tot_gspace, & ! n_in
1951 ehartree=ehartree, &
1952 vhartree=v_hartree_gspace, & ! v_H[n_in]
1953 h_stress=h_stress, &
1954 aux_density=rhodn_tot_gspace) ! n_out
1955
1956 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1957 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1958
1959 IF (debug_stress) THEN
1960 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1961 CALL para_env%sum(stdeb)
1962 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1963 'STRESS| GREEN 1st v_H[n_in]*n_out ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1964 END IF
1965
1966 ! activate stress calculation
1967 virial%pv_calculate = .true.
1968
1969 NULLIFY (v_rspace, v_tau_rspace)
1970 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
1971 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
1972
1973 ! Stress tensor XC-functional GGA contribution
1974 virial%pv_exc = virial%pv_exc - virial%pv_xc
1975 virial%pv_virial = virial%pv_virial - virial%pv_xc
1976
1977 IF (debug_stress) THEN
1978 stdeb = -1.0_dp*fconv*virial%pv_xc
1979 CALL para_env%sum(stdeb)
1980 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1981 'STRESS| GGA 1st E_xc[Pin] ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1982 END IF
1983
1984 IF (ASSOCIATED(v_rspace)) THEN
1985 DO ispin = 1, nspins
1986 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
1987 END DO
1988 DEALLOCATE (v_rspace)
1989 END IF
1990 IF (ASSOCIATED(v_tau_rspace)) THEN
1991 DO ispin = 1, nspins
1992 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1993 END DO
1994 DEALLOCATE (v_tau_rspace)
1995 END IF
1996 CALL pw_zero(rhodn_tot_gspace)
1997
1998 END IF
1999
2000 ! rho_out - rho_in
2001 DO ispin = 1, nspins
2002 CALL pw_axpy(rho_r(ispin), rhoout_r(ispin), -1.0_dp)
2003 CALL pw_axpy(rho_g(ispin), rhoout_g(ispin), -1.0_dp)
2004 CALL pw_axpy(rhoout_g(ispin), rhodn_tot_gspace)
2005 IF (dft_control%use_kinetic_energy_density) CALL pw_axpy(tau_r(ispin), tauout_r(ispin), -1.0_dp)
2006 END DO
2007
2008 ! calculate associated hartree potential
2009 IF (use_virial) THEN
2010
2011 ! Stress tensor - 2nd derivative Volume and Green function contribution
2012 h_stress(:, :) = 0.0_dp
2013 CALL pw_poisson_solve(poisson_env, &
2014 density=rhodn_tot_gspace, & ! delta_n
2015 ehartree=dehartree, &
2016 vhartree=v_hartree_gspace, & ! v_H[delta_n]
2017 h_stress=h_stress, &
2018 aux_density=rho_tot_gspace) ! n_in
2019
2020 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
2021
2022 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
2023 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
2024
2025 IF (debug_stress) THEN
2026 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
2027 CALL para_env%sum(stdeb)
2028 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2029 'STRESS| GREEN 2nd V_H[dP]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2030 END IF
2031
2032 ELSE
2033 ! v_H[dn]
2034 CALL pw_poisson_solve(poisson_env, rhodn_tot_gspace, dehartree, &
2035 v_hartree_gspace)
2036 END IF
2037
2038 CALL pw_transfer(v_hartree_gspace, dv_hartree_rspace)
2039 CALL pw_scale(dv_hartree_rspace, dv_hartree_rspace%pw_grid%dvol)
2040 ! Getting nuclear force contribution from the core charge density
2041 ! Vh(rho_in + rho_c) + Vh(rho_out - rho_in)
2042 CALL pw_transfer(v_hartree_rspace, vtot_rspace)
2043 CALL pw_axpy(dv_hartree_rspace, vtot_rspace)
2044 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
2045 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
2046 CALL integrate_v_core_rspace(vtot_rspace, qs_env)
2047 IF (debug_forces) THEN
2048 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
2049 CALL para_env%sum(fodeb)
2050 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vtot*dncore", fodeb
2051 END IF
2052 IF (debug_stress .AND. use_virial) THEN
2053 stdeb = fconv*(virial%pv_ehartree - stdeb)
2054 CALL para_env%sum(stdeb)
2055 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2056 'STRESS| Vtot*dncore', one_third_sum_diag(stdeb), det_3x3(stdeb)
2057 END IF
2058 !
2059 ! Pulay force from Tr P_in (V_H(drho)+ Fxc(rho_in)*drho)
2060 ! RHS of CPKS equations: (V_H(drho)+ Fxc(rho_in)*drho)*C0
2061 ! Fxc*drho term
2062 xc_section => ec_env%xc_section
2063
2064 IF (use_virial) virial%pv_xc = 0.0_dp
2065 NULLIFY (v_xc, v_xc_tau)
2066 CALL create_kernel(qs_env, &
2067 vxc=v_xc, &
2068 vxc_tau=v_xc_tau, &
2069 rho=rho, &
2070 rho1_r=rhoout_r, &
2071 rho1_g=rhoout_g, &
2072 tau1_r=tauout_r, &
2073 xc_section=xc_section, &
2074 compute_virial=use_virial, &
2075 virial_xc=virial%pv_xc)
2076
2077 IF (use_virial) THEN
2078 ! Stress-tensor XC-functional 2nd GGA terms
2079 virial%pv_exc = virial%pv_exc + virial%pv_xc
2080 virial%pv_virial = virial%pv_virial + virial%pv_xc
2081 END IF
2082 IF (debug_stress .AND. use_virial) THEN
2083 stdeb = 1.0_dp*fconv*virial%pv_xc
2084 CALL para_env%sum(stdeb)
2085 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2086 'STRESS| GGA 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2087 END IF
2088 !
2089 CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s_kp=matrix_s)
2090 NULLIFY (ec_env%matrix_hz)
2091 CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
2092 DO ispin = 1, nspins
2093 ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
2094 CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, template=matrix_s(1, 1)%matrix)
2095 CALL dbcsr_copy(ec_env%matrix_hz(ispin)%matrix, matrix_s(1, 1)%matrix)
2096 CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
2097 END DO
2098 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2099 ! vtot = v_xc(ispin) + dv_hartree
2100 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2101 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2102
2103 ! Stress-tensor 2nd derivative integral contribution
2104 IF (use_virial) THEN
2105 pv_loc = virial%pv_virial
2106 END IF
2107
2108 DO ispin = 1, nspins
2109 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2110 CALL pw_axpy(dv_hartree_rspace, v_xc(ispin))
2111 CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
2112 hmat=ec_env%matrix_hz(ispin), &
2113 pmat=matrix_p(ispin, 1), &
2114 qs_env=qs_env, &
2115 calculate_forces=.true.)
2116 END DO
2117
2118 IF (debug_forces) THEN
2119 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2120 CALL para_env%sum(fodeb)
2121 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKdrho", fodeb
2122 END IF
2123 IF (debug_stress .AND. use_virial) THEN
2124 stdeb = fconv*(virial%pv_virial - stdeb)
2125 CALL para_env%sum(stdeb)
2126 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2127 'STRESS| INT 2nd f_Hxc[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2128 END IF
2129
2130 IF (ASSOCIATED(v_xc_tau)) THEN
2131 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2132 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2133
2134 DO ispin = 1, nspins
2135 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2136 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
2137 hmat=ec_env%matrix_hz(ispin), &
2138 pmat=matrix_p(ispin, 1), &
2139 qs_env=qs_env, &
2140 compute_tau=.true., &
2141 calculate_forces=.true.)
2142 END DO
2143
2144 IF (debug_forces) THEN
2145 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2146 CALL para_env%sum(fodeb)
2147 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKtaudtau", fodeb
2148 END IF
2149 IF (debug_stress .AND. use_virial) THEN
2150 stdeb = fconv*(virial%pv_virial - stdeb)
2151 CALL para_env%sum(stdeb)
2152 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2153 'STRESS| INT 2nd f_xctau[dP]*Pin ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2154 END IF
2155 END IF
2156 ! Stress-tensor 2nd derivative integral contribution
2157 IF (use_virial) THEN
2158 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2159 END IF
2160
2161 ! initialize srcm matrix
2162 NULLIFY (scrm)
2163 CALL dbcsr_allocate_matrix_set(scrm, nspins)
2164 DO ispin = 1, nspins
2165 ALLOCATE (scrm(ispin)%matrix)
2166 CALL dbcsr_create(scrm(ispin)%matrix, template=ec_env%matrix_ks(ispin, 1)%matrix)
2167 CALL dbcsr_copy(scrm(ispin)%matrix, ec_env%matrix_ks(ispin, 1)%matrix)
2168 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
2169 END DO
2170
2171 ! v_rspace and v_tau_rspace are generated from the auxbas pool
2172 NULLIFY (v_rspace, v_tau_rspace)
2173
2174 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=ec_env%xc_section, &
2175 vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=eexc, just_energy=.false.)
2176
2177 IF (use_virial) THEN
2178 eexc = 0.0_dp
2179 IF (ASSOCIATED(v_rspace)) THEN
2180 DO ispin = 1, nspins
2181 ! 2nd deriv xc-volume term
2182 eexc = eexc + pw_integral_ab(rhoout_r(ispin), v_rspace(ispin))
2183 END DO
2184 END IF
2185 IF (ASSOCIATED(v_tau_rspace)) THEN
2186 DO ispin = 1, nspins
2187 ! 2nd deriv xc-volume term
2188 eexc = eexc + pw_integral_ab(tauout_r(ispin), v_tau_rspace(ispin))
2189 END DO
2190 END IF
2191 END IF
2192
2193 IF (.NOT. ASSOCIATED(v_rspace)) THEN
2194 ALLOCATE (v_rspace(nspins))
2195 DO ispin = 1, nspins
2196 CALL auxbas_pw_pool%create_pw(v_rspace(ispin))
2197 CALL pw_zero(v_rspace(ispin))
2198 END DO
2199 END IF
2200
2201 ! Stress-tensor contribution derivative of integrand
2202 ! int v_Hxc[n^în]*n^out
2203 IF (use_virial) THEN
2204 pv_loc = virial%pv_virial
2205 END IF
2206
2207 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2208 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2209 DO ispin = 1, nspins
2210 ! Add v_hartree + v_xc = v_rspace
2211 CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
2212 CALL pw_axpy(v_hartree_rspace, v_rspace(ispin))
2213 ! integrate over potential <a|V|b>
2214 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
2215 hmat=scrm(ispin), &
2216 pmat=ec_env%matrix_p(ispin, 1), &
2217 qs_env=qs_env, &
2218 calculate_forces=.true., &
2219 basis_type="HARRIS", &
2220 task_list_external=ec_env%task_list)
2221 END DO
2222
2223 IF (debug_forces) THEN
2224 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2225 CALL para_env%sum(fodeb)
2226 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc ", fodeb
2227 END IF
2228 IF (debug_stress .AND. use_virial) THEN
2229 stdeb = fconv*(virial%pv_virial - stdeb)
2230 CALL para_env%sum(stdeb)
2231 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2232 'STRESS| INT Pout*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2233 END IF
2234
2235 ! Stress-tensor
2236 IF (use_virial) THEN
2237 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2238 END IF
2239
2240 IF (ASSOCIATED(v_tau_rspace)) THEN
2241 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2242 DO ispin = 1, nspins
2243 ! integrate over Tau-potential <nabla.a|V|nabla.b>
2244 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
2245 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
2246 hmat=scrm(ispin), &
2247 pmat=ec_env%matrix_p(ispin, 1), &
2248 qs_env=qs_env, &
2249 calculate_forces=.true., &
2250 compute_tau=.true., &
2251 basis_type="HARRIS", &
2252 task_list_external=ec_env%task_list)
2253 END DO
2254 IF (debug_forces) THEN
2255 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2256 CALL para_env%sum(fodeb)
2257 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dVhxc_tau ", fodeb
2258 END IF
2259 END IF
2260
2261!------------------------------------------------------------------------------
2262! HFX direct force
2263!------------------------------------------------------------------------------
2264
2265 ! If hybrid functional
2266 ec_hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%XC%HF")
2267 CALL section_vals_get(ec_hfx_sections, explicit=do_ec_hfx)
2268
2269 IF (do_ec_hfx) THEN
2270
2271 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2272 IF (use_virial) virial%pv_fock_4c = 0.0_dp
2273
2274 CALL calculate_exx(qs_env=qs_env, &
2275 unit_nr=iounit, &
2276 hfx_sections=ec_hfx_sections, &
2277 x_data=ec_env%x_data, &
2278 do_gw=.false., &
2279 do_admm=ec_env%do_ec_admm, &
2280 calc_forces=.true., &
2281 reuse_hfx=ec_env%reuse_hfx, &
2282 do_im_time=.false., &
2283 e_ex_from_gw=dummy_real, &
2284 e_admm_from_gw=dummy_real2, &
2285 t3=dummy_real)
2286
2287 IF (use_virial) THEN
2288 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2289 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2290 virial%pv_calculate = .false.
2291 END IF
2292 IF (debug_forces) THEN
2293 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2294 CALL para_env%sum(fodeb)
2295 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*hfx ", fodeb
2296 END IF
2297 IF (debug_stress .AND. use_virial) THEN
2298 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2299 CALL para_env%sum(stdeb)
2300 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2301 'STRESS| Pout*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2302 END IF
2303
2304 END IF
2305
2306!------------------------------------------------------------------------------
2307
2308 ! delete scrm matrix
2309 CALL dbcsr_deallocate_matrix_set(scrm)
2310
2311 ! return pw grids
2312 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
2313 DO ispin = 1, nspins
2314 CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
2315 IF (ASSOCIATED(v_tau_rspace)) THEN
2316 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
2317 END IF
2318 END DO
2319 IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
2320
2321 ! Core overlap
2322 IF (debug_forces) fodeb(1:3) = force(1)%core_overlap(1:3, 1)
2323 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ecore_overlap
2324 CALL calculate_ecore_overlap(qs_env, para_env, .true., e_overlap_core=eovrl)
2325 IF (debug_forces) THEN
2326 fodeb(1:3) = force(1)%core_overlap(1:3, 1) - fodeb(1:3)
2327 CALL para_env%sum(fodeb)
2328 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: CoreOverlap", fodeb
2329 END IF
2330 IF (debug_stress .AND. use_virial) THEN
2331 stdeb = fconv*(stdeb - virial%pv_ecore_overlap)
2332 CALL para_env%sum(stdeb)
2333 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2334 'STRESS| CoreOverlap ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2335 END IF
2336
2337 IF (debug_forces) THEN
2338 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2339 ALLOCATE (ftot(3, natom))
2340 CALL total_qs_force(ftot, force, atomic_kind_set)
2341 fodeb(1:3) = ftot(1:3, 1)
2342 DEALLOCATE (ftot)
2343 CALL para_env%sum(fodeb)
2344 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Force Explicit", fodeb
2345 END IF
2346
2347 DEALLOCATE (v_rspace)
2348 !
2349 CALL auxbas_pw_pool%give_back_pw(dv_hartree_rspace)
2350 CALL auxbas_pw_pool%give_back_pw(vtot_rspace)
2351 DO ispin = 1, nspins
2352 CALL auxbas_pw_pool%give_back_pw(rhoout_r(ispin))
2353 CALL auxbas_pw_pool%give_back_pw(rhoout_g(ispin))
2354 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2355 END DO
2356 DEALLOCATE (rhoout_r, rhoout_g, v_xc)
2357 IF (ASSOCIATED(tauout_r)) THEN
2358 DO ispin = 1, nspins
2359 CALL auxbas_pw_pool%give_back_pw(tauout_r(ispin))
2360 END DO
2361 DEALLOCATE (tauout_r)
2362 END IF
2363 IF (ASSOCIATED(v_xc_tau)) THEN
2364 DO ispin = 1, nspins
2365 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2366 END DO
2367 DEALLOCATE (v_xc_tau)
2368 END IF
2369 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
2370 CALL auxbas_pw_pool%give_back_pw(rhodn_tot_gspace)
2371
2372 ! Stress tensor - volume terms need to be stored,
2373 ! for a sign correction in QS at the end of qs_force
2374 IF (use_virial) THEN
2375 IF (qs_env%energy_correction) THEN
2376 ec_env%ehartree = ehartree + dehartree
2377 ec_env%exc = exc + eexc
2378 END IF
2379 END IF
2380
2381 IF (debug_stress .AND. use_virial) THEN
2382 ! In total: -1.0*E_H
2383 stdeb = -1.0_dp*fconv*ehartree
2384 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2385 'STRESS| VOL 1st v_H[n_in]*n_out', one_third_sum_diag(stdeb), det_3x3(stdeb)
2386
2387 stdeb = -1.0_dp*fconv*exc
2388 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2389 'STRESS| VOL 1st E_XC[n_in]', one_third_sum_diag(stdeb), det_3x3(stdeb)
2390
2391 stdeb = -1.0_dp*fconv*dehartree
2392 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2393 'STRESS| VOL 2nd v_H[dP]*n_in', one_third_sum_diag(stdeb), det_3x3(stdeb)
2394
2395 stdeb = -1.0_dp*fconv*eexc
2396 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2397 'STRESS| VOL 2nd v_XC[n_in]*dP', one_third_sum_diag(stdeb), det_3x3(stdeb)
2398
2399 ! For debugging, create a second virial environment,
2400 ! apply volume terms immediately
2401 block
2402 TYPE(virial_type) :: virdeb
2403 virdeb = virial
2404
2405 CALL para_env%sum(virdeb%pv_overlap)
2406 CALL para_env%sum(virdeb%pv_ekinetic)
2407 CALL para_env%sum(virdeb%pv_ppl)
2408 CALL para_env%sum(virdeb%pv_ppnl)
2409 CALL para_env%sum(virdeb%pv_ecore_overlap)
2410 CALL para_env%sum(virdeb%pv_ehartree)
2411 CALL para_env%sum(virdeb%pv_exc)
2412 CALL para_env%sum(virdeb%pv_exx)
2413 CALL para_env%sum(virdeb%pv_vdw)
2414 CALL para_env%sum(virdeb%pv_mp2)
2415 CALL para_env%sum(virdeb%pv_nlcc)
2416 CALL para_env%sum(virdeb%pv_gapw)
2417 CALL para_env%sum(virdeb%pv_lrigpw)
2418 CALL para_env%sum(virdeb%pv_virial)
2419 CALL symmetrize_virial(virdeb)
2420
2421 ! apply stress-tensor 1st and 2nd volume terms
2422 DO i = 1, 3
2423 virdeb%pv_ehartree(i, i) = virdeb%pv_ehartree(i, i) - 2.0_dp*(ehartree + dehartree)
2424 virdeb%pv_virial(i, i) = virdeb%pv_virial(i, i) - exc - eexc &
2425 - 2.0_dp*(ehartree + dehartree)
2426 virdeb%pv_exc(i, i) = virdeb%pv_exc(i, i) - exc - eexc
2427 ! The factor 2 is a hack. It compensates the plus sign in h_stress/pw_poisson_solve.
2428 ! The sign in pw_poisson_solve is correct for FIST, but not for QS.
2429 ! There should be a more elegant solution to that ...
2430 END DO
2431
2432 CALL para_env%sum(sttot)
2433 stdeb = fconv*(virdeb%pv_virial - sttot)
2434 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2435 'STRESS| Explicit electronic stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2436
2437 stdeb = fconv*(virdeb%pv_virial)
2438 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2439 'STRESS| Explicit total stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2440
2441 CALL write_stress_tensor_components(virdeb, iounit, cell, unit_string)
2442 CALL write_stress_tensor(virdeb%pv_virial, iounit, cell, unit_string, .false.)
2443
2444 END block
2445 END IF
2446
2447 CALL timestop(handle)
2448
2449 END SUBROUTINE ec_build_ks_matrix_force
2450
2451! **************************************************************************************************
2452!> \brief Solve KS equation for a given matrix
2453!> \param qs_env ...
2454!> \param ec_env ...
2455!> \par History
2456!> 03.2014 created [JGH]
2457!> \author JGH
2458! **************************************************************************************************
2459 SUBROUTINE ec_ks_solver(qs_env, ec_env)
2460
2461 TYPE(qs_environment_type), POINTER :: qs_env
2462 TYPE(energy_correction_type), POINTER :: ec_env
2463
2464 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_ks_solver'
2465
2466 CHARACTER(LEN=default_string_length) :: headline
2467 INTEGER :: handle, ispin, nspins
2468 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, pmat, smat, wmat
2469 TYPE(dft_control_type), POINTER :: dft_control
2470
2471 CALL timeset(routinen, handle)
2472
2473 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2474 nspins = dft_control%nspins
2475
2476 ! create density matrix
2477 IF (.NOT. ASSOCIATED(ec_env%matrix_p)) THEN
2478 headline = "DENSITY MATRIX"
2479 CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
2480 DO ispin = 1, nspins
2481 ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
2482 CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=trim(headline), &
2483 template=ec_env%matrix_s(1, 1)%matrix)
2484 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, ec_env%sab_orb)
2485 END DO
2486 END IF
2487 ! create energy weighted density matrix
2488 IF (.NOT. ASSOCIATED(ec_env%matrix_w)) THEN
2489 headline = "ENERGY WEIGHTED DENSITY MATRIX"
2490 CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
2491 DO ispin = 1, nspins
2492 ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
2493 CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=trim(headline), &
2494 template=ec_env%matrix_s(1, 1)%matrix)
2495 CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, ec_env%sab_orb)
2496 END DO
2497 END IF
2498
2499 IF (ec_env%mao) THEN
2500 CALL mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2501 ELSE
2502 ksmat => ec_env%matrix_ks
2503 smat => ec_env%matrix_s
2504 pmat => ec_env%matrix_p
2505 wmat => ec_env%matrix_w
2506 END IF
2507
2508 SELECT CASE (ec_env%ks_solver)
2509 CASE (ec_diagonalization)
2510 CALL ec_diag_solver(qs_env, ksmat, smat, pmat, wmat)
2511 CASE (ec_ot_diag)
2512 CALL ec_ot_diag_solver(qs_env, ec_env, ksmat, smat, pmat, wmat)
2513 CASE (ec_matrix_sign, ec_matrix_trs4, ec_matrix_tc2)
2514 CALL ec_ls_init(qs_env, ksmat, smat)
2515 CALL ec_ls_solver(qs_env, pmat, wmat, ec_ls_method=ec_env%ks_solver)
2516 CASE DEFAULT
2517 cpassert(.false.)
2518 END SELECT
2519
2520 IF (ec_env%mao) THEN
2521 CALL mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2522 END IF
2523
2524 CALL timestop(handle)
2525
2526 END SUBROUTINE ec_ks_solver
2527
2528! **************************************************************************************************
2529!> \brief Create matrices with MAO sizes
2530!> \param ec_env ...
2531!> \param ksmat ...
2532!> \param smat ...
2533!> \param pmat ...
2534!> \param wmat ...
2535!> \par History
2536!> 08.2016 created [JGH]
2537!> \author JGH
2538! **************************************************************************************************
2539 SUBROUTINE mao_create_matrices(ec_env, ksmat, smat, pmat, wmat)
2540
2541 TYPE(energy_correction_type), POINTER :: ec_env
2542 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat
2543
2544 CHARACTER(LEN=*), PARAMETER :: routinen = 'mao_create_matrices'
2545
2546 INTEGER :: handle, ispin, nspins
2547 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes
2548 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2549 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
2550 TYPE(dbcsr_type) :: cgmat
2551
2552 CALL timeset(routinen, handle)
2553
2554 mao_coef => ec_env%mao_coef
2555
2556 NULLIFY (ksmat, smat, pmat, wmat)
2557 nspins = SIZE(ec_env%matrix_ks, 1)
2558 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
2559 CALL dbcsr_allocate_matrix_set(ksmat, nspins, 1)
2560 CALL dbcsr_allocate_matrix_set(smat, nspins, 1)
2561 DO ispin = 1, nspins
2562 ALLOCATE (ksmat(ispin, 1)%matrix)
2563 CALL dbcsr_create(ksmat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO KS mat", &
2564 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2565 col_blk_size=col_blk_sizes)
2566 ALLOCATE (smat(ispin, 1)%matrix)
2567 CALL dbcsr_create(smat(ispin, 1)%matrix, dist=dbcsr_dist, name="MAO S mat", &
2568 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, &
2569 col_blk_size=col_blk_sizes)
2570 END DO
2571 !
2572 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
2573 DO ispin = 1, nspins
2574 CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_s(1, 1)%matrix, mao_coef(ispin)%matrix, &
2575 0.0_dp, cgmat)
2576 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, smat(ispin, 1)%matrix)
2577 CALL dbcsr_multiply("N", "N", 1.0_dp, ec_env%matrix_ks(1, 1)%matrix, mao_coef(ispin)%matrix, &
2578 0.0_dp, cgmat)
2579 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, ksmat(ispin, 1)%matrix)
2580 END DO
2581 CALL dbcsr_release(cgmat)
2582
2583 CALL dbcsr_allocate_matrix_set(pmat, nspins, 1)
2584 DO ispin = 1, nspins
2585 ALLOCATE (pmat(ispin, 1)%matrix)
2586 CALL dbcsr_create(pmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name="MAO P mat")
2587 CALL cp_dbcsr_alloc_block_from_nbl(pmat(ispin, 1)%matrix, ec_env%sab_orb)
2588 END DO
2589
2590 CALL dbcsr_allocate_matrix_set(wmat, nspins, 1)
2591 DO ispin = 1, nspins
2592 ALLOCATE (wmat(ispin, 1)%matrix)
2593 CALL dbcsr_create(wmat(ispin, 1)%matrix, template=smat(1, 1)%matrix, name="MAO W mat")
2594 CALL cp_dbcsr_alloc_block_from_nbl(wmat(ispin, 1)%matrix, ec_env%sab_orb)
2595 END DO
2596
2597 CALL timestop(handle)
2598
2599 END SUBROUTINE mao_create_matrices
2600
2601! **************************************************************************************************
2602!> \brief Release matrices with MAO sizes
2603!> \param ec_env ...
2604!> \param ksmat ...
2605!> \param smat ...
2606!> \param pmat ...
2607!> \param wmat ...
2608!> \par History
2609!> 08.2016 created [JGH]
2610!> \author JGH
2611! **************************************************************************************************
2612 SUBROUTINE mao_release_matrices(ec_env, ksmat, smat, pmat, wmat)
2613
2614 TYPE(energy_correction_type), POINTER :: ec_env
2615 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, smat, pmat, wmat
2616
2617 CHARACTER(LEN=*), PARAMETER :: routinen = 'mao_release_matrices'
2618
2619 INTEGER :: handle, ispin, nspins
2620 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef
2621 TYPE(dbcsr_type) :: cgmat
2622
2623 CALL timeset(routinen, handle)
2624
2625 mao_coef => ec_env%mao_coef
2626 nspins = SIZE(mao_coef, 1)
2627
2628 ! save pmat/wmat in full basis format
2629 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix)
2630 DO ispin = 1, nspins
2631 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, pmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2632 CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2633 ec_env%matrix_p(ispin, 1)%matrix, retain_sparsity=.true.)
2634 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef(ispin)%matrix, wmat(ispin, 1)%matrix, 0.0_dp, cgmat)
2635 CALL dbcsr_multiply("N", "T", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, &
2636 ec_env%matrix_w(ispin, 1)%matrix, retain_sparsity=.true.)
2637 END DO
2638 CALL dbcsr_release(cgmat)
2639
2640 CALL dbcsr_deallocate_matrix_set(ksmat)
2641 CALL dbcsr_deallocate_matrix_set(smat)
2642 CALL dbcsr_deallocate_matrix_set(pmat)
2643 CALL dbcsr_deallocate_matrix_set(wmat)
2644
2645 CALL timestop(handle)
2646
2647 END SUBROUTINE mao_release_matrices
2648
2649! **************************************************************************************************
2650!> \brief Solve KS equation using diagonalization
2651!> \param qs_env ...
2652!> \param matrix_ks ...
2653!> \param matrix_s ...
2654!> \param matrix_p ...
2655!> \param matrix_w ...
2656!> \par History
2657!> 03.2014 created [JGH]
2658!> \author JGH
2659! **************************************************************************************************
2660 SUBROUTINE ec_diag_solver(qs_env, matrix_ks, matrix_s, matrix_p, matrix_w)
2661
2662 TYPE(qs_environment_type), POINTER :: qs_env
2663 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_p, matrix_w
2664
2665 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_diag_solver'
2666
2667 INTEGER :: handle, ispin, nmo(2), nsize, nspins
2668 REAL(kind=dp) :: eps_filter, focc(2)
2669 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
2670 TYPE(cp_blacs_env_type), POINTER :: blacs_env
2671 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2672 TYPE(cp_fm_type) :: fm_ks, fm_mo, fm_ortho
2673 TYPE(dbcsr_type), POINTER :: buf1_dbcsr, buf2_dbcsr, buf3_dbcsr, &
2674 ortho_dbcsr, ref_matrix
2675 TYPE(dft_control_type), POINTER :: dft_control
2676 TYPE(mp_para_env_type), POINTER :: para_env
2677
2678 CALL timeset(routinen, handle)
2679
2680 NULLIFY (blacs_env, para_env)
2681 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
2682
2683 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
2684 eps_filter = dft_control%qs_control%eps_filter_matrix
2685 nspins = dft_control%nspins
2686
2687 nmo = 0
2688 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
2689 focc = 1._dp
2690 IF (nspins == 1) THEN
2691 focc = 2._dp
2692 nmo(1) = nmo(1)/2
2693 END IF
2694
2695 CALL dbcsr_get_info(matrix_ks(1, 1)%matrix, nfullrows_total=nsize)
2696 ALLOCATE (eigenvalues(nsize))
2697
2698 NULLIFY (fm_struct, ref_matrix)
2699 CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nsize, &
2700 ncol_global=nsize, para_env=para_env)
2701 CALL cp_fm_create(fm_ortho, fm_struct)
2702 CALL cp_fm_create(fm_ks, fm_struct)
2703 CALL cp_fm_create(fm_mo, fm_struct)
2704 CALL cp_fm_struct_release(fm_struct)
2705
2706 ! factorization
2707 ref_matrix => matrix_s(1, 1)%matrix
2708 NULLIFY (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2709 CALL dbcsr_init_p(ortho_dbcsr)
2710 CALL dbcsr_create(ortho_dbcsr, template=ref_matrix, &
2711 matrix_type=dbcsr_type_no_symmetry)
2712 CALL dbcsr_init_p(buf1_dbcsr)
2713 CALL dbcsr_create(buf1_dbcsr, template=ref_matrix, &
2714 matrix_type=dbcsr_type_no_symmetry)
2715 CALL dbcsr_init_p(buf2_dbcsr)
2716 CALL dbcsr_create(buf2_dbcsr, template=ref_matrix, &
2717 matrix_type=dbcsr_type_no_symmetry)
2718 CALL dbcsr_init_p(buf3_dbcsr)
2719 CALL dbcsr_create(buf3_dbcsr, template=ref_matrix, &
2720 matrix_type=dbcsr_type_no_symmetry)
2721
2722 ref_matrix => matrix_s(1, 1)%matrix
2723 CALL copy_dbcsr_to_fm(ref_matrix, fm_ortho)
2724 CALL cp_fm_cholesky_decompose(fm_ortho)
2725 CALL cp_fm_triangular_invert(fm_ortho)
2726 CALL cp_fm_set_all(fm_ks, 0.0_dp)
2727 CALL cp_fm_to_fm_triangular(fm_ortho, fm_ks, "U")
2728 CALL copy_fm_to_dbcsr(fm_ks, ortho_dbcsr)
2729 DO ispin = 1, nspins
2730 ! calculate ZHZ(T)
2731 CALL dbcsr_desymmetrize(matrix_ks(ispin, 1)%matrix, buf1_dbcsr)
2732 CALL dbcsr_multiply("N", "N", 1.0_dp, buf1_dbcsr, ortho_dbcsr, &
2733 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2734 CALL dbcsr_multiply("T", "N", 1.0_dp, ortho_dbcsr, buf2_dbcsr, &
2735 0.0_dp, buf1_dbcsr, filter_eps=eps_filter)
2736 ! copy to fm format
2737 CALL copy_dbcsr_to_fm(buf1_dbcsr, fm_ks)
2738 CALL choose_eigv_solver(fm_ks, fm_mo, eigenvalues)
2739 ! back transform of mos c = Z(T)*c
2740 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2741 CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2742 0.0_dp, buf2_dbcsr, filter_eps=eps_filter)
2743 ! density matrix
2744 CALL dbcsr_set(matrix_p(ispin, 1)%matrix, 0.0_dp)
2745 CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf2_dbcsr, &
2746 1.0_dp, matrix_p(ispin, 1)%matrix, &
2747 retain_sparsity=.true., last_k=nmo(ispin))
2748
2749 ! energy weighted density matrix
2750 CALL dbcsr_set(matrix_w(ispin, 1)%matrix, 0.0_dp)
2751 CALL cp_fm_column_scale(fm_mo, eigenvalues)
2752 CALL copy_fm_to_dbcsr(fm_mo, buf1_dbcsr)
2753 CALL dbcsr_multiply("N", "N", 1.0_dp, ortho_dbcsr, buf1_dbcsr, &
2754 0.0_dp, buf3_dbcsr, filter_eps=eps_filter)
2755 CALL dbcsr_multiply("N", "T", focc(ispin), buf2_dbcsr, buf3_dbcsr, &
2756 1.0_dp, matrix_w(ispin, 1)%matrix, &
2757 retain_sparsity=.true., last_k=nmo(ispin))
2758 END DO
2759
2760 CALL cp_fm_release(fm_ks)
2761 CALL cp_fm_release(fm_mo)
2762 CALL cp_fm_release(fm_ortho)
2763 CALL dbcsr_release(ortho_dbcsr)
2764 CALL dbcsr_release(buf1_dbcsr)
2765 CALL dbcsr_release(buf2_dbcsr)
2766 CALL dbcsr_release(buf3_dbcsr)
2767 DEALLOCATE (ortho_dbcsr, buf1_dbcsr, buf2_dbcsr, buf3_dbcsr)
2768 DEALLOCATE (eigenvalues)
2769
2770 CALL timestop(handle)
2771
2772 END SUBROUTINE ec_diag_solver
2773
2774! **************************************************************************************************
2775!> \brief Calculate the energy correction
2776!> \param ec_env ...
2777!> \param unit_nr ...
2778!> \author Creation (03.2014,JGH)
2779! **************************************************************************************************
2780 SUBROUTINE ec_energy(ec_env, unit_nr)
2781 TYPE(energy_correction_type) :: ec_env
2782 INTEGER, INTENT(IN) :: unit_nr
2783
2784 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_energy'
2785
2786 INTEGER :: handle, ispin, nspins
2787 REAL(kind=dp) :: eband, trace
2788
2789 CALL timeset(routinen, handle)
2790
2791 nspins = SIZE(ec_env%matrix_p, 1)
2792 DO ispin = 1, nspins
2793 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace)
2794 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T65,F16.10)') 'Tr[PS] ', trace
2795 END DO
2796
2797 ! Total energy depends on energy correction method
2798 SELECT CASE (ec_env%energy_functional)
2799 CASE (ec_functional_harris)
2800
2801 ! Get energy of "band structure" term
2802 eband = 0.0_dp
2803 DO ispin = 1, nspins
2804 CALL dbcsr_dot(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%matrix_p(ispin, 1)%matrix, trace)
2805 eband = eband + trace
2806 END DO
2807 ec_env%eband = eband + ec_env%efield_nuclear
2808
2809 ! Add Harris functional "correction" terms
2810 ec_env%etotal = ec_env%eband + ec_env%ehartree + ec_env%exc - ec_env%vhxc + ec_env%edispersion - ec_env%ex
2811 IF (unit_nr > 0) THEN
2812 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Eband ", ec_env%eband
2813 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree
2814 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc
2815 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex
2816 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Evhxc ", ec_env%vhxc
2817 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion
2818 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Harris Functional ", ec_env%etotal
2819 END IF
2820
2821 CASE (ec_functional_dc)
2822
2823 ! Core hamiltonian energy
2824 CALL calculate_ptrace(ec_env%matrix_h, ec_env%matrix_p, ec_env%ecore, SIZE(ec_env%matrix_p, 1))
2825
2826 ec_env%ecore = ec_env%ecore + ec_env%efield_nuclear
2827 ec_env%etotal = ec_env%ecore + ec_env%ehartree + ec_env%exc + ec_env%edispersion &
2828 + ec_env%ex + ec_env%exc_aux_fit
2829
2830 IF (unit_nr > 0) THEN
2831 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ecore ", ec_env%ecore
2832 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ehartree ", ec_env%ehartree
2833 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc ", ec_env%exc
2834 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Ex ", ec_env%ex
2835 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Exc_aux_fit", ec_env%exc_aux_fit
2836 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Edisp ", ec_env%edispersion
2837 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal
2838 END IF
2839
2840 CASE (ec_functional_ext)
2841
2842 ec_env%etotal = ec_env%ex
2843 IF (unit_nr > 0) THEN
2844 WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal
2845 END IF
2846
2847 CASE DEFAULT
2848
2849 cpassert(.false.)
2850
2851 END SELECT
2852
2853 CALL timestop(handle)
2854
2855 END SUBROUTINE ec_energy
2856
2857! **************************************************************************************************
2858!> \brief builds either the full neighborlist or neighborlists of molecular
2859!> \brief subsets, depending on parameter values
2860!> \param qs_env ...
2861!> \param ec_env ...
2862!> \par History
2863!> 2012.07 created [Martin Haeufel]
2864!> 2016.07 Adapted for Harris functional [JGH]
2865!> \author Martin Haeufel
2866! **************************************************************************************************
2867 SUBROUTINE ec_build_neighborlist(qs_env, ec_env)
2868 TYPE(qs_environment_type), POINTER :: qs_env
2869 TYPE(energy_correction_type), POINTER :: ec_env
2870
2871 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_build_neighborlist'
2872
2873 INTEGER :: handle, ikind, nkind, zat
2874 LOGICAL :: all_potential_present, &
2875 gth_potential_present, &
2876 sgp_potential_present, &
2877 skip_load_balance_distributed
2878 LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, default_present, &
2879 orb_present, ppl_present, ppnl_present
2880 REAL(dp) :: subcells
2881 REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_radius, c_radius, orb_radius, &
2882 ppl_radius, ppnl_radius
2883 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
2884 TYPE(all_potential_type), POINTER :: all_potential
2885 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2886 TYPE(cell_type), POINTER :: cell
2887 TYPE(dft_control_type), POINTER :: dft_control
2888 TYPE(distribution_1d_type), POINTER :: distribution_1d
2889 TYPE(distribution_2d_type), POINTER :: distribution_2d
2890 TYPE(gth_potential_type), POINTER :: gth_potential
2891 TYPE(gto_basis_set_type), POINTER :: basis_set
2892 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
2893 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
2894 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2895 POINTER :: sab_cn, sab_vdw
2896 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2897 TYPE(qs_dispersion_type), POINTER :: dispersion_env
2898 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2899 TYPE(qs_kind_type), POINTER :: qs_kind
2900 TYPE(qs_ks_env_type), POINTER :: ks_env
2901 TYPE(sgp_potential_type), POINTER :: sgp_potential
2902
2903 CALL timeset(routinen, handle)
2904
2905 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
2906 CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present, &
2907 gth_potential_present=gth_potential_present, &
2908 sgp_potential_present=sgp_potential_present)
2909 nkind = SIZE(qs_kind_set)
2910 ALLOCATE (c_radius(nkind), default_present(nkind))
2911 ALLOCATE (orb_radius(nkind), all_radius(nkind), ppl_radius(nkind), ppnl_radius(nkind))
2912 ALLOCATE (orb_present(nkind), all_present(nkind), ppl_present(nkind), ppnl_present(nkind))
2913 ALLOCATE (pair_radius(nkind, nkind))
2914 ALLOCATE (atom2d(nkind))
2915
2916 CALL get_qs_env(qs_env, &
2917 atomic_kind_set=atomic_kind_set, &
2918 cell=cell, &
2919 distribution_2d=distribution_2d, &
2920 local_particles=distribution_1d, &
2921 particle_set=particle_set, &
2922 molecule_set=molecule_set)
2923
2924 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
2925 molecule_set, .false., particle_set)
2926
2927 DO ikind = 1, nkind
2928 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
2929 qs_kind => qs_kind_set(ikind)
2930 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="HARRIS")
2931 IF (ASSOCIATED(basis_set)) THEN
2932 orb_present(ikind) = .true.
2933 CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=orb_radius(ikind))
2934 ELSE
2935 orb_present(ikind) = .false.
2936 orb_radius(ikind) = 0.0_dp
2937 END IF
2938 CALL get_qs_kind(qs_kind, all_potential=all_potential, &
2939 gth_potential=gth_potential, sgp_potential=sgp_potential)
2940 IF (gth_potential_present .OR. sgp_potential_present) THEN
2941 IF (ASSOCIATED(gth_potential)) THEN
2942 CALL get_potential(potential=gth_potential, &
2943 ppl_present=ppl_present(ikind), &
2944 ppl_radius=ppl_radius(ikind), &
2945 ppnl_present=ppnl_present(ikind), &
2946 ppnl_radius=ppnl_radius(ikind))
2947 ELSE IF (ASSOCIATED(sgp_potential)) THEN
2948 CALL get_potential(potential=sgp_potential, &
2949 ppl_present=ppl_present(ikind), &
2950 ppl_radius=ppl_radius(ikind), &
2951 ppnl_present=ppnl_present(ikind), &
2952 ppnl_radius=ppnl_radius(ikind))
2953 ELSE
2954 ppl_present(ikind) = .false.
2955 ppl_radius(ikind) = 0.0_dp
2956 ppnl_present(ikind) = .false.
2957 ppnl_radius(ikind) = 0.0_dp
2958 END IF
2959 END IF
2960 ! Check the presence of an all electron potential or ERFC potential
2961 IF (all_potential_present .OR. sgp_potential_present) THEN
2962 all_present(ikind) = .false.
2963 all_radius(ikind) = 0.0_dp
2964 IF (ASSOCIATED(all_potential)) THEN
2965 all_present(ikind) = .true.
2966 CALL get_potential(potential=all_potential, core_charge_radius=all_radius(ikind))
2967 ELSE IF (ASSOCIATED(sgp_potential)) THEN
2968 IF (sgp_potential%ecp_local) THEN
2969 all_present(ikind) = .true.
2970 CALL get_potential(potential=sgp_potential, core_charge_radius=all_radius(ikind))
2971 END IF
2972 END IF
2973 END IF
2974 END DO
2975
2976 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
2977
2978 ! overlap
2979 CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
2980 CALL build_neighbor_lists(ec_env%sab_orb, particle_set, atom2d, cell, pair_radius, &
2981 subcells=subcells, nlname="sab_orb")
2982 ! pseudopotential/AE
2983 IF (all_potential_present .OR. sgp_potential_present) THEN
2984 IF (any(all_present)) THEN
2985 CALL pair_radius_setup(orb_present, all_present, orb_radius, all_radius, pair_radius)
2986 CALL build_neighbor_lists(ec_env%sac_ae, particle_set, atom2d, cell, pair_radius, &
2987 subcells=subcells, operator_type="ABC", nlname="sac_ae")
2988 END IF
2989 END IF
2990
2991 IF (gth_potential_present .OR. sgp_potential_present) THEN
2992 IF (any(ppl_present)) THEN
2993 CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
2994 CALL build_neighbor_lists(ec_env%sac_ppl, particle_set, atom2d, cell, pair_radius, &
2995 subcells=subcells, operator_type="ABC", nlname="sac_ppl")
2996 END IF
2997
2998 IF (any(ppnl_present)) THEN
2999 CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
3000 CALL build_neighbor_lists(ec_env%sap_ppnl, particle_set, atom2d, cell, pair_radius, &
3001 subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
3002 END IF
3003 END IF
3004
3005 ! Build the neighbor lists for the vdW pair potential
3006 c_radius(:) = 0.0_dp
3007 dispersion_env => ec_env%dispersion_env
3008 sab_vdw => dispersion_env%sab_vdw
3009 sab_cn => dispersion_env%sab_cn
3010 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
3011 c_radius(:) = dispersion_env%rc_disp
3012 default_present = .true. !include all atoms in vdW (even without basis)
3013 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
3014 CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
3015 subcells=subcells, operator_type="PP", nlname="sab_vdw")
3016 dispersion_env%sab_vdw => sab_vdw
3017 IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
3018 dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
3019 ! Build the neighbor lists for coordination numbers as needed by the DFT-D3 method
3020 DO ikind = 1, nkind
3021 CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
3022 c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
3023 END DO
3024 CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
3025 CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
3026 subcells=subcells, operator_type="PP", nlname="sab_cn")
3027 dispersion_env%sab_cn => sab_cn
3028 END IF
3029 END IF
3030
3031 ! Release work storage
3032 CALL atom2d_cleanup(atom2d)
3033 DEALLOCATE (atom2d)
3034 DEALLOCATE (orb_present, default_present, all_present, ppl_present, ppnl_present)
3035 DEALLOCATE (orb_radius, all_radius, ppl_radius, ppnl_radius, c_radius)
3036 DEALLOCATE (pair_radius)
3037
3038 ! Task list
3039 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
3040 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
3041 IF (ASSOCIATED(ec_env%task_list)) CALL deallocate_task_list(ec_env%task_list)
3042 CALL allocate_task_list(ec_env%task_list)
3043 CALL generate_qs_task_list(ks_env, ec_env%task_list, &
3044 reorder_rs_grid_ranks=.false., soft_valid=.false., &
3045 skip_load_balance_distributed=skip_load_balance_distributed, &
3046 basis_type="HARRIS", sab_orb_external=ec_env%sab_orb)
3047
3048 CALL timestop(handle)
3049
3050 END SUBROUTINE ec_build_neighborlist
3051
3052! **************************************************************************************************
3053!> \brief ...
3054!> \param qs_env ...
3055!> \param ec_env ...
3056! **************************************************************************************************
3057 SUBROUTINE ec_properties(qs_env, ec_env)
3058 TYPE(qs_environment_type), POINTER :: qs_env
3059 TYPE(energy_correction_type), POINTER :: ec_env
3060
3061 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_properties'
3062
3063 CHARACTER(LEN=8), DIMENSION(3) :: rlab
3064 CHARACTER(LEN=default_path_length) :: filename, my_pos_voro
3065 CHARACTER(LEN=default_string_length) :: description
3066 INTEGER :: akind, handle, i, ia, iatom, idir, ikind, iounit, ispin, maxmom, nspins, &
3067 reference, should_print_bqb, should_print_voro, unit_nr, unit_nr_voro
3068 LOGICAL :: append_voro, magnetic, periodic, &
3069 voro_print_txt
3070 REAL(kind=dp) :: charge, dd, focc, tmp
3071 REAL(kind=dp), DIMENSION(3) :: cdip, pdip, rcc, rdip, ria, tdip
3072 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
3073 TYPE(atomic_kind_type), POINTER :: atomic_kind
3074 TYPE(cell_type), POINTER :: cell
3075 TYPE(cp_logger_type), POINTER :: logger
3076 TYPE(cp_result_type), POINTER :: results
3077 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
3078 TYPE(dft_control_type), POINTER :: dft_control
3079 TYPE(distribution_1d_type), POINTER :: local_particles
3080 TYPE(mp_para_env_type), POINTER :: para_env
3081 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3082 TYPE(pw_env_type), POINTER :: pw_env
3083 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
3084 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3085 TYPE(pw_r3d_rs_type) :: rho_elec_rspace
3086 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3087 TYPE(section_vals_type), POINTER :: ec_section, print_key, print_key_bqb, &
3088 print_key_voro
3089
3090 CALL timeset(routinen, handle)
3091
3092 rlab(1) = "X"
3093 rlab(2) = "Y"
3094 rlab(3) = "Z"
3095
3096 logger => cp_get_default_logger()
3097 IF (logger%para_env%is_source()) THEN
3098 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
3099 ELSE
3100 iounit = -1
3101 END IF
3102
3103 NULLIFY (dft_control)
3104 CALL get_qs_env(qs_env, dft_control=dft_control)
3105 nspins = dft_control%nspins
3106
3107 ec_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION")
3108 print_key => section_vals_get_subs_vals(section_vals=ec_section, &
3109 subsection_name="PRINT%MOMENTS")
3110
3111 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
3112
3113 maxmom = section_get_ival(section_vals=ec_section, &
3114 keyword_name="PRINT%MOMENTS%MAX_MOMENT")
3115 periodic = section_get_lval(section_vals=ec_section, &
3116 keyword_name="PRINT%MOMENTS%PERIODIC")
3117 reference = section_get_ival(section_vals=ec_section, &
3118 keyword_name="PRINT%MOMENTS%REFERENCE")
3119 magnetic = section_get_lval(section_vals=ec_section, &
3120 keyword_name="PRINT%MOMENTS%MAGNETIC")
3121 NULLIFY (ref_point)
3122 CALL section_vals_val_get(ec_section, "PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
3123 unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=ec_section, &
3124 print_key_path="PRINT%MOMENTS", extension=".dat", &
3125 middle_name="moments", log_filename=.false.)
3126
3127 IF (iounit > 0) THEN
3128 IF (unit_nr /= iounit .AND. unit_nr > 0) THEN
3129 INQUIRE (unit=unit_nr, name=filename)
3130 WRITE (unit=iounit, fmt="(/,T2,A,2(/,T3,A),/)") &
3131 "MOMENTS", "The electric/magnetic moments are written to file:", &
3132 trim(filename)
3133 ELSE
3134 WRITE (unit=iounit, fmt="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
3135 END IF
3136 END IF
3137
3138 IF (periodic) THEN
3139 cpabort("Periodic moments not implemented with EC")
3140 ELSE
3141 cpassert(maxmom < 2)
3142 cpassert(.NOT. magnetic)
3143 IF (maxmom == 1) THEN
3144 CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
3145 ! reference point
3146 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3147 ! nuclear contribution
3148 cdip = 0.0_dp
3149 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
3150 qs_kind_set=qs_kind_set, local_particles=local_particles)
3151 DO ikind = 1, SIZE(local_particles%n_el)
3152 DO ia = 1, local_particles%n_el(ikind)
3153 iatom = local_particles%list(ikind)%array(ia)
3154 ! fold atomic positions back into unit cell
3155 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
3156 ria = ria - rcc
3157 atomic_kind => particle_set(iatom)%atomic_kind
3158 CALL get_atomic_kind(atomic_kind, kind_number=akind)
3159 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
3160 cdip(1:3) = cdip(1:3) - charge*ria(1:3)
3161 END DO
3162 END DO
3163 CALL para_env%sum(cdip)
3164 !
3165 ! direct density contribution
3166 CALL ec_efield_integrals(qs_env, ec_env, rcc)
3167 !
3168 pdip = 0.0_dp
3169 DO ispin = 1, nspins
3170 DO idir = 1, 3
3171 CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, &
3172 ec_env%efield%dipmat(idir)%matrix, tmp)
3173 pdip(idir) = pdip(idir) + tmp
3174 END DO
3175 END DO
3176 !
3177 ! response contribution
3178 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
3179 NULLIFY (moments)
3180 CALL dbcsr_allocate_matrix_set(moments, 4)
3181 DO i = 1, 4
3182 ALLOCATE (moments(i)%matrix)
3183 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
3184 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
3185 END DO
3186 CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
3187 !
3188 focc = 2.0_dp
3189 IF (nspins == 2) focc = 1.0_dp
3190 rdip = 0.0_dp
3191 DO ispin = 1, nspins
3192 DO idir = 1, 3
3193 CALL dbcsr_dot(ec_env%matrix_z(ispin)%matrix, moments(idir)%matrix, tmp)
3194 rdip(idir) = rdip(idir) + tmp
3195 END DO
3196 END DO
3197 CALL dbcsr_deallocate_matrix_set(moments)
3198 !
3199 tdip = -(rdip + pdip + cdip)
3200 IF (unit_nr > 0) THEN
3201 WRITE (unit_nr, "(T3,A)") "Dipoles are based on the traditional operator."
3202 dd = sqrt(sum(tdip(1:3)**2))*debye
3203 WRITE (unit_nr, "(T3,A)") "Dipole moment [Debye]"
3204 WRITE (unit_nr, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
3205 (trim(rlab(i)), "=", tdip(i)*debye, i=1, 3), "Total=", dd
3206 END IF
3207 END IF
3208 END IF
3209
3210 CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
3211 basis_section=ec_section, print_key_path="PRINT%MOMENTS")
3212 CALL get_qs_env(qs_env=qs_env, results=results)
3213 description = "[DIPOLE]"
3214 CALL cp_results_erase(results=results, description=description)
3215 CALL put_results(results=results, description=description, values=tdip(1:3))
3216 END IF
3217
3218 ! Do a Voronoi Integration or write a compressed BQB File
3219 print_key_voro => section_vals_get_subs_vals(ec_section, "PRINT%VORONOI")
3220 print_key_bqb => section_vals_get_subs_vals(ec_section, "PRINT%E_DENSITY_BQB")
3221 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_voro), cp_p_file)) THEN
3222 should_print_voro = 1
3223 ELSE
3224 should_print_voro = 0
3225 END IF
3226 IF (btest(cp_print_key_should_output(logger%iter_info, print_key_bqb), cp_p_file)) THEN
3227 should_print_bqb = 1
3228 ELSE
3229 should_print_bqb = 0
3230 END IF
3231 IF ((should_print_voro /= 0) .OR. (should_print_bqb /= 0)) THEN
3232
3233 CALL get_qs_env(qs_env=qs_env, &
3234 pw_env=pw_env)
3235 CALL pw_env_get(pw_env=pw_env, &
3236 auxbas_pw_pool=auxbas_pw_pool, &
3237 pw_pools=pw_pools)
3238 CALL auxbas_pw_pool%create_pw(pw=rho_elec_rspace)
3239
3240 IF (dft_control%nspins > 1) THEN
3241
3242 ! add Pout and Pz
3243 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3244 CALL pw_axpy(ec_env%rhoout_r(2), rho_elec_rspace)
3245
3246 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3247 CALL pw_axpy(ec_env%rhoz_r(2), rho_elec_rspace)
3248 ELSE
3249
3250 ! add Pout and Pz
3251 CALL pw_copy(ec_env%rhoout_r(1), rho_elec_rspace)
3252 CALL pw_axpy(ec_env%rhoz_r(1), rho_elec_rspace)
3253 END IF ! nspins
3254
3255 IF (should_print_voro /= 0) THEN
3256 CALL section_vals_val_get(print_key_voro, "OUTPUT_TEXT", l_val=voro_print_txt)
3257 IF (voro_print_txt) THEN
3258 append_voro = section_get_lval(ec_section, "PRINT%VORONOI%APPEND")
3259 my_pos_voro = "REWIND"
3260 IF (append_voro) THEN
3261 my_pos_voro = "APPEND"
3262 END IF
3263 unit_nr_voro = cp_print_key_unit_nr(logger, ec_section, "PRINT%VORONOI", extension=".voronoi", &
3264 file_position=my_pos_voro, log_filename=.false.)
3265 ELSE
3266 unit_nr_voro = 0
3267 END IF
3268 ELSE
3269 unit_nr_voro = 0
3270 END IF
3271
3272 CALL entry_voronoi_or_bqb(should_print_voro, should_print_bqb, print_key_voro, print_key_bqb, &
3273 unit_nr_voro, qs_env, rho_elec_rspace)
3274
3275 CALL auxbas_pw_pool%give_back_pw(rho_elec_rspace)
3276
3277 IF (unit_nr_voro > 0) THEN
3278 CALL cp_print_key_finished_output(unit_nr_voro, logger, ec_section, "PRINT%VORONOI")
3279 END IF
3280
3281 END IF
3282
3283 CALL timestop(handle)
3284
3285 END SUBROUTINE ec_properties
3286
3287! **************************************************************************************************
3288!> \brief Solve the Harris functional by linear scaling density purification scheme,
3289!> instead of the diagonalization performed in ec_diag_solver
3290!>
3291!> \param qs_env ...
3292!> \param matrix_ks Harris Kohn-Sham matrix
3293!> \param matrix_s Overlap matrix in Harris functional basis
3294!> \par History
3295!> 09.2020 created
3296!> \author F.Belleflamme
3297! **************************************************************************************************
3298 SUBROUTINE ec_ls_init(qs_env, matrix_ks, matrix_s)
3299 TYPE(qs_environment_type), POINTER :: qs_env
3300 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
3301
3302 CHARACTER(len=*), PARAMETER :: routinen = 'ec_ls_init'
3303
3304 INTEGER :: handle, ispin, nspins
3305 TYPE(dft_control_type), POINTER :: dft_control
3306 TYPE(energy_correction_type), POINTER :: ec_env
3307 TYPE(ls_scf_env_type), POINTER :: ls_env
3308
3309 CALL timeset(routinen, handle)
3310
3311 CALL get_qs_env(qs_env=qs_env, &
3312 dft_control=dft_control, &
3313 ec_env=ec_env)
3314 nspins = dft_control%nspins
3315 ls_env => ec_env%ls_env
3316
3317 ! create the matrix template for use in the ls procedures
3318 CALL matrix_ls_create(matrix_ls=ls_env%matrix_s, matrix_qs=matrix_s(1, 1)%matrix, &
3319 ls_mstruct=ls_env%ls_mstruct)
3320
3321 IF (ALLOCATED(ls_env%matrix_p)) THEN
3322 DO ispin = 1, SIZE(ls_env%matrix_p)
3323 CALL dbcsr_release(ls_env%matrix_p(ispin))
3324 END DO
3325 ELSE
3326 ALLOCATE (ls_env%matrix_p(nspins))
3327 END IF
3328
3329 DO ispin = 1, nspins
3330 CALL dbcsr_create(ls_env%matrix_p(ispin), template=ls_env%matrix_s, &
3331 matrix_type=dbcsr_type_no_symmetry)
3332 END DO
3333
3334 ALLOCATE (ls_env%matrix_ks(nspins))
3335 DO ispin = 1, nspins
3336 CALL dbcsr_create(ls_env%matrix_ks(ispin), template=ls_env%matrix_s, &
3337 matrix_type=dbcsr_type_no_symmetry)
3338 END DO
3339
3340 ! Set up S matrix and needed functions of S
3341 CALL ls_scf_init_matrix_s(matrix_s(1, 1)%matrix, ls_env)
3342
3343 ! Bring KS matrix from QS to LS form
3344 ! EC KS-matrix already calculated
3345 DO ispin = 1, nspins
3346 CALL matrix_qs_to_ls(matrix_ls=ls_env%matrix_ks(ispin), &
3347 matrix_qs=matrix_ks(ispin, 1)%matrix, &
3348 ls_mstruct=ls_env%ls_mstruct, &
3349 covariant=.true.)
3350 END DO
3351
3352 CALL timestop(handle)
3353
3354 END SUBROUTINE ec_ls_init
3355
3356! **************************************************************************************************
3357!> \brief Solve the Harris functional by linear scaling density purification scheme,
3358!> instead of the diagonalization performed in ec_diag_solver
3359!>
3360!> \param qs_env ...
3361!> \param matrix_p Harris dentiy matrix, calculated here
3362!> \param matrix_w Harris energy weighted density matrix, calculated here
3363!> \param ec_ls_method which purification scheme should be used
3364!> \par History
3365!> 12.2019 created [JGH]
3366!> 08.2020 refactoring [fbelle]
3367!> \author Fabian Belleflamme
3368! **************************************************************************************************
3369
3370 SUBROUTINE ec_ls_solver(qs_env, matrix_p, matrix_w, ec_ls_method)
3371 TYPE(qs_environment_type), POINTER :: qs_env
3372 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
3373 INTEGER, INTENT(IN) :: ec_ls_method
3374
3375 CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_ls_solver'
3376
3377 INTEGER :: handle, ispin, nelectron_spin_real, &
3378 nspins
3379 INTEGER, DIMENSION(2) :: nmo
3380 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: wmat
3381 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_ks_deviation
3382 TYPE(dft_control_type), POINTER :: dft_control
3383 TYPE(energy_correction_type), POINTER :: ec_env
3384 TYPE(ls_scf_env_type), POINTER :: ls_env
3385 TYPE(mp_para_env_type), POINTER :: para_env
3386
3387 CALL timeset(routinen, handle)
3388
3389 NULLIFY (para_env)
3390 CALL get_qs_env(qs_env, &
3391 dft_control=dft_control, &
3392 para_env=para_env)
3393 nspins = dft_control%nspins
3394 ec_env => qs_env%ec_env
3395 ls_env => ec_env%ls_env
3396
3397 nmo = 0
3398 CALL get_qs_env(qs_env=qs_env, nelectron_spin=nmo)
3399 IF (nspins == 1) nmo(1) = nmo(1)/2
3400 ls_env%homo_spin(:) = 0.0_dp
3401 ls_env%lumo_spin(:) = 0.0_dp
3402
3403 ALLOCATE (matrix_ks_deviation(nspins))
3404 DO ispin = 1, nspins
3405 CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_env%matrix_ks(ispin))
3406 CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
3407 END DO
3408
3409 ! F = S^-1/2 * F * S^-1/2
3410 IF (ls_env%has_s_preconditioner) THEN
3411 DO ispin = 1, nspins
3412 CALL apply_matrix_preconditioner(ls_env%matrix_ks(ispin), "forward", &
3413 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3414
3415 CALL dbcsr_filter(ls_env%matrix_ks(ispin), ls_env%eps_filter)
3416 END DO
3417 END IF
3418
3419 DO ispin = 1, nspins
3420 nelectron_spin_real = ls_env%nelectron_spin(ispin)
3421 IF (ls_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
3422
3423 SELECT CASE (ec_ls_method)
3424 CASE (ec_matrix_sign)
3425 CALL density_matrix_sign(ls_env%matrix_p(ispin), &
3426 ls_env%mu_spin(ispin), &
3427 ls_env%fixed_mu, &
3428 ls_env%sign_method, &
3429 ls_env%sign_order, &
3430 ls_env%matrix_ks(ispin), &
3431 ls_env%matrix_s, &
3432 ls_env%matrix_s_inv, &
3433 nelectron_spin_real, &
3434 ec_env%eps_default)
3435
3436 CASE (ec_matrix_trs4)
3437 CALL density_matrix_trs4( &
3438 ls_env%matrix_p(ispin), &
3439 ls_env%matrix_ks(ispin), &
3440 ls_env%matrix_s_sqrt_inv, &
3441 nelectron_spin_real, &
3442 ec_env%eps_default, &
3443 ls_env%homo_spin(ispin), &
3444 ls_env%lumo_spin(ispin), &
3445 ls_env%mu_spin(ispin), &
3446 matrix_ks_deviation=matrix_ks_deviation(ispin), &
3447 dynamic_threshold=ls_env%dynamic_threshold, &
3448 eps_lanczos=ls_env%eps_lanczos, &
3449 max_iter_lanczos=ls_env%max_iter_lanczos)
3450
3451 CASE (ec_matrix_tc2)
3452 CALL density_matrix_tc2( &
3453 ls_env%matrix_p(ispin), &
3454 ls_env%matrix_ks(ispin), &
3455 ls_env%matrix_s_sqrt_inv, &
3456 nelectron_spin_real, &
3457 ec_env%eps_default, &
3458 ls_env%homo_spin(ispin), &
3459 ls_env%lumo_spin(ispin), &
3460 non_monotonic=ls_env%non_monotonic, &
3461 eps_lanczos=ls_env%eps_lanczos, &
3462 max_iter_lanczos=ls_env%max_iter_lanczos)
3463
3464 END SELECT
3465
3466 END DO
3467
3468 ! de-orthonormalize
3469 IF (ls_env%has_s_preconditioner) THEN
3470 DO ispin = 1, nspins
3471 ! P = S^-1/2 * P_tilde * S^-1/2 (forward)
3472 CALL apply_matrix_preconditioner(ls_env%matrix_p(ispin), "forward", &
3473 ls_env%matrix_bs_sqrt, ls_env%matrix_bs_sqrt_inv)
3474
3475 CALL dbcsr_filter(ls_env%matrix_p(ispin), ls_env%eps_filter)
3476 END DO
3477 END IF
3478
3479 ! Closed-shell
3480 IF (nspins == 1) CALL dbcsr_scale(ls_env%matrix_p(1), 2.0_dp)
3481
3482 IF (ls_env%report_all_sparsities) CALL post_scf_sparsities(ls_env)
3483
3484 ! ls_scf_dm_to_ks
3485 ! Density matrix from LS to EC
3486 DO ispin = 1, nspins
3487 CALL matrix_ls_to_qs(matrix_qs=matrix_p(ispin, 1)%matrix, &
3488 matrix_ls=ls_env%matrix_p(ispin), &
3489 ls_mstruct=ls_env%ls_mstruct, &
3490 covariant=.false.)
3491 END DO
3492
3493 wmat => matrix_w(:, 1)
3494 CALL calculate_w_matrix_ls(wmat, ec_env%ls_env)
3495
3496 ! clean up
3497 CALL dbcsr_release(ls_env%matrix_s)
3498 IF (ls_env%has_s_preconditioner) THEN
3499 CALL dbcsr_release(ls_env%matrix_bs_sqrt)
3500 CALL dbcsr_release(ls_env%matrix_bs_sqrt_inv)
3501 END IF
3502 IF (ls_env%needs_s_inv) THEN
3503 CALL dbcsr_release(ls_env%matrix_s_inv)
3504 END IF
3505 IF (ls_env%use_s_sqrt) THEN
3506 CALL dbcsr_release(ls_env%matrix_s_sqrt)
3507 CALL dbcsr_release(ls_env%matrix_s_sqrt_inv)
3508 END IF
3509
3510 DO ispin = 1, SIZE(ls_env%matrix_ks)
3511 CALL dbcsr_release(ls_env%matrix_ks(ispin))
3512 END DO
3513 DEALLOCATE (ls_env%matrix_ks)
3514
3515 DO ispin = 1, nspins
3516 CALL dbcsr_release(matrix_ks_deviation(ispin))
3517 END DO
3518 DEALLOCATE (matrix_ks_deviation)
3519
3520 CALL timestop(handle)
3521
3522 END SUBROUTINE ec_ls_solver
3523
3524! **************************************************************************************************
3525!> \brief Use OT-diagonalziation to obtain density matrix from Harris Kohn-Sham matrix
3526!> Initial guess of density matrix is either the atomic block initial guess from SCF
3527!> or the ground-state density matrix. The latter only works if the same basis is used
3528!>
3529!> \param qs_env ...
3530!> \param ec_env ...
3531!> \param matrix_ks Harris Kohn-Sham matrix
3532!> \param matrix_s Overlap matrix in Harris functional basis
3533!> \param matrix_p Harris dentiy matrix, calculated here
3534!> \param matrix_w Harris energy weighted density matrix, calculated here
3535!>
3536!> \par History
3537!> 09.2020 created
3538!> \author F.Belleflamme
3539! **************************************************************************************************
3540 SUBROUTINE ec_ot_diag_solver(qs_env, ec_env, matrix_ks, matrix_s, matrix_p, matrix_w)
3541 TYPE(qs_environment_type), POINTER :: qs_env
3542 TYPE(energy_correction_type), POINTER :: ec_env
3543 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
3544 POINTER :: matrix_ks, matrix_s
3545 TYPE(dbcsr_p_type), DIMENSION(:, :), &
3546 INTENT(INOUT), POINTER :: matrix_p, matrix_w
3547
3548 CHARACTER(len=*), PARAMETER :: routinen = 'ec_ot_diag_solver'
3549
3550 INTEGER :: handle, homo, ikind, iounit, ispin, &
3551 max_iter, nao, nkind, nmo, nspins
3552 INTEGER, DIMENSION(2) :: nelectron_spin
3553 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
3554 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3555 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3556 TYPE(cp_fm_type) :: sv
3557 TYPE(cp_fm_type), POINTER :: mo_coeff
3558 TYPE(cp_logger_type), POINTER :: logger
3559 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_rmpv
3560 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
3561 TYPE(dft_control_type), POINTER :: dft_control
3562 TYPE(gto_basis_set_type), POINTER :: basis_set, harris_basis
3563 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3564 TYPE(mp_para_env_type), POINTER :: para_env
3565 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3566 TYPE(preconditioner_type), POINTER :: local_preconditioner
3567 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3568 TYPE(qs_kind_type), POINTER :: qs_kind
3569 TYPE(qs_rho_type), POINTER :: rho
3570
3571 CALL timeset(routinen, handle)
3572
3573 logger => cp_get_default_logger()
3574 iounit = cp_logger_get_default_unit_nr(logger)
3575
3576 cpassert(ASSOCIATED(qs_env))
3577 cpassert(ASSOCIATED(ec_env))
3578 cpassert(ASSOCIATED(matrix_ks))
3579 cpassert(ASSOCIATED(matrix_s))
3580 cpassert(ASSOCIATED(matrix_p))
3581 cpassert(ASSOCIATED(matrix_w))
3582
3583 CALL get_qs_env(qs_env=qs_env, &
3584 atomic_kind_set=atomic_kind_set, &
3585 blacs_env=blacs_env, &
3586 dft_control=dft_control, &
3587 nelectron_spin=nelectron_spin, &
3588 para_env=para_env, &
3589 particle_set=particle_set, &
3590 qs_kind_set=qs_kind_set)
3591 nspins = dft_control%nspins
3592
3593 ! Maximum number of OT iterations for diagonalization
3594 max_iter = 200
3595
3596 ! If linear scaling, need to allocate and init MO set
3597 ! set it to qs_env%mos
3598 IF (dft_control%qs_control%do_ls_scf) THEN
3599 CALL ec_mos_init(qs_env, matrix_s(1, 1)%matrix)
3600 END IF
3601
3602 CALL get_qs_env(qs_env, mos=mos)
3603
3604 ! Inital guess to use
3605 NULLIFY (p_rmpv)
3606
3607 ! Using ether ground-state DM or ATOMIC GUESS requires
3608 ! Harris functional to use the same basis set
3609 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind)
3610 CALL uppercase(ec_env%basis)
3611 ! Harris basis only differs from ground-state basis if explicitly added
3612 ! thus only two cases that need to be tested
3613 ! 1) explicit Harris basis present?
3614 IF (ec_env%basis == "HARRIS") THEN
3615 DO ikind = 1, nkind
3616 qs_kind => qs_kind_set(ikind)
3617 ! Basis sets of ground-state
3618 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB")
3619 ! Basis sets of energy correction
3620 CALL get_qs_kind(qs_kind=qs_kind, basis_set=harris_basis, basis_type="HARRIS")
3621
3622 IF (basis_set%name .NE. harris_basis%name) THEN
3623 cpabort("OT-Diag initial guess: Harris and ground state need to use the same basis")
3624 END IF
3625 END DO
3626 END IF
3627 ! 2) Harris uses MAOs
3628 IF (ec_env%mao) THEN
3629 cpabort("OT-Diag initial guess: not implemented for MAOs")
3630 END IF
3631
3632 ! Initital guess obtained for OT Diagonalization
3633 SELECT CASE (ec_env%ec_initial_guess)
3634 CASE (ec_ot_atomic)
3635
3636 p_rmpv => matrix_p(:, 1)
3637
3638 CALL calculate_atomic_block_dm(p_rmpv, matrix_s(1, 1)%matrix, atomic_kind_set, qs_kind_set, &
3639 nspins, nelectron_spin, iounit, para_env)
3640
3641 CASE (ec_ot_gs)
3642
3643 CALL get_qs_env(qs_env, rho=rho)
3644 CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3645 p_rmpv => rho_ao(:, 1)
3646
3647 CASE DEFAULT
3648 cpabort("Unknown inital guess for OT-Diagonalization (Harris functional)")
3649 END SELECT
3650
3651 DO ispin = 1, nspins
3652 CALL get_mo_set(mo_set=mos(ispin), &
3653 mo_coeff=mo_coeff, &
3654 nmo=nmo, &
3655 nao=nao, &
3656 homo=homo)
3657
3658 ! Calculate first MOs
3659 CALL cp_fm_set_all(mo_coeff, 0.0_dp)
3660 CALL cp_fm_init_random(mo_coeff, nmo)
3661
3662 CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
3663 ! multiply times PS
3664 ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
3665 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff, sv, nmo)
3666 CALL cp_dbcsr_sm_fm_multiply(p_rmpv(ispin)%matrix, sv, mo_coeff, homo)
3667 CALL cp_fm_release(sv)
3668 ! and ortho the result
3669 ! If DFBT or SE, then needs has_unit_metrix option
3670 CALL make_basis_sm(mo_coeff, nmo, matrix_s(1, 1)%matrix)
3671 END DO
3672
3673 ! Preconditioner
3674 NULLIFY (local_preconditioner)
3675 ALLOCATE (local_preconditioner)
3676 CALL init_preconditioner(local_preconditioner, para_env=para_env, &
3677 blacs_env=blacs_env)
3678 DO ispin = 1, nspins
3679 CALL make_preconditioner(local_preconditioner, &
3680 precon_type=ot_precond_full_single_inverse, &
3681 solver_type=ot_precond_solver_default, &
3682 matrix_h=matrix_ks(ispin, 1)%matrix, &
3683 matrix_s=matrix_s(ispin, 1)%matrix, &
3684 convert_precond_to_dbcsr=.true., &
3685 mo_set=mos(ispin), energy_gap=0.2_dp)
3686
3687 CALL get_mo_set(mos(ispin), &
3688 mo_coeff=mo_coeff, &
3689 eigenvalues=eigenvalues, &
3690 nmo=nmo, &
3691 homo=homo)
3692 CALL ot_eigensolver(matrix_h=matrix_ks(ispin, 1)%matrix, &
3693 matrix_s=matrix_s(1, 1)%matrix, &
3694 matrix_c_fm=mo_coeff, &
3695 preconditioner=local_preconditioner, &
3696 eps_gradient=ec_env%eps_default, &
3697 iter_max=max_iter, &
3698 silent=.false.)
3699 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin, 1)%matrix, &
3700 evals_arg=eigenvalues, do_rotation=.true.)
3701
3702 ! Deallocate preconditioner
3703 CALL destroy_preconditioner(local_preconditioner)
3704 DEALLOCATE (local_preconditioner)
3705
3706 !fm->dbcsr
3707 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
3708 mos(ispin)%mo_coeff_b)
3709 END DO
3710
3711 ! Calculate density matrix from MOs
3712 DO ispin = 1, nspins
3713 CALL calculate_density_matrix(mos(ispin), matrix_p(ispin, 1)%matrix)
3714
3715 CALL calculate_w_matrix(mos(ispin), matrix_w(ispin, 1)%matrix)
3716 END DO
3717
3718 ! Get rid of MO environment again
3719 IF (dft_control%qs_control%do_ls_scf) THEN
3720 DO ispin = 1, nspins
3721 CALL deallocate_mo_set(mos(ispin))
3722 END DO
3723 IF (ASSOCIATED(qs_env%mos)) THEN
3724 DO ispin = 1, SIZE(qs_env%mos)
3725 CALL deallocate_mo_set(qs_env%mos(ispin))
3726 END DO
3727 DEALLOCATE (qs_env%mos)
3728 END IF
3729 END IF
3730
3731 CALL timestop(handle)
3732
3733 END SUBROUTINE ec_ot_diag_solver
3734
3735END MODULE energy_corrections
3736
Contains ADMM methods which only require the density matrix.
subroutine, public admm_dm_calc_rho_aux(qs_env)
Entry methods: Calculates auxiliary density matrix from primary one.
Contains ADMM methods which require molecular orbitals.
subroutine, public admm_mo_calc_rho_aux(qs_env)
...
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public belleflamme2023
Handles all functions related to the CELL.
Definition cell_types.F:15
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
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_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_triangular_invert(matrix_a, uplo_tr)
inverts a triangular matrix
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:236
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
subroutine, public cp_fm_to_fm_triangular(msource, mtarget, uplo)
copy just a triangular matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
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
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
set of type/routines to handle the storage of results in force_envs
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
set of type/routines to handle the storage of results in force_envs
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1178
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
lower level routines for linear scaling SCF
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged, iounit)
compute the density matrix using a trace-resetting algorithm
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
subroutine, public apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv)
apply a preconditioner either forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs)) backward (rest...
subroutine, public density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, matrix_s_sqrt_inv)
compute the density matrix with a trace that is close to nelectron. take a mu as input,...
subroutine, public density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, non_monotonic, eps_lanczos, max_iter_lanczos, iounit)
compute the density matrix using a non monotonic trace conserving algorithm based on SIAM DOI....
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
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_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 ...
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Routines for a linear scaling quickstep SCF run based on the density matrix.
Definition dm_ls_scf.F:15
subroutine, public post_scf_sparsities(ls_scf_env)
Report on the sparsities of various interesting matrices.
Definition dm_ls_scf.F:1066
subroutine, public calculate_w_matrix_ls(matrix_w, ls_scf_env)
Compute matrix_w as needed for the forces.
Definition dm_ls_scf.F:1223
Calculates the energy contribution and the mo_derivative of a static electric field (nonperiodic)
subroutine, public ec_efield_local_operator(qs_env, ec_env, calculate_forces)
...
subroutine, public ec_efield_integrals(qs_env, ec_env, rpoint)
...
Types needed for a for a Energy Correction.
Routines for an external energy correction on top of a Kohn-Sham calculation.
Definition ec_external.F:14
subroutine, public ec_ext_energy(qs_env, ec_env, calculate_forces)
External energy method.
Definition ec_external.F:91
Routines used for Harris functional Kohn-Sham calculation.
Definition ec_methods.F:15
subroutine, public ec_mos_init(qs_env, matrix_s)
Allocate and initiate molecular orbitals environment.
Definition ec_methods.F:162
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
Definition ec_methods.F:88
Routines for an energy correction on top of a Kohn-Sham calculation.
subroutine, public energy_correction(qs_env, ec_init, calculate_forces)
Energy Correction to a Kohn-Sham simulation Available energy corrections: (1) Harris energy functiona...
Definition of the atomic potential types.
Routines to calculate EXX in RPA and energy correction methods.
Definition hfx_exx.F:16
subroutine, public calculate_exx(qs_env, unit_nr, hfx_sections, x_data, do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, e_ex_from_gw, e_admm_from_gw, t3)
...
Definition hfx_exx.F:106
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
Definition hfx_exx.F:325
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public ec_functional_harris
integer, parameter, public ec_functional_dc
integer, parameter, public vdw_pairpot_dftd3
integer, parameter, public ec_ot_atomic
integer, parameter, public ec_ot_diag
integer, parameter, public ec_ot_gs
integer, parameter, public ot_precond_solver_default
integer, parameter, public ot_precond_full_single_inverse
integer, parameter, public ec_matrix_tc2
integer, parameter, public ec_diagonalization
integer, parameter, public ec_matrix_trs4
integer, parameter, public ec_functional_ext
integer, parameter, public ec_matrix_sign
integer, parameter, public xc_vdw_fun_pairpot
integer, parameter, public vdw_pairpot_dftd3bj
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
integer function, public section_get_ival(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
logical function, public section_get_lval(section_vals, keyword_name)
...
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
integer, parameter, public default_path_length
Definition kinds.F:58
Calculate MAO's and analyze wavefunctions.
Definition mao_basis.F:15
subroutine, public mao_generate_basis(qs_env, mao_coef, ref_basis_set, pmat_external, smat_external, molecular, max_iter, eps_grad, nmao_external, eps1_mao, iolevel, unit_nr)
...
Definition mao_basis.F:75
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
Interface to the message passing library MPI.
Define the data structure for the molecule information.
Calculates the moment integrals <a|r^m|b>
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Define the data structure for the particle information.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public pascal
Definition physcon.F:174
real(kind=dp), parameter, public bohr
Definition physcon.F:147
real(kind=dp), parameter, public debye
Definition physcon.F:201
types of preconditioners
subroutine, public init_preconditioner(preconditioner_env, para_env, blacs_env)
...
subroutine, public destroy_preconditioner(preconditioner_env)
...
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, matrix_t, mo_set, energy_gap, convert_precond_to_dbcsr, chol_type)
...
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
functions related to the poisson solver on regular grids
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.
subroutine, public calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, e_overlap_core, atecc)
Calculate the overlap energy of the core charge distribution.
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
subroutine, public kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, matrix_name, calculate_forces, nderivative, sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
Calculate kinetic energy matrix and possible relativistic correction.
subroutine, public core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
...
collects routines that calculate density matrices
Calculation of dispersion using pair potentials.
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition of disperson types for DFT calculations.
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, sab_cneo, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Calculation of kinetic energy matrix and forces.
Definition qs_kinetic.F:15
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition qs_kinetic.F:101
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
collects routines that perform operations directly related to MOs
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:560
Define the neighbor list data types and the corresponding functionality.
Generate the atomic neighbor lists.
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
subroutine, public build_oce_matrices(intac, calculate_forces, nder, qs_kind_set, particle_set, sap_oce, eps_fit)
Set up the sparse matrix for the coefficients of one center expansions This routine uses the same log...
subroutine, public allocate_oce_set(oce_set, nkind)
Allocate and initialize the matrix set of oce coefficients.
subroutine, public create_oce_set(oce_set)
...
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
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...
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:98
Calculate the CPKS equation and the resulting forces.
subroutine, public response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
...
subroutine, public response_calculation(qs_env, ec_env)
Initializes solver of linear response equation for energy correction.
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's)
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
generate the tasks lists used by collocate and integrate routines
subroutine, public generate_qs_task_list(ks_env, task_list, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
types for task lists
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
The module to read/write TREX IO files for interfacing CP2K with other programs.
subroutine, public write_trexio(qs_env, trexio_section, energy_derivative)
Write a trexio file.
subroutine, public write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
Print stress tensor to output file.
subroutine, public write_stress_tensor_components(virial, iw, cell, unit_string)
...
pure real(kind=dp) function, public one_third_sum_diag(a)
...
subroutine, public zero_virial(virial, reset)
...
subroutine, public symmetrize_virial(virial)
Symmetrize the virial components.
Interface for Voronoi Integration and output of BQB files.
subroutine, public entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
Does a Voronoi integration of density or stores the density to compressed BQB format.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
contains arbitrary information which need to be stored
structure to store local (to a processor) ordered lists of integers.
distributes pairs on a 2d grid of processors
Contains information on the energy correction functional for KG.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
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 ...
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.