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