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