Line data Source code
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 external energy correction on top of a Kohn-Sham calculation
10 : !> \par History
11 : !> 10.2024 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE ec_external
15 : USE cell_types, ONLY: cell_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
18 : dbcsr_copy,&
19 : dbcsr_create,&
20 : dbcsr_p_type,&
21 : dbcsr_set,&
22 : dbcsr_type,&
23 : dbcsr_type_symmetric
24 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
26 : cp_dbcsr_sm_fm_multiply,&
27 : dbcsr_allocate_matrix_set,&
28 : dbcsr_deallocate_matrix_set
29 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30 : cp_fm_struct_release,&
31 : cp_fm_struct_type
32 : USE cp_fm_types, ONLY: cp_fm_create,&
33 : cp_fm_get_diag,&
34 : cp_fm_get_info,&
35 : cp_fm_maxabsval,&
36 : cp_fm_release,&
37 : cp_fm_set_all,&
38 : cp_fm_to_fm,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_get_default_logger,&
41 : cp_logger_get_default_unit_nr,&
42 : cp_logger_type
43 : USE ec_env_types, ONLY: energy_correction_type
44 : USE kinds, ONLY: default_string_length,&
45 : dp
46 : USE mathlib, ONLY: det_3x3
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE parallel_gemm_api, ONLY: parallel_gemm
49 : USE physcon, ONLY: pascal
50 : USE qs_core_energies, ONLY: calculate_ptrace
51 : USE qs_core_matrices, ONLY: core_matrices,&
52 : kinetic_energy_matrix
53 : USE qs_environment_types, ONLY: get_qs_env,&
54 : qs_environment_type
55 : USE qs_force_types, ONLY: qs_force_type
56 : USE qs_ks_types, ONLY: qs_ks_env_type
57 : USE qs_mo_types, ONLY: deallocate_mo_set,&
58 : get_mo_set,&
59 : mo_set_type
60 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
61 : USE qs_overlap, ONLY: build_overlap_matrix
62 : USE qs_rho_types, ONLY: qs_rho_get,&
63 : qs_rho_type
64 : USE trexio_utils, ONLY: read_trexio
65 : USE virial_methods, ONLY: one_third_sum_diag
66 : USE virial_types, ONLY: virial_type
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : ! *** Global parameters ***
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_external'
76 :
77 : PUBLIC :: ec_ext_energy
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief External energy method
83 : !> \param qs_env ...
84 : !> \param ec_env ...
85 : !> \param calculate_forces ...
86 : !> \par History
87 : !> 10.2024 created
88 : !> \author JGH
89 : ! **************************************************************************************************
90 44 : SUBROUTINE ec_ext_energy(qs_env, ec_env, calculate_forces)
91 : TYPE(qs_environment_type), POINTER :: qs_env
92 : TYPE(energy_correction_type), POINTER :: ec_env
93 : LOGICAL, INTENT(IN) :: calculate_forces
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'ec_ext_energy'
96 :
97 : INTEGER :: handle, ispin, nocc, nspins, unit_nr
98 : REAL(KIND=dp) :: focc
99 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
100 44 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ, mo_ref
101 : TYPE(cp_fm_type), POINTER :: mo_coeff
102 : TYPE(cp_logger_type), POINTER :: logger
103 44 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
104 : TYPE(dft_control_type), POINTER :: dft_control
105 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
106 :
107 44 : CALL timeset(routineN, handle)
108 :
109 44 : CALL get_qs_env(qs_env, dft_control=dft_control)
110 44 : nspins = dft_control%nspins
111 :
112 44 : logger => cp_get_default_logger()
113 44 : IF (logger%para_env%is_source()) THEN
114 22 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
115 : ELSE
116 22 : unit_nr = -1
117 : END IF
118 :
119 44 : ec_env%etotal = 0.0_dp
120 44 : ec_env%ex = 0.0_dp
121 44 : ec_env%exc = 0.0_dp
122 44 : ec_env%ehartree = 0.0_dp
123 :
124 44 : CALL get_qs_env(qs_env, mos=mos)
125 352 : ALLOCATE (cpmos(nspins), mo_ref(nspins), mo_occ(nspins))
126 :
127 44 : CALL cp_fm_release(ec_env%mo_occ)
128 44 : CALL cp_fm_release(ec_env%cpmos)
129 44 : IF (ec_env%debug_external) THEN
130 : !
131 88 : DO ispin = 1, nspins
132 44 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
133 44 : NULLIFY (fm_struct)
134 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
135 44 : template_fmstruct=mo_coeff%matrix_struct)
136 44 : CALL cp_fm_create(cpmos(ispin), fm_struct)
137 44 : CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
138 44 : CALL cp_fm_create(mo_ref(ispin), fm_struct)
139 44 : CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp)
140 44 : CALL cp_fm_create(mo_occ(ispin), fm_struct)
141 44 : CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
142 132 : CALL cp_fm_struct_release(fm_struct)
143 : END DO
144 : !
145 44 : ec_env%mo_occ => mo_ref
146 44 : CALL ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr)
147 : !
148 44 : IF (calculate_forces) THEN
149 12 : focc = 2.0_dp
150 12 : IF (nspins == 1) focc = 4.0_dp
151 24 : DO ispin = 1, nspins
152 12 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
153 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_h(1, 1)%matrix, ec_env%mo_occ(ispin), &
154 : cpmos(ispin), nocc, &
155 24 : alpha=focc, beta=0.0_dp)
156 : END DO
157 : END IF
158 44 : ec_env%cpmos => cpmos
159 : ELSE
160 0 : DO ispin = 1, nspins
161 0 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
162 0 : NULLIFY (fm_struct)
163 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
164 0 : template_fmstruct=mo_coeff%matrix_struct)
165 0 : CALL cp_fm_create(cpmos(ispin), fm_struct)
166 0 : CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
167 0 : CALL cp_fm_create(mo_occ(ispin), fm_struct)
168 0 : CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
169 0 : CALL cp_fm_create(mo_ref(ispin), fm_struct)
170 0 : CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp)
171 0 : CALL cp_fm_struct_release(fm_struct)
172 : END DO
173 :
174 : ! get external information
175 0 : CALL ec_ext_interface(qs_env, ec_env%exresp_fn, mo_occ, mo_ref, cpmos, calculate_forces, unit_nr)
176 0 : ec_env%mo_occ => mo_ref
177 0 : ec_env%cpmos => cpmos
178 : END IF
179 :
180 44 : IF (calculate_forces) THEN
181 : ! check for orbital rotations
182 12 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
183 24 : DO ispin = 1, nspins
184 : CALL align_vectors(ec_env%cpmos(ispin), ec_env%mo_occ(ispin), mo_occ(ispin), &
185 24 : matrix_s(1)%matrix, unit_nr)
186 : END DO
187 : ! set up matrices for response
188 12 : CALL ec_ext_setup(qs_env, ec_env, .TRUE., unit_nr)
189 : ! orthogonality force
190 : CALL matrix_r_forces(qs_env, ec_env%cpmos, ec_env%mo_occ, &
191 12 : ec_env%matrix_w(1, 1)%matrix, unit_nr)
192 : ELSE
193 32 : CALL ec_ext_setup(qs_env, ec_env, .FALSE., unit_nr)
194 : END IF
195 :
196 44 : CALL cp_fm_release(mo_occ)
197 :
198 44 : CALL timestop(handle)
199 :
200 44 : END SUBROUTINE ec_ext_energy
201 :
202 : ! **************************************************************************************************
203 :
204 : ! **************************************************************************************************
205 : !> \brief ...
206 : !> \param qs_env ...
207 : !> \param trexio_fn ...
208 : !> \param mo_occ ...
209 : !> \param mo_ref ...
210 : !> \param cpmos ...
211 : !> \param calculate_forces ...
212 : !> \param unit_nr ...
213 : ! **************************************************************************************************
214 0 : SUBROUTINE ec_ext_interface(qs_env, trexio_fn, mo_occ, mo_ref, cpmos, calculate_forces, unit_nr)
215 : TYPE(qs_environment_type), POINTER :: qs_env
216 : CHARACTER(LEN=*) :: trexio_fn
217 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_occ, mo_ref, cpmos
218 : LOGICAL, INTENT(IN) :: calculate_forces
219 : INTEGER, INTENT(IN) :: unit_nr
220 :
221 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_interface'
222 :
223 : INTEGER :: handle, ispin, nao, nmos, nocc(2), nspins
224 : REAL(KIND=dp) :: focc
225 : TYPE(cp_fm_type), POINTER :: mo_coeff
226 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: energy_derivative
227 : TYPE(dft_control_type), POINTER :: dft_control
228 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_trexio
229 : TYPE(mp_para_env_type), POINTER :: para_env
230 :
231 0 : CALL timeset(routineN, handle)
232 :
233 0 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, para_env=para_env)
234 :
235 0 : nspins = dft_control%nspins
236 0 : nocc = 0
237 0 : DO ispin = 1, nspins
238 0 : CALL cp_fm_get_info(mo_occ(ispin), nrow_global=nao, ncol_global=nocc(ispin))
239 : END DO
240 :
241 0 : IF (unit_nr > 0) THEN
242 0 : WRITE (unit_nr, '(T2,A)') " Read EXTERNAL Response from file: "//TRIM(trexio_fn)
243 : END IF
244 0 : ALLOCATE (mos_trexio(nspins))
245 0 : IF (calculate_forces) THEN
246 0 : NULLIFY (energy_derivative)
247 0 : CALL dbcsr_allocate_matrix_set(energy_derivative, nspins)
248 : !
249 : CALL read_trexio(qs_env, trexio_filename=trexio_fn, &
250 : mo_set_trexio=mos_trexio, &
251 0 : energy_derivative=energy_derivative)
252 : !
253 0 : focc = 2.0_dp
254 0 : IF (nspins == 1) focc = 4.0_dp
255 0 : DO ispin = 1, nspins
256 0 : CALL get_mo_set(mo_set=mos_trexio(ispin), mo_coeff=mo_coeff, homo=nmos)
257 : CALL cp_dbcsr_sm_fm_multiply(energy_derivative(ispin)%matrix, mo_coeff, &
258 0 : cpmos(ispin), ncol=nmos, alpha=focc, beta=0.0_dp)
259 : END DO
260 : !
261 0 : CALL dbcsr_deallocate_matrix_set(energy_derivative)
262 : ELSE
263 : CALL read_trexio(qs_env, trexio_filename=trexio_fn, &
264 0 : mo_set_trexio=mos_trexio)
265 : END IF
266 : !
267 0 : DO ispin = 1, nspins
268 0 : CALL get_mo_set(mo_set=mos_trexio(ispin), mo_coeff=mo_coeff, homo=nmos)
269 0 : CALL cp_fm_to_fm(mo_coeff, mo_ref(ispin), nmos)
270 0 : CALL deallocate_mo_set(mos_trexio(ispin))
271 : END DO
272 0 : DEALLOCATE (mos_trexio)
273 :
274 0 : CALL timestop(handle)
275 :
276 0 : END SUBROUTINE ec_ext_interface
277 :
278 : ! **************************************************************************************************
279 :
280 : ! **************************************************************************************************
281 : !> \brief ...
282 : !> \param qs_env ...
283 : !> \param ec_env ...
284 : !> \param calculate_forces ...
285 : !> \param unit_nr ...
286 : ! **************************************************************************************************
287 44 : SUBROUTINE ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr)
288 : TYPE(qs_environment_type), POINTER :: qs_env
289 : TYPE(energy_correction_type), POINTER :: ec_env
290 : LOGICAL, INTENT(IN) :: calculate_forces
291 : INTEGER, INTENT(IN) :: unit_nr
292 :
293 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_debug'
294 :
295 : CHARACTER(LEN=default_string_length) :: headline
296 : INTEGER :: handle, ispin, nocc, nspins
297 : TYPE(cp_fm_type), POINTER :: mo_coeff
298 44 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s
299 : TYPE(dft_control_type), POINTER :: dft_control
300 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
301 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
302 44 : POINTER :: sab_orb
303 : TYPE(qs_rho_type), POINTER :: rho
304 :
305 44 : CALL timeset(routineN, handle)
306 :
307 : CALL get_qs_env(qs_env=qs_env, &
308 : dft_control=dft_control, &
309 : sab_orb=sab_orb, &
310 : rho=rho, &
311 : matrix_s_kp=matrix_s, &
312 44 : matrix_h_kp=matrix_h)
313 :
314 44 : nspins = dft_control%nspins
315 44 : CALL get_qs_env(qs_env, mos=mos)
316 88 : DO ispin = 1, nspins
317 44 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
318 88 : CALL cp_fm_to_fm(mo_coeff, ec_env%mo_occ(ispin), nocc)
319 : END DO
320 :
321 : ! Core Hamiltonian matrix
322 44 : IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h)
323 44 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1)
324 44 : headline = "CORE HAMILTONIAN MATRIX"
325 44 : ALLOCATE (ec_env%matrix_h(1, 1)%matrix)
326 : CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=TRIM(headline), &
327 44 : template=matrix_h(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
328 44 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb)
329 44 : CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix)
330 :
331 : ! Get density matrix of reference calculation
332 44 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
333 : ! Use Core energy as model energy
334 44 : CALL calculate_ptrace(ec_env%matrix_h, matrix_p, ec_env%ex, nspins)
335 :
336 44 : IF (calculate_forces) THEN
337 : ! force of model energy
338 12 : CALL ec_debug_force(qs_env, matrix_p, unit_nr)
339 : END IF
340 :
341 44 : CALL timestop(handle)
342 :
343 44 : END SUBROUTINE ec_ext_debug
344 :
345 : ! **************************************************************************************************
346 : !> \brief ...
347 : !> \param qs_env ...
348 : !> \param matrix_p ...
349 : !> \param unit_nr ...
350 : ! **************************************************************************************************
351 12 : SUBROUTINE ec_debug_force(qs_env, matrix_p, unit_nr)
352 : TYPE(qs_environment_type), POINTER :: qs_env
353 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
354 : INTEGER, INTENT(IN) :: unit_nr
355 :
356 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force'
357 :
358 : INTEGER :: handle, iounit, nder, nimages
359 : LOGICAL :: calculate_forces, debug_forces, &
360 : debug_stress, use_virial
361 : REAL(KIND=dp) :: fconv
362 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
363 : TYPE(cell_type), POINTER :: cell
364 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
365 : TYPE(dft_control_type), POINTER :: dft_control
366 : TYPE(mp_para_env_type), POINTER :: para_env
367 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
368 12 : POINTER :: sab_orb
369 : TYPE(virial_type), POINTER :: virial
370 :
371 12 : CALL timeset(routineN, handle)
372 :
373 12 : debug_forces = .TRUE.
374 12 : debug_stress = .TRUE.
375 12 : iounit = unit_nr
376 :
377 12 : calculate_forces = .TRUE.
378 :
379 : ! no k-points possible
380 : CALL get_qs_env(qs_env=qs_env, &
381 : cell=cell, &
382 : dft_control=dft_control, &
383 : para_env=para_env, &
384 12 : virial=virial)
385 12 : nimages = dft_control%nimages
386 12 : IF (nimages /= 1) THEN
387 0 : CPABORT("K-points not implemented")
388 : END IF
389 :
390 : ! check for virial
391 12 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
392 :
393 12 : fconv = 1.0E-9_dp*pascal/cell%deth
394 12 : IF (debug_stress .AND. use_virial) THEN
395 0 : sttot = virial%pv_virial
396 : END IF
397 :
398 : ! initialize src matrix
399 12 : NULLIFY (scrm)
400 12 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
401 12 : ALLOCATE (scrm(1, 1)%matrix)
402 12 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix)
403 12 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
404 12 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
405 :
406 : ! kinetic energy
407 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm, matrix_p=matrix_p, &
408 : calculate_forces=calculate_forces, &
409 12 : debug_forces=debug_forces, debug_stress=debug_stress)
410 :
411 12 : nder = 1
412 : CALL core_matrices(qs_env, scrm, matrix_p, calculate_forces, nder, &
413 12 : debug_forces=debug_forces, debug_stress=debug_stress)
414 :
415 12 : IF (debug_stress .AND. use_virial) THEN
416 0 : stdeb = fconv*(virial%pv_virial - sttot)
417 0 : CALL para_env%sum(stdeb)
418 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
419 0 : 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
420 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") ' '
421 : END IF
422 :
423 : ! delete scr matrix
424 12 : CALL dbcsr_deallocate_matrix_set(scrm)
425 :
426 12 : CALL timestop(handle)
427 :
428 12 : END SUBROUTINE ec_debug_force
429 :
430 : ! **************************************************************************************************
431 :
432 : ! **************************************************************************************************
433 : !> \brief ...
434 : !> \param qs_env ...
435 : !> \param ec_env ...
436 : !> \param calc_forces ...
437 : !> \param unit_nr ...
438 : ! **************************************************************************************************
439 44 : SUBROUTINE ec_ext_setup(qs_env, ec_env, calc_forces, unit_nr)
440 : TYPE(qs_environment_type), POINTER :: qs_env
441 : TYPE(energy_correction_type), POINTER :: ec_env
442 : LOGICAL, INTENT(IN) :: calc_forces
443 : INTEGER, INTENT(IN) :: unit_nr
444 :
445 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_setup'
446 :
447 : CHARACTER(LEN=default_string_length) :: headline
448 : INTEGER :: handle, ispin, nao, nocc, nspins
449 : REAL(KIND=dp) :: a_max, c_max
450 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
451 : TYPE(cp_fm_type) :: emat, ksmo, smo
452 44 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks, matrix_p, matrix_s
453 : TYPE(dft_control_type), POINTER :: dft_control
454 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
455 44 : POINTER :: sab_orb
456 : TYPE(qs_rho_type), POINTER :: rho
457 :
458 44 : CALL timeset(routineN, handle)
459 :
460 : CALL get_qs_env(qs_env=qs_env, &
461 : dft_control=dft_control, &
462 : sab_orb=sab_orb, &
463 : rho=rho, &
464 : matrix_s_kp=matrix_s, &
465 : matrix_ks_kp=matrix_ks, &
466 44 : matrix_h_kp=matrix_h)
467 :
468 44 : nspins = dft_control%nspins
469 :
470 : ! KS Hamiltonian matrix
471 44 : IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks)
472 44 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1)
473 44 : headline = "HAMILTONIAN MATRIX"
474 88 : DO ispin = 1, nspins
475 44 : ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix)
476 : CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), &
477 44 : template=matrix_ks(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
478 44 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, sab_orb)
479 88 : CALL dbcsr_copy(ec_env%matrix_ks(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix)
480 : END DO
481 :
482 : ! Overlap matrix
483 44 : IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s)
484 44 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_s, 1, 1)
485 44 : headline = "OVERLAP MATRIX"
486 44 : ALLOCATE (ec_env%matrix_s(1, 1)%matrix)
487 : CALL dbcsr_create(ec_env%matrix_s(1, 1)%matrix, name=TRIM(headline), &
488 44 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
489 44 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_s(1, 1)%matrix, sab_orb)
490 44 : CALL dbcsr_copy(ec_env%matrix_s(1, 1)%matrix, matrix_s(1, 1)%matrix)
491 :
492 : ! density matrix
493 : ! Get density matrix of reference calculation
494 44 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
495 44 : IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p)
496 44 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1)
497 44 : headline = "DENSITY MATRIX"
498 88 : DO ispin = 1, nspins
499 44 : ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix)
500 : CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), &
501 44 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
502 44 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, sab_orb)
503 88 : CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix)
504 : END DO
505 :
506 44 : IF (calc_forces) THEN
507 : ! energy weighted density matrix
508 : ! for security, we recalculate W here (stored in qs_env)
509 12 : IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w)
510 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1)
511 12 : headline = "ENERGY WEIGHTED DENSITY MATRIX"
512 24 : DO ispin = 1, nspins
513 12 : ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix)
514 : CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), &
515 12 : template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric)
516 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, sab_orb)
517 24 : CALL dbcsr_set(ec_env%matrix_w(ispin, 1)%matrix, 0.0_dp)
518 : END DO
519 :
520 : ! hz matrix
521 12 : IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz)
522 12 : CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins)
523 12 : headline = "Hz MATRIX"
524 24 : DO ispin = 1, nspins
525 12 : ALLOCATE (ec_env%matrix_hz(ispin)%matrix)
526 : CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, name=TRIM(headline), &
527 12 : template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
528 12 : CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_hz(ispin)%matrix, sab_orb)
529 24 : CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp)
530 : END DO
531 :
532 : ! Test for consistency of orbitals and KS matrix
533 24 : DO ispin = 1, nspins
534 12 : mat_struct => ec_env%mo_occ(ispin)%matrix_struct
535 12 : CALL cp_fm_create(ksmo, mat_struct)
536 12 : CALL cp_fm_get_info(ksmo, nrow_global=nao, ncol_global=nocc)
537 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%mo_occ(ispin), &
538 12 : ksmo, nocc, alpha=1.0_dp, beta=0.0_dp)
539 12 : CALL cp_fm_create(smo, mat_struct)
540 : CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_s(1, 1)%matrix, ec_env%mo_occ(ispin), &
541 12 : smo, nocc, alpha=1.0_dp, beta=0.0_dp)
542 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
543 12 : template_fmstruct=mat_struct)
544 12 : CALL cp_fm_create(emat, fm_struct)
545 12 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, ec_env%mo_occ(ispin), ksmo, 0.0_dp, emat)
546 12 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, -1.0_dp, smo, emat, 1.0_dp, ksmo)
547 12 : CALL cp_fm_maxabsval(ksmo, a_max)
548 12 : CALL cp_fm_struct_release(fm_struct)
549 12 : CALL cp_fm_release(smo)
550 12 : CALL cp_fm_release(ksmo)
551 12 : CALL cp_fm_release(emat)
552 12 : CALL cp_fm_maxabsval(ec_env%mo_occ(ispin), c_max)
553 72 : IF (unit_nr > 0) THEN
554 6 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO coeficients", ispin, c_max
555 6 : WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO gradients", ispin, a_max
556 : END IF
557 : END DO
558 : END IF
559 :
560 44 : CALL timestop(handle)
561 :
562 44 : END SUBROUTINE ec_ext_setup
563 :
564 : ! **************************************************************************************************
565 : !> \brief ...
566 : !> \param cpmos ...
567 : !> \param mo_ref ...
568 : !> \param mo_occ ...
569 : !> \param matrix_s ...
570 : !> \param unit_nr ...
571 : ! **************************************************************************************************
572 12 : SUBROUTINE align_vectors(cpmos, mo_ref, mo_occ, matrix_s, unit_nr)
573 : TYPE(cp_fm_type), INTENT(IN) :: cpmos, mo_ref, mo_occ
574 : TYPE(dbcsr_type), POINTER :: matrix_s
575 : INTEGER, INTENT(IN) :: unit_nr
576 :
577 : CHARACTER(LEN=*), PARAMETER :: routineN = 'align_vectors'
578 :
579 : INTEGER :: handle, i, nao, nocc, scg
580 : REAL(KIND=dp) :: a_max
581 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag
582 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
583 : TYPE(cp_fm_type) :: emat, smo
584 :
585 12 : CALL timeset(routineN, handle)
586 :
587 12 : mat_struct => cpmos%matrix_struct
588 12 : CALL cp_fm_create(smo, mat_struct)
589 12 : CALL cp_fm_get_info(smo, nrow_global=nao, ncol_global=nocc)
590 12 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_occ, smo, nocc, alpha=1.0_dp, beta=0.0_dp)
591 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, &
592 12 : template_fmstruct=mat_struct)
593 12 : CALL cp_fm_create(emat, fm_struct)
594 12 : CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, mo_ref, smo, 0.0_dp, emat)
595 12 : CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, cpmos, emat, 0.0_dp, smo)
596 12 : CALL cp_fm_to_fm(smo, cpmos)
597 12 : CALL cp_fm_to_fm(mo_occ, mo_ref)
598 : !
599 36 : ALLOCATE (diag(nocc))
600 12 : CALL cp_fm_get_diag(emat, diag)
601 58 : a_max = nocc - SUM(diag)
602 12 : scg = 0
603 58 : DO i = 1, nocc
604 58 : IF (ABS(diag(i) + 1.0_dp) < 0.001) scg = scg + 1
605 : END DO
606 12 : IF (unit_nr > 0) THEN
607 6 : WRITE (unit_nr, "(T3,A,T61,F20.8)") "External:: Orbital rotation index", a_max
608 6 : WRITE (unit_nr, "(T3,A,T71,I10)") "External:: Number of orbital phase changes", scg
609 : END IF
610 :
611 12 : DEALLOCATE (diag)
612 12 : CALL cp_fm_struct_release(fm_struct)
613 12 : CALL cp_fm_release(smo)
614 12 : CALL cp_fm_release(emat)
615 :
616 12 : CALL timestop(handle)
617 :
618 24 : END SUBROUTINE align_vectors
619 :
620 : ! **************************************************************************************************
621 : !> \brief ...
622 : !> \param qs_env ...
623 : !> \param matrix_w ...
624 : !> \param unit_nr ...
625 : ! **************************************************************************************************
626 0 : SUBROUTINE matrix_w_forces(qs_env, matrix_w, unit_nr)
627 : TYPE(qs_environment_type), POINTER :: qs_env
628 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w
629 : INTEGER, INTENT(IN) :: unit_nr
630 :
631 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_w_forces'
632 :
633 : INTEGER :: handle, iounit, nder, nimages
634 : LOGICAL :: debug_forces, debug_stress, use_virial
635 : REAL(KIND=dp) :: fconv
636 : REAL(KIND=dp), DIMENSION(3) :: fodeb
637 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot
638 : TYPE(cell_type), POINTER :: cell
639 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm
640 : TYPE(dft_control_type), POINTER :: dft_control
641 : TYPE(mp_para_env_type), POINTER :: para_env
642 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
643 0 : POINTER :: sab_orb
644 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
645 : TYPE(qs_ks_env_type), POINTER :: ks_env
646 : TYPE(virial_type), POINTER :: virial
647 :
648 0 : CALL timeset(routineN, handle)
649 :
650 0 : debug_forces = .TRUE.
651 0 : debug_stress = .TRUE.
652 :
653 0 : iounit = unit_nr
654 :
655 : ! no k-points possible
656 : CALL get_qs_env(qs_env=qs_env, &
657 : cell=cell, &
658 : dft_control=dft_control, &
659 : force=force, &
660 : ks_env=ks_env, &
661 : sab_orb=sab_orb, &
662 : para_env=para_env, &
663 0 : virial=virial)
664 0 : nimages = dft_control%nimages
665 0 : IF (nimages /= 1) THEN
666 0 : CPABORT("K-points not implemented")
667 : END IF
668 :
669 : ! check for virial
670 0 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
671 :
672 0 : fconv = 1.0E-9_dp*pascal/cell%deth
673 0 : IF (debug_stress .AND. use_virial) THEN
674 0 : sttot = virial%pv_virial
675 : END IF
676 :
677 : ! initialize src matrix
678 0 : NULLIFY (scrm)
679 0 : CALL dbcsr_allocate_matrix_set(scrm, 1, 1)
680 0 : ALLOCATE (scrm(1, 1)%matrix)
681 0 : CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_w(1, 1)%matrix)
682 0 : CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb)
683 :
684 0 : nder = 1
685 0 : IF (SIZE(matrix_w, 1) == 2) THEN
686 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
687 0 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
688 : END IF
689 :
690 : ! Overlap and kinetic energy matrices
691 0 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
692 0 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
693 : CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, &
694 : matrix_name="OVERLAP MATRIX", &
695 : basis_type_a="ORB", &
696 : basis_type_b="ORB", &
697 : sab_nl=sab_orb, calculate_forces=.TRUE., &
698 0 : matrixkp_p=matrix_w)
699 :
700 : IF (debug_forces) THEN
701 0 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
702 0 : CALL para_env%sum(fodeb)
703 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb
704 : END IF
705 0 : IF (debug_stress .AND. use_virial) THEN
706 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
707 0 : CALL para_env%sum(stdeb)
708 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
709 0 : 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb)
710 : END IF
711 0 : IF (SIZE(matrix_w, 1) == 2) THEN
712 : CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, &
713 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
714 : END IF
715 :
716 : ! delete scrm matrix
717 0 : CALL dbcsr_deallocate_matrix_set(scrm)
718 :
719 0 : CALL timestop(handle)
720 :
721 0 : END SUBROUTINE matrix_w_forces
722 :
723 : ! **************************************************************************************************
724 : !> \brief ...
725 : !> \param qs_env ...
726 : !> \param cpmos ...
727 : !> \param mo_occ ...
728 : !> \param matrix_r ...
729 : !> \param unit_nr ...
730 : ! **************************************************************************************************
731 12 : SUBROUTINE matrix_r_forces(qs_env, cpmos, mo_occ, matrix_r, unit_nr)
732 : TYPE(qs_environment_type), POINTER :: qs_env
733 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ
734 : TYPE(dbcsr_type), POINTER :: matrix_r
735 : INTEGER, INTENT(IN) :: unit_nr
736 :
737 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_r_forces'
738 :
739 : INTEGER :: handle, iounit, ispin, nao, nocc, nspins
740 : LOGICAL :: debug_forces, debug_stress, use_virial
741 : REAL(KIND=dp) :: fconv, focc
742 : REAL(KIND=dp), DIMENSION(3) :: fodeb
743 : REAL(KIND=dp), DIMENSION(3, 3) :: stdeb
744 : TYPE(cell_type), POINTER :: cell
745 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct
746 : TYPE(cp_fm_type) :: chcmat, rcvec
747 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: scrm
748 : TYPE(dft_control_type), POINTER :: dft_control
749 : TYPE(mp_para_env_type), POINTER :: para_env
750 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
751 12 : POINTER :: sab_orb
752 12 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
753 : TYPE(qs_ks_env_type), POINTER :: ks_env
754 : TYPE(virial_type), POINTER :: virial
755 :
756 12 : CALL timeset(routineN, handle)
757 :
758 12 : debug_forces = .TRUE.
759 12 : debug_stress = .TRUE.
760 12 : iounit = unit_nr
761 :
762 12 : nspins = SIZE(mo_occ)
763 12 : focc = 1.0_dp
764 12 : IF (nspins == 1) focc = 2.0_dp
765 12 : focc = 0.25_dp*focc
766 :
767 12 : CALL dbcsr_set(matrix_r, 0.0_dp)
768 24 : DO ispin = 1, nspins
769 12 : CALL cp_fm_get_info(cpmos(ispin), matrix_struct=fm_struct, nrow_global=nao, ncol_global=nocc)
770 12 : CALL cp_fm_create(rcvec, fm_struct)
771 12 : CALL cp_fm_struct_create(mat_struct, nrow_global=nocc, ncol_global=nocc, template_fmstruct=fm_struct)
772 12 : CALL cp_fm_create(chcmat, mat_struct)
773 12 : CALL parallel_gemm("T", "N", nocc, nocc, nao, 1.0_dp, mo_occ(ispin), cpmos(ispin), 0.0_dp, chcmat)
774 12 : CALL parallel_gemm("N", "N", nao, nocc, nocc, 1.0_dp, mo_occ(ispin), chcmat, 0.0_dp, rcvec)
775 12 : CALL cp_dbcsr_plus_fm_fm_t(matrix_r, matrix_v=rcvec, matrix_g=mo_occ(ispin), ncol=nocc, alpha=focc)
776 12 : CALL cp_fm_struct_release(mat_struct)
777 12 : CALL cp_fm_release(rcvec)
778 48 : CALL cp_fm_release(chcmat)
779 : END DO
780 :
781 : CALL get_qs_env(qs_env=qs_env, &
782 : cell=cell, &
783 : dft_control=dft_control, &
784 : force=force, &
785 : ks_env=ks_env, &
786 : sab_orb=sab_orb, &
787 : para_env=para_env, &
788 12 : virial=virial)
789 : ! check for virial
790 12 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
791 12 : fconv = 1.0E-9_dp*pascal/cell%deth
792 :
793 48 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
794 12 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
795 12 : NULLIFY (scrm)
796 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
797 : matrix_name="OVERLAP MATRIX", &
798 : basis_type_a="ORB", basis_type_b="ORB", &
799 : sab_nl=sab_orb, calculate_forces=.TRUE., &
800 12 : matrix_p=matrix_r)
801 : IF (debug_forces) THEN
802 48 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
803 12 : CALL para_env%sum(fodeb)
804 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
805 : END IF
806 12 : IF (debug_stress .AND. use_virial) THEN
807 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
808 0 : CALL para_env%sum(stdeb)
809 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
810 0 : 'STRESS| Wz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
811 : END IF
812 12 : CALL dbcsr_deallocate_matrix_set(scrm)
813 :
814 12 : CALL timestop(handle)
815 :
816 12 : END SUBROUTINE matrix_r_forces
817 :
818 : END MODULE ec_external
|
