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