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