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