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