1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ***************************************************************************** 7!> \brief Utility subroutine for qs energy calculation 8!> \par History 9!> none 10!> \author MK (29.10.2002) 11! ***************************************************************************** 12MODULE qs_energy_matrix_w 13 USE cp_control_types, ONLY: dft_control_type 14 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 15 cp_fm_struct_release,& 16 cp_fm_struct_type 17 USE cp_fm_types, ONLY: cp_fm_create,& 18 cp_fm_p_type,& 19 cp_fm_release,& 20 cp_fm_type 21 USE dbcsr_api, ONLY: dbcsr_p_type,& 22 dbcsr_set 23 USE kinds, ONLY: dp 24 USE kpoint_methods, ONLY: kpoint_density_matrices,& 25 kpoint_density_transform 26 USE kpoint_types, ONLY: kpoint_type 27 USE qs_environment_types, ONLY: get_qs_env,& 28 qs_environment_type 29 USE qs_ks_methods, ONLY: calculate_w_matrix,& 30 calculate_w_matrix_ot 31 USE qs_mo_types, ONLY: get_mo_set,& 32 mo_set_p_type,& 33 mo_set_type 34 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 35 USE qs_rho_types, ONLY: qs_rho_get,& 36 qs_rho_type 37 USE scf_control_types, ONLY: scf_control_type 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 42 PRIVATE 43 44! *** Global parameters *** 45 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_matrix_w' 47 48 PUBLIC :: qs_energies_compute_w 49 50CONTAINS 51 52! ***************************************************************************** 53!> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w 54!> into separate subroutine 55!> \param qs_env ... 56!> \param calc_forces ... 57!> \par History 58!> 05.2013 created [Florian Schiffmann] 59! ************************************************************************************************** 60 61 SUBROUTINE qs_energies_compute_w(qs_env, calc_forces) 62 TYPE(qs_environment_type), POINTER :: qs_env 63 LOGICAL, INTENT(IN) :: calc_forces 64 65 CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_compute_w', & 66 routineP = moduleN//':'//routineN 67 68 INTEGER :: handle, is, ispin, nao, nspin 69 LOGICAL :: do_kpoints, has_unit_metric 70 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fmwork 71 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct 72 TYPE(cp_fm_type), POINTER :: mo_coeff 73 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, & 74 matrix_w_mp2, mo_derivs, rho_ao 75 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrix_w_kp 76 TYPE(dft_control_type), POINTER :: dft_control 77 TYPE(kpoint_type), POINTER :: kpoints 78 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 79 TYPE(mo_set_type), POINTER :: mo_set 80 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 81 POINTER :: sab_nl 82 TYPE(qs_rho_type), POINTER :: rho 83 TYPE(scf_control_type), POINTER :: scf_control 84 85 CALL timeset(routineN, handle) 86 87 ! if calculate forces, time to compute the w matrix 88 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) 89 90 IF (calc_forces .AND. .NOT. has_unit_metric) THEN 91 CALL get_qs_env(qs_env, do_kpoints=do_kpoints) 92 93 IF (do_kpoints) THEN 94 95 CALL get_qs_env(qs_env, & 96 matrix_w_kp=matrix_w_kp, & 97 matrix_s_kp=matrix_s_kp, & 98 sab_orb=sab_nl, & 99 mos=mos, & 100 kpoints=kpoints) 101 102 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao) 103 CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, & 104 template_fmstruct=mo_coeff%matrix_struct) 105 106 ALLOCATE (fmwork(2)) 107 DO is = 1, SIZE(fmwork) 108 NULLIFY (fmwork(is)%matrix) 109 CALL cp_fm_create(fmwork(is)%matrix, matrix_struct=ao_ao_fmstruct) 110 END DO 111 CALL cp_fm_struct_release(ao_ao_fmstruct) 112 113 ! energy weighted density matrices in k-space 114 CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.) 115 ! energy weighted density matrices in real space 116 CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., & 117 matrix_s_kp(1, 1)%matrix, sab_nl, fmwork) 118 119 DO is = 1, SIZE(fmwork) 120 CALL cp_fm_release(fmwork(is)%matrix) 121 END DO 122 DEALLOCATE (fmwork) 123 124 ELSE 125 126 NULLIFY (dft_control, rho_ao) 127 CALL get_qs_env(qs_env, & 128 matrix_w=matrix_w, & 129 matrix_ks=matrix_ks, & 130 matrix_s=matrix_s, & 131 matrix_w_mp2=matrix_w_mp2, & 132 mo_derivs=mo_derivs, & 133 scf_control=scf_control, & 134 mos=mos, & 135 rho=rho, & 136 dft_control=dft_control) 137 138 CALL qs_rho_get(rho, rho_ao=rho_ao) 139 140 nspin = SIZE(mos) 141 DO ispin = 1, nspin 142 mo_set => mos(ispin)%mo_set 143 IF (dft_control%roks) THEN 144 IF (scf_control%use_ot) THEN 145 IF (ispin > 1) THEN 146 ! not very elegant, indeed ... 147 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp) 148 ELSE 149 CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & 150 matrix_w(ispin)%matrix, matrix_s(1)%matrix) 151 END IF 152 ELSE 153 CALL calculate_w_matrix(mo_set=mo_set, & 154 matrix_ks=matrix_ks(ispin)%matrix, & 155 matrix_p=rho_ao(ispin)%matrix, & 156 matrix_w=matrix_w(ispin)%matrix) 157 END IF 158 ELSE 159 IF (scf_control%use_ot) THEN 160 CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & 161 matrix_w(ispin)%matrix, matrix_s(1)%matrix) 162 ELSE 163 CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix) 164 END IF 165 END IF 166 END DO 167 168 END IF 169 170 END IF 171 CALL timestop(handle) 172 173 END SUBROUTINE qs_energies_compute_w 174 175END MODULE qs_energy_matrix_w 176