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 : !> \par History
10 : !> Add CP2K error reporting, new add_force routine [07.2014,JGH]
11 : !> \author MK (03.06.2002)
12 : ! **************************************************************************************************
13 : MODULE qs_force_types
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE cp_log_handling, ONLY: cp_get_default_logger,&
18 : cp_logger_get_default_io_unit,&
19 : cp_logger_type
20 : USE kinds, ONLY: dp
21 : USE message_passing, ONLY: mp_para_env_type
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_force_types'
26 : PRIVATE
27 :
28 : TYPE qs_force_type
29 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: all_potential => NULL(), &
30 : cneo_potential => NULL(), &
31 : core_overlap => NULL(), &
32 : gth_ppl => NULL(), &
33 : gth_nlcc => NULL(), &
34 : gth_ppnl => NULL(), &
35 : kinetic => NULL(), &
36 : overlap => NULL(), &
37 : overlap_admm => NULL(), &
38 : rho_core => NULL(), &
39 : rho_elec => NULL(), &
40 : rho_lri_elec => NULL(), &
41 : rho_cneo_nuc => NULL(), &
42 : vhxc_atom => NULL(), &
43 : g0s_Vh_elec => NULL(), &
44 : repulsive => NULL(), &
45 : dispersion => NULL(), &
46 : gcp => NULL(), &
47 : other => NULL(), &
48 : ch_pulay => NULL(), &
49 : fock_4c => NULL(), &
50 : ehrenfest => NULL(), &
51 : efield => NULL(), &
52 : eev => NULL(), &
53 : mp2_non_sep => NULL(), &
54 : total => NULL()
55 : END TYPE qs_force_type
56 :
57 : PUBLIC :: qs_force_type
58 :
59 : PUBLIC :: allocate_qs_force, &
60 : add_qs_force, &
61 : deallocate_qs_force, &
62 : replicate_qs_force, &
63 : sum_qs_force, &
64 : get_qs_force, &
65 : put_qs_force, &
66 : total_qs_force, &
67 : zero_qs_force, &
68 : write_forces_debug
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief Allocate a Quickstep force data structure.
74 : !> \param qs_force ...
75 : !> \param natom_of_kind ...
76 : !> \date 05.06.2002
77 : !> \author MK
78 : !> \version 1.0
79 : ! **************************************************************************************************
80 4351 : SUBROUTINE allocate_qs_force(qs_force, natom_of_kind)
81 :
82 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
83 : INTEGER, DIMENSION(:), INTENT(IN) :: natom_of_kind
84 :
85 : INTEGER :: ikind, n, nkind
86 :
87 4351 : IF (ASSOCIATED(qs_force)) CALL deallocate_qs_force(qs_force)
88 :
89 4351 : nkind = SIZE(natom_of_kind)
90 :
91 21505 : ALLOCATE (qs_force(nkind))
92 :
93 12803 : DO ikind = 1, nkind
94 8452 : n = natom_of_kind(ikind)
95 25356 : ALLOCATE (qs_force(ikind)%all_potential(3, n))
96 16904 : ALLOCATE (qs_force(ikind)%cneo_potential(3, n))
97 16904 : ALLOCATE (qs_force(ikind)%core_overlap(3, n))
98 16904 : ALLOCATE (qs_force(ikind)%gth_ppl(3, n))
99 16904 : ALLOCATE (qs_force(ikind)%gth_nlcc(3, n))
100 16904 : ALLOCATE (qs_force(ikind)%gth_ppnl(3, n))
101 16904 : ALLOCATE (qs_force(ikind)%kinetic(3, n))
102 16904 : ALLOCATE (qs_force(ikind)%overlap(3, n))
103 16904 : ALLOCATE (qs_force(ikind)%overlap_admm(3, n))
104 16904 : ALLOCATE (qs_force(ikind)%rho_core(3, n))
105 16904 : ALLOCATE (qs_force(ikind)%rho_elec(3, n))
106 16904 : ALLOCATE (qs_force(ikind)%rho_lri_elec(3, n))
107 16904 : ALLOCATE (qs_force(ikind)%rho_cneo_nuc(3, n))
108 16904 : ALLOCATE (qs_force(ikind)%vhxc_atom(3, n))
109 16904 : ALLOCATE (qs_force(ikind)%g0s_Vh_elec(3, n))
110 16904 : ALLOCATE (qs_force(ikind)%repulsive(3, n))
111 16904 : ALLOCATE (qs_force(ikind)%dispersion(3, n))
112 16904 : ALLOCATE (qs_force(ikind)%gcp(3, n))
113 16904 : ALLOCATE (qs_force(ikind)%other(3, n))
114 16904 : ALLOCATE (qs_force(ikind)%ch_pulay(3, n))
115 16904 : ALLOCATE (qs_force(ikind)%ehrenfest(3, n))
116 16904 : ALLOCATE (qs_force(ikind)%efield(3, n))
117 16904 : ALLOCATE (qs_force(ikind)%eev(3, n))
118 : ! Always initialize ch_pulay to zero..
119 101432 : qs_force(ikind)%ch_pulay = 0.0_dp
120 16904 : ALLOCATE (qs_force(ikind)%fock_4c(3, n))
121 16904 : ALLOCATE (qs_force(ikind)%mp2_non_sep(3, n))
122 21255 : ALLOCATE (qs_force(ikind)%total(3, n))
123 : END DO
124 :
125 4351 : END SUBROUTINE allocate_qs_force
126 :
127 : ! **************************************************************************************************
128 : !> \brief Deallocate a Quickstep force data structure.
129 : !> \param qs_force ...
130 : !> \date 05.06.2002
131 : !> \author MK
132 : !> \version 1.0
133 : ! **************************************************************************************************
134 4351 : SUBROUTINE deallocate_qs_force(qs_force)
135 :
136 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
137 :
138 : INTEGER :: ikind, nkind
139 :
140 4351 : CPASSERT(ASSOCIATED(qs_force))
141 :
142 4351 : nkind = SIZE(qs_force)
143 :
144 12803 : DO ikind = 1, nkind
145 :
146 8452 : IF (ASSOCIATED(qs_force(ikind)%all_potential)) THEN
147 8452 : DEALLOCATE (qs_force(ikind)%all_potential)
148 : END IF
149 :
150 8452 : IF (ASSOCIATED(qs_force(ikind)%cneo_potential)) THEN
151 8452 : DEALLOCATE (qs_force(ikind)%cneo_potential)
152 : END IF
153 :
154 8452 : IF (ASSOCIATED(qs_force(ikind)%core_overlap)) THEN
155 8452 : DEALLOCATE (qs_force(ikind)%core_overlap)
156 : END IF
157 :
158 8452 : IF (ASSOCIATED(qs_force(ikind)%gth_ppl)) THEN
159 8452 : DEALLOCATE (qs_force(ikind)%gth_ppl)
160 : END IF
161 :
162 8452 : IF (ASSOCIATED(qs_force(ikind)%gth_nlcc)) THEN
163 8452 : DEALLOCATE (qs_force(ikind)%gth_nlcc)
164 : END IF
165 :
166 8452 : IF (ASSOCIATED(qs_force(ikind)%gth_ppnl)) THEN
167 8452 : DEALLOCATE (qs_force(ikind)%gth_ppnl)
168 : END IF
169 :
170 8452 : IF (ASSOCIATED(qs_force(ikind)%kinetic)) THEN
171 8452 : DEALLOCATE (qs_force(ikind)%kinetic)
172 : END IF
173 :
174 8452 : IF (ASSOCIATED(qs_force(ikind)%overlap)) THEN
175 8452 : DEALLOCATE (qs_force(ikind)%overlap)
176 : END IF
177 :
178 8452 : IF (ASSOCIATED(qs_force(ikind)%overlap_admm)) THEN
179 8452 : DEALLOCATE (qs_force(ikind)%overlap_admm)
180 : END IF
181 :
182 8452 : IF (ASSOCIATED(qs_force(ikind)%rho_core)) THEN
183 8452 : DEALLOCATE (qs_force(ikind)%rho_core)
184 : END IF
185 :
186 8452 : IF (ASSOCIATED(qs_force(ikind)%rho_elec)) THEN
187 8452 : DEALLOCATE (qs_force(ikind)%rho_elec)
188 : END IF
189 8452 : IF (ASSOCIATED(qs_force(ikind)%rho_lri_elec)) THEN
190 8452 : DEALLOCATE (qs_force(ikind)%rho_lri_elec)
191 : END IF
192 :
193 8452 : IF (ASSOCIATED(qs_force(ikind)%rho_cneo_nuc)) THEN
194 8452 : DEALLOCATE (qs_force(ikind)%rho_cneo_nuc)
195 : END IF
196 :
197 8452 : IF (ASSOCIATED(qs_force(ikind)%vhxc_atom)) THEN
198 8452 : DEALLOCATE (qs_force(ikind)%vhxc_atom)
199 : END IF
200 :
201 8452 : IF (ASSOCIATED(qs_force(ikind)%g0s_Vh_elec)) THEN
202 8452 : DEALLOCATE (qs_force(ikind)%g0s_Vh_elec)
203 : END IF
204 :
205 8452 : IF (ASSOCIATED(qs_force(ikind)%repulsive)) THEN
206 8452 : DEALLOCATE (qs_force(ikind)%repulsive)
207 : END IF
208 :
209 8452 : IF (ASSOCIATED(qs_force(ikind)%dispersion)) THEN
210 8452 : DEALLOCATE (qs_force(ikind)%dispersion)
211 : END IF
212 :
213 8452 : IF (ASSOCIATED(qs_force(ikind)%gcp)) THEN
214 8452 : DEALLOCATE (qs_force(ikind)%gcp)
215 : END IF
216 :
217 8452 : IF (ASSOCIATED(qs_force(ikind)%other)) THEN
218 8452 : DEALLOCATE (qs_force(ikind)%other)
219 : END IF
220 :
221 8452 : IF (ASSOCIATED(qs_force(ikind)%total)) THEN
222 8452 : DEALLOCATE (qs_force(ikind)%total)
223 : END IF
224 :
225 8452 : IF (ASSOCIATED(qs_force(ikind)%ch_pulay)) THEN
226 8452 : DEALLOCATE (qs_force(ikind)%ch_pulay)
227 : END IF
228 :
229 8452 : IF (ASSOCIATED(qs_force(ikind)%fock_4c)) THEN
230 8452 : DEALLOCATE (qs_force(ikind)%fock_4c)
231 : END IF
232 :
233 8452 : IF (ASSOCIATED(qs_force(ikind)%mp2_non_sep)) THEN
234 8452 : DEALLOCATE (qs_force(ikind)%mp2_non_sep)
235 : END IF
236 :
237 8452 : IF (ASSOCIATED(qs_force(ikind)%ehrenfest)) THEN
238 8452 : DEALLOCATE (qs_force(ikind)%ehrenfest)
239 : END IF
240 :
241 8452 : IF (ASSOCIATED(qs_force(ikind)%efield)) THEN
242 8452 : DEALLOCATE (qs_force(ikind)%efield)
243 : END IF
244 :
245 12803 : IF (ASSOCIATED(qs_force(ikind)%eev)) THEN
246 8452 : DEALLOCATE (qs_force(ikind)%eev)
247 : END IF
248 : END DO
249 :
250 4351 : DEALLOCATE (qs_force)
251 :
252 4351 : END SUBROUTINE deallocate_qs_force
253 :
254 : ! **************************************************************************************************
255 : !> \brief Initialize a Quickstep force data structure.
256 : !> \param qs_force ...
257 : !> \date 15.07.2002
258 : !> \author MK
259 : !> \version 1.0
260 : ! **************************************************************************************************
261 12781 : SUBROUTINE zero_qs_force(qs_force)
262 :
263 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
264 :
265 : INTEGER :: ikind
266 :
267 12781 : CPASSERT(ASSOCIATED(qs_force))
268 :
269 37477 : DO ikind = 1, SIZE(qs_force)
270 319980 : qs_force(ikind)%all_potential(:, :) = 0.0_dp
271 319980 : qs_force(ikind)%cneo_potential(:, :) = 0.0_dp
272 319980 : qs_force(ikind)%core_overlap(:, :) = 0.0_dp
273 319980 : qs_force(ikind)%gth_ppl(:, :) = 0.0_dp
274 319980 : qs_force(ikind)%gth_nlcc(:, :) = 0.0_dp
275 319980 : qs_force(ikind)%gth_ppnl(:, :) = 0.0_dp
276 319980 : qs_force(ikind)%kinetic(:, :) = 0.0_dp
277 319980 : qs_force(ikind)%overlap(:, :) = 0.0_dp
278 319980 : qs_force(ikind)%overlap_admm(:, :) = 0.0_dp
279 319980 : qs_force(ikind)%rho_core(:, :) = 0.0_dp
280 319980 : qs_force(ikind)%rho_elec(:, :) = 0.0_dp
281 319980 : qs_force(ikind)%rho_lri_elec(:, :) = 0.0_dp
282 319980 : qs_force(ikind)%rho_cneo_nuc(:, :) = 0.0_dp
283 319980 : qs_force(ikind)%vhxc_atom(:, :) = 0.0_dp
284 319980 : qs_force(ikind)%g0s_Vh_elec(:, :) = 0.0_dp
285 319980 : qs_force(ikind)%repulsive(:, :) = 0.0_dp
286 319980 : qs_force(ikind)%dispersion(:, :) = 0.0_dp
287 319980 : qs_force(ikind)%gcp(:, :) = 0.0_dp
288 319980 : qs_force(ikind)%other(:, :) = 0.0_dp
289 319980 : qs_force(ikind)%fock_4c(:, :) = 0.0_dp
290 319980 : qs_force(ikind)%ehrenfest(:, :) = 0.0_dp
291 319980 : qs_force(ikind)%efield(:, :) = 0.0_dp
292 319980 : qs_force(ikind)%eev(:, :) = 0.0_dp
293 319980 : qs_force(ikind)%mp2_non_sep(:, :) = 0.0_dp
294 332761 : qs_force(ikind)%total(:, :) = 0.0_dp
295 : END DO
296 :
297 12781 : END SUBROUTINE zero_qs_force
298 :
299 : ! **************************************************************************************************
300 : !> \brief Sum up two qs_force entities qs_force_out = qs_force_out + qs_force_in
301 : !> \param qs_force_out ...
302 : !> \param qs_force_in ...
303 : !> \author JGH
304 : ! **************************************************************************************************
305 1432 : SUBROUTINE sum_qs_force(qs_force_out, qs_force_in)
306 :
307 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force_out, qs_force_in
308 :
309 : INTEGER :: ikind
310 :
311 1432 : CPASSERT(ASSOCIATED(qs_force_out))
312 1432 : CPASSERT(ASSOCIATED(qs_force_in))
313 :
314 4348 : DO ikind = 1, SIZE(qs_force_out)
315 : qs_force_out(ikind)%all_potential(:, :) = qs_force_out(ikind)%all_potential(:, :) + &
316 43800 : qs_force_in(ikind)%all_potential(:, :)
317 : qs_force_out(ikind)%cneo_potential(:, :) = qs_force_out(ikind)%cneo_potential(:, :) + &
318 43800 : qs_force_in(ikind)%cneo_potential(:, :)
319 : qs_force_out(ikind)%core_overlap(:, :) = qs_force_out(ikind)%core_overlap(:, :) + &
320 43800 : qs_force_in(ikind)%core_overlap(:, :)
321 : qs_force_out(ikind)%gth_ppl(:, :) = qs_force_out(ikind)%gth_ppl(:, :) + &
322 43800 : qs_force_in(ikind)%gth_ppl(:, :)
323 : qs_force_out(ikind)%gth_nlcc(:, :) = qs_force_out(ikind)%gth_nlcc(:, :) + &
324 43800 : qs_force_in(ikind)%gth_nlcc(:, :)
325 : qs_force_out(ikind)%gth_ppnl(:, :) = qs_force_out(ikind)%gth_ppnl(:, :) + &
326 43800 : qs_force_in(ikind)%gth_ppnl(:, :)
327 : qs_force_out(ikind)%kinetic(:, :) = qs_force_out(ikind)%kinetic(:, :) + &
328 43800 : qs_force_in(ikind)%kinetic(:, :)
329 : qs_force_out(ikind)%overlap(:, :) = qs_force_out(ikind)%overlap(:, :) + &
330 43800 : qs_force_in(ikind)%overlap(:, :)
331 : qs_force_out(ikind)%overlap_admm(:, :) = qs_force_out(ikind)%overlap_admm(:, :) + &
332 43800 : qs_force_in(ikind)%overlap_admm(:, :)
333 : qs_force_out(ikind)%rho_core(:, :) = qs_force_out(ikind)%rho_core(:, :) + &
334 43800 : qs_force_in(ikind)%rho_core(:, :)
335 : qs_force_out(ikind)%rho_elec(:, :) = qs_force_out(ikind)%rho_elec(:, :) + &
336 43800 : qs_force_in(ikind)%rho_elec(:, :)
337 : qs_force_out(ikind)%rho_lri_elec(:, :) = qs_force_out(ikind)%rho_lri_elec(:, :) + &
338 43800 : qs_force_in(ikind)%rho_lri_elec(:, :)
339 : qs_force_out(ikind)%rho_cneo_nuc(:, :) = qs_force_out(ikind)%rho_cneo_nuc(:, :) + &
340 43800 : qs_force_in(ikind)%rho_cneo_nuc(:, :)
341 : qs_force_out(ikind)%vhxc_atom(:, :) = qs_force_out(ikind)%vhxc_atom(:, :) + &
342 43800 : qs_force_in(ikind)%vhxc_atom(:, :)
343 : qs_force_out(ikind)%g0s_Vh_elec(:, :) = qs_force_out(ikind)%g0s_Vh_elec(:, :) + &
344 43800 : qs_force_in(ikind)%g0s_Vh_elec(:, :)
345 : qs_force_out(ikind)%repulsive(:, :) = qs_force_out(ikind)%repulsive(:, :) + &
346 43800 : qs_force_in(ikind)%repulsive(:, :)
347 : qs_force_out(ikind)%dispersion(:, :) = qs_force_out(ikind)%dispersion(:, :) + &
348 43800 : qs_force_in(ikind)%dispersion(:, :)
349 : qs_force_out(ikind)%gcp(:, :) = qs_force_out(ikind)%gcp(:, :) + &
350 43800 : qs_force_in(ikind)%gcp(:, :)
351 : qs_force_out(ikind)%other(:, :) = qs_force_out(ikind)%other(:, :) + &
352 43800 : qs_force_in(ikind)%other(:, :)
353 : qs_force_out(ikind)%fock_4c(:, :) = qs_force_out(ikind)%fock_4c(:, :) + &
354 43800 : qs_force_in(ikind)%fock_4c(:, :)
355 : qs_force_out(ikind)%ehrenfest(:, :) = qs_force_out(ikind)%ehrenfest(:, :) + &
356 43800 : qs_force_in(ikind)%ehrenfest(:, :)
357 : qs_force_out(ikind)%efield(:, :) = qs_force_out(ikind)%efield(:, :) + &
358 43800 : qs_force_in(ikind)%efield(:, :)
359 : qs_force_out(ikind)%eev(:, :) = qs_force_out(ikind)%eev(:, :) + &
360 43800 : qs_force_in(ikind)%eev(:, :)
361 : qs_force_out(ikind)%mp2_non_sep(:, :) = qs_force_out(ikind)%mp2_non_sep(:, :) + &
362 43800 : qs_force_in(ikind)%mp2_non_sep(:, :)
363 : qs_force_out(ikind)%total(:, :) = qs_force_out(ikind)%total(:, :) + &
364 45232 : qs_force_in(ikind)%total(:, :)
365 : END DO
366 :
367 1432 : END SUBROUTINE sum_qs_force
368 :
369 : ! **************************************************************************************************
370 : !> \brief Replicate and sum up the force
371 : !> \param qs_force ...
372 : !> \param para_env ...
373 : !> \date 25.05.2016
374 : !> \author JHU
375 : !> \version 1.0
376 : ! **************************************************************************************************
377 10625 : SUBROUTINE replicate_qs_force(qs_force, para_env)
378 :
379 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
380 : TYPE(mp_para_env_type), POINTER :: para_env
381 :
382 : INTEGER :: ikind
383 :
384 : ! *** replicate forces ***
385 31381 : DO ikind = 1, SIZE(qs_force)
386 558780 : CALL para_env%sum(qs_force(ikind)%overlap)
387 558780 : CALL para_env%sum(qs_force(ikind)%overlap_admm)
388 558780 : CALL para_env%sum(qs_force(ikind)%kinetic)
389 558780 : CALL para_env%sum(qs_force(ikind)%gth_ppl)
390 558780 : CALL para_env%sum(qs_force(ikind)%gth_nlcc)
391 558780 : CALL para_env%sum(qs_force(ikind)%gth_ppnl)
392 558780 : CALL para_env%sum(qs_force(ikind)%all_potential)
393 558780 : CALL para_env%sum(qs_force(ikind)%cneo_potential)
394 558780 : CALL para_env%sum(qs_force(ikind)%core_overlap)
395 558780 : CALL para_env%sum(qs_force(ikind)%rho_core)
396 558780 : CALL para_env%sum(qs_force(ikind)%rho_elec)
397 558780 : CALL para_env%sum(qs_force(ikind)%rho_lri_elec)
398 558780 : CALL para_env%sum(qs_force(ikind)%rho_cneo_nuc)
399 558780 : CALL para_env%sum(qs_force(ikind)%vhxc_atom)
400 558780 : CALL para_env%sum(qs_force(ikind)%g0s_Vh_elec)
401 558780 : CALL para_env%sum(qs_force(ikind)%fock_4c)
402 558780 : CALL para_env%sum(qs_force(ikind)%mp2_non_sep)
403 558780 : CALL para_env%sum(qs_force(ikind)%repulsive)
404 558780 : CALL para_env%sum(qs_force(ikind)%dispersion)
405 558780 : CALL para_env%sum(qs_force(ikind)%gcp)
406 558780 : CALL para_env%sum(qs_force(ikind)%ehrenfest)
407 :
408 : qs_force(ikind)%total(:, :) = qs_force(ikind)%total(:, :) + &
409 : qs_force(ikind)%core_overlap(:, :) + &
410 : qs_force(ikind)%gth_ppl(:, :) + &
411 : qs_force(ikind)%gth_nlcc(:, :) + &
412 : qs_force(ikind)%gth_ppnl(:, :) + &
413 : qs_force(ikind)%all_potential(:, :) + &
414 : qs_force(ikind)%cneo_potential(:, :) + &
415 : qs_force(ikind)%kinetic(:, :) + &
416 : qs_force(ikind)%overlap(:, :) + &
417 : qs_force(ikind)%overlap_admm(:, :) + &
418 : qs_force(ikind)%rho_core(:, :) + &
419 : qs_force(ikind)%rho_elec(:, :) + &
420 : qs_force(ikind)%rho_lri_elec(:, :) + &
421 : qs_force(ikind)%rho_cneo_nuc(:, :) + &
422 : qs_force(ikind)%vhxc_atom(:, :) + &
423 : qs_force(ikind)%g0s_Vh_elec(:, :) + &
424 : qs_force(ikind)%fock_4c(:, :) + &
425 : qs_force(ikind)%mp2_non_sep(:, :) + &
426 : qs_force(ikind)%repulsive(:, :) + &
427 : qs_force(ikind)%dispersion(:, :) + &
428 : qs_force(ikind)%gcp(:, :) + &
429 : qs_force(ikind)%ehrenfest(:, :) + &
430 : qs_force(ikind)%efield(:, :) + &
431 300393 : qs_force(ikind)%eev(:, :)
432 : END DO
433 :
434 10625 : END SUBROUTINE replicate_qs_force
435 :
436 : ! **************************************************************************************************
437 : !> \brief Add force to a force_type variable.
438 : !> \param force Input force, dimension (3,natom)
439 : !> \param qs_force The force type variable to be used
440 : !> \param forcetype ...
441 : !> \param atomic_kind_set ...
442 : !> \par History
443 : !> 07.2014 JGH
444 : !> \author JGH
445 : ! **************************************************************************************************
446 1112 : SUBROUTINE add_qs_force(force, qs_force, forcetype, atomic_kind_set)
447 :
448 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: force
449 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
450 : CHARACTER(LEN=*), INTENT(IN) :: forcetype
451 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
452 :
453 : INTEGER :: ia, iatom, ikind, natom_kind
454 : TYPE(atomic_kind_type), POINTER :: atomic_kind
455 :
456 : ! ------------------------------------------------------------------------
457 :
458 1112 : CPASSERT(ASSOCIATED(qs_force))
459 :
460 1112 : SELECT CASE (forcetype)
461 : CASE ("overlap_admm")
462 3102 : DO ikind = 1, SIZE(atomic_kind_set, 1)
463 1990 : atomic_kind => atomic_kind_set(ikind)
464 1990 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
465 6282 : DO ia = 1, natom_kind
466 3180 : iatom = atomic_kind%atom_list(ia)
467 14710 : qs_force(ikind)%overlap_admm(:, ia) = qs_force(ikind)%overlap_admm(:, ia) + force(:, iatom)
468 : END DO
469 : END DO
470 : CASE DEFAULT
471 1112 : CPABORT("")
472 : END SELECT
473 :
474 1112 : END SUBROUTINE add_qs_force
475 :
476 : ! **************************************************************************************************
477 : !> \brief Put force to a force_type variable.
478 : !> \param force Input force, dimension (3,natom)
479 : !> \param qs_force The force type variable to be used
480 : !> \param forcetype ...
481 : !> \param atomic_kind_set ...
482 : !> \par History
483 : !> 09.2019 JGH
484 : !> \author JGH
485 : ! **************************************************************************************************
486 0 : SUBROUTINE put_qs_force(force, qs_force, forcetype, atomic_kind_set)
487 :
488 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: force
489 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
490 : CHARACTER(LEN=*), INTENT(IN) :: forcetype
491 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
492 :
493 : INTEGER :: ia, iatom, ikind, natom_kind
494 : TYPE(atomic_kind_type), POINTER :: atomic_kind
495 :
496 : ! ------------------------------------------------------------------------
497 :
498 0 : SELECT CASE (forcetype)
499 : CASE ("dispersion")
500 0 : DO ikind = 1, SIZE(atomic_kind_set, 1)
501 0 : atomic_kind => atomic_kind_set(ikind)
502 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
503 0 : DO ia = 1, natom_kind
504 0 : iatom = atomic_kind%atom_list(ia)
505 0 : qs_force(ikind)%dispersion(:, ia) = force(:, iatom)
506 : END DO
507 : END DO
508 : CASE DEFAULT
509 0 : CPABORT("")
510 : END SELECT
511 :
512 0 : END SUBROUTINE put_qs_force
513 :
514 : ! **************************************************************************************************
515 : !> \brief Get force from a force_type variable.
516 : !> \param force Input force, dimension (3,natom)
517 : !> \param qs_force The force type variable to be used
518 : !> \param forcetype ...
519 : !> \param atomic_kind_set ...
520 : !> \par History
521 : !> 09.2019 JGH
522 : !> \author JGH
523 : ! **************************************************************************************************
524 0 : SUBROUTINE get_qs_force(force, qs_force, forcetype, atomic_kind_set)
525 :
526 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force
527 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
528 : CHARACTER(LEN=*), INTENT(IN) :: forcetype
529 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
530 :
531 : INTEGER :: ia, iatom, ikind, natom_kind
532 : TYPE(atomic_kind_type), POINTER :: atomic_kind
533 :
534 : ! ------------------------------------------------------------------------
535 :
536 0 : SELECT CASE (forcetype)
537 : CASE ("dispersion")
538 0 : DO ikind = 1, SIZE(atomic_kind_set, 1)
539 0 : atomic_kind => atomic_kind_set(ikind)
540 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
541 0 : DO ia = 1, natom_kind
542 0 : iatom = atomic_kind%atom_list(ia)
543 0 : force(:, iatom) = qs_force(ikind)%dispersion(:, ia)
544 : END DO
545 : END DO
546 : CASE DEFAULT
547 0 : CPABORT("")
548 : END SELECT
549 :
550 0 : END SUBROUTINE get_qs_force
551 :
552 : ! **************************************************************************************************
553 : !> \brief Get current total force
554 : !> \param force Input force, dimension (3,natom)
555 : !> \param qs_force The force type variable to be used
556 : !> \param atomic_kind_set ...
557 : !> \par History
558 : !> 09.2019 JGH
559 : !> \author JGH
560 : ! **************************************************************************************************
561 556 : SUBROUTINE total_qs_force(force, qs_force, atomic_kind_set)
562 :
563 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force
564 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
565 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
566 :
567 : INTEGER :: ia, iatom, ikind, natom_kind
568 : TYPE(atomic_kind_type), POINTER :: atomic_kind
569 :
570 : ! ------------------------------------------------------------------------
571 :
572 7180 : force(:, :) = 0.0_dp
573 1656 : DO ikind = 1, SIZE(atomic_kind_set, 1)
574 1100 : atomic_kind => atomic_kind_set(ikind)
575 1100 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
576 3312 : DO ia = 1, natom_kind
577 1656 : iatom = atomic_kind%atom_list(ia)
578 : force(:, iatom) = qs_force(ikind)%core_overlap(:, ia) + &
579 : qs_force(ikind)%gth_ppl(:, ia) + &
580 : qs_force(ikind)%gth_nlcc(:, ia) + &
581 : qs_force(ikind)%gth_ppnl(:, ia) + &
582 : qs_force(ikind)%all_potential(:, ia) + &
583 : qs_force(ikind)%cneo_potential(:, ia) + &
584 : qs_force(ikind)%kinetic(:, ia) + &
585 : qs_force(ikind)%overlap(:, ia) + &
586 : qs_force(ikind)%overlap_admm(:, ia) + &
587 : qs_force(ikind)%rho_core(:, ia) + &
588 : qs_force(ikind)%rho_elec(:, ia) + &
589 : qs_force(ikind)%rho_lri_elec(:, ia) + &
590 : qs_force(ikind)%rho_cneo_nuc(:, ia) + &
591 : qs_force(ikind)%vhxc_atom(:, ia) + &
592 : qs_force(ikind)%g0s_Vh_elec(:, ia) + &
593 : qs_force(ikind)%fock_4c(:, ia) + &
594 : qs_force(ikind)%mp2_non_sep(:, ia) + &
595 : qs_force(ikind)%repulsive(:, ia) + &
596 : qs_force(ikind)%dispersion(:, ia) + &
597 : qs_force(ikind)%gcp(:, ia) + &
598 : qs_force(ikind)%ehrenfest(:, ia) + &
599 : qs_force(ikind)%efield(:, ia) + &
600 7724 : qs_force(ikind)%eev(:, ia)
601 : END DO
602 : END DO
603 :
604 556 : END SUBROUTINE total_qs_force
605 :
606 : ! **************************************************************************************************
607 : !> \brief Write a Quickstep force data for 1 atom
608 : !> \param qs_force ...
609 : !> \param ikind ...
610 : !> \param iatom ...
611 : !> \param iunit ...
612 : !> \date 05.06.2002
613 : !> \author MK/JGH
614 : !> \version 1.0
615 : ! **************************************************************************************************
616 0 : SUBROUTINE write_forces_debug(qs_force, ikind, iatom, iunit)
617 :
618 : TYPE(qs_force_type), DIMENSION(:), POINTER :: qs_force
619 : INTEGER, INTENT(IN), OPTIONAL :: ikind, iatom, iunit
620 :
621 : CHARACTER(LEN=35) :: fmtstr2
622 : CHARACTER(LEN=48) :: fmtstr1
623 : INTEGER :: iounit, jatom, jkind
624 : REAL(KIND=dp), DIMENSION(3) :: total
625 : TYPE(cp_logger_type), POINTER :: logger
626 :
627 0 : IF (PRESENT(iunit)) THEN
628 0 : iounit = iunit
629 : ELSE
630 0 : NULLIFY (logger)
631 0 : logger => cp_get_default_logger()
632 0 : iounit = cp_logger_get_default_io_unit(logger)
633 : END IF
634 0 : IF (PRESENT(ikind)) THEN
635 0 : jkind = ikind
636 : ELSE
637 0 : jkind = 1
638 : END IF
639 0 : IF (PRESENT(iatom)) THEN
640 0 : jatom = iatom
641 : ELSE
642 0 : jatom = 1
643 : END IF
644 :
645 0 : IF (iounit > 0) THEN
646 :
647 0 : fmtstr1 = "(/,T2,A,/,T3,A,T11,A,T23,A,T40,A1,2(17X,A1))"
648 0 : fmtstr2 = "((T2,I5,4X,I4,T18,A,T34,3F18.12))"
649 :
650 : WRITE (UNIT=iounit, FMT=fmtstr1) &
651 0 : "FORCES [a.u.]", "Atom", "Kind", "Component", "X", "Y", "Z"
652 :
653 : total(1:3) = qs_force(jkind)%overlap(1:3, jatom) &
654 : + qs_force(jkind)%overlap_admm(1:3, jatom) &
655 : + qs_force(jkind)%kinetic(1:3, jatom) &
656 : + qs_force(jkind)%gth_ppl(1:3, jatom) &
657 : + qs_force(jkind)%gth_ppnl(1:3, jatom) &
658 : + qs_force(jkind)%gth_nlcc(1:3, jatom) &
659 : + qs_force(jkind)%all_potential(1:3, jatom) &
660 : + qs_force(jkind)%cneo_potential(1:3, jatom) &
661 : + qs_force(jkind)%rho_cneo_nuc(1:3, jatom) &
662 : + qs_force(jkind)%core_overlap(1:3, jatom) &
663 : + qs_force(jkind)%rho_core(1:3, jatom) &
664 : + qs_force(jkind)%rho_elec(1:3, jatom) &
665 : + qs_force(jkind)%rho_lri_elec(1:3, jatom) &
666 : + qs_force(jkind)%vhxc_atom(1:3, jatom) &
667 : + qs_force(jkind)%g0s_Vh_elec(1:3, jatom) &
668 : + qs_force(jkind)%dispersion(1:3, jatom) &
669 : + qs_force(jkind)%repulsive(1:3, jatom) &
670 : + qs_force(jkind)%gcp(1:3, jatom) &
671 : + qs_force(jkind)%efield(1:3, jatom) &
672 : + qs_force(jkind)%eev(1:3, jatom) &
673 : + qs_force(jkind)%ehrenfest(1:3, jatom) &
674 : + qs_force(jkind)%fock_4c(1:3, jatom) &
675 0 : + qs_force(jkind)%mp2_non_sep(1:3, jatom)
676 :
677 : WRITE (UNIT=iounit, FMT=fmtstr2) &
678 0 : jatom, jkind, " overlap", qs_force(jkind)%overlap(1:3, jatom), &
679 0 : jatom, jkind, " overlap_admm", qs_force(jkind)%overlap_admm(1:3, jatom), &
680 0 : jatom, jkind, " kinetic", qs_force(jkind)%kinetic(1:3, jatom), &
681 0 : jatom, jkind, " gth_ppl", qs_force(jkind)%gth_ppl(1:3, jatom), &
682 0 : jatom, jkind, " gth_ppnl", qs_force(jkind)%gth_ppnl(1:3, jatom), &
683 0 : jatom, jkind, " gth_nlcc", qs_force(jkind)%gth_nlcc(1:3, jatom), &
684 0 : jatom, jkind, " all_potential", qs_force(jkind)%all_potential(1:3, jatom), &
685 0 : jatom, jkind, "cneo_potential", qs_force(jkind)%cneo_potential(1:3, jatom), &
686 0 : jatom, jkind, " rho_cneo_nuc", qs_force(jkind)%rho_cneo_nuc(1:3, jatom), &
687 0 : jatom, jkind, " core_overlap", qs_force(jkind)%core_overlap(1:3, jatom), &
688 0 : jatom, jkind, " rho_core", qs_force(jkind)%rho_core(1:3, jatom), &
689 0 : jatom, jkind, " rho_elec", qs_force(jkind)%rho_elec(1:3, jatom), &
690 0 : jatom, jkind, " rho_lri_elec", qs_force(jkind)%rho_lri_elec(1:3, jatom), &
691 0 : jatom, jkind, " vhxc_atom", qs_force(jkind)%vhxc_atom(1:3, jatom), &
692 0 : jatom, jkind, " g0s_Vh_elec", qs_force(jkind)%g0s_Vh_elec(1:3, jatom), &
693 0 : jatom, jkind, " dispersion", qs_force(jkind)%dispersion(1:3, jatom), &
694 0 : jatom, jkind, " repulsive", qs_force(jkind)%repulsive(1:3, jatom), &
695 0 : jatom, jkind, " gcp", qs_force(jkind)%gcp(1:3, jatom), &
696 0 : jatom, jkind, " efield", qs_force(jkind)%efield(1:3, jatom), &
697 0 : jatom, jkind, " eev", qs_force(jkind)%eev(1:3, jatom), &
698 0 : jatom, jkind, " ehrenfest", qs_force(jkind)%ehrenfest(1:3, jatom), &
699 0 : jatom, jkind, " fock_4c", qs_force(jkind)%fock_4c(1:3, jatom), &
700 0 : jatom, jkind, " mp2_non_sep", qs_force(jkind)%mp2_non_sep(1:3, jatom), &
701 0 : jatom, jkind, " total", total(1:3)
702 :
703 : END IF
704 :
705 0 : END SUBROUTINE write_forces_debug
706 :
707 0 : END MODULE qs_force_types
|
