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 : !> Torsions added(DG)05-Dec-2000
11 : !> Variable names changed(DG)05-Dec-2000
12 : !> CJM SEPT-12-2002: int_env is now passed
13 : !> CJM NOV-30-2003: only uses fist_env
14 : !> \author CJM & JGH
15 : ! **************************************************************************************************
16 : MODULE fist_force
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE atprop_types, ONLY: atprop_type
20 : USE cell_methods, ONLY: init_cell
21 : USE cell_types, ONLY: cell_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_type
24 : USE cp_output_handling, ONLY: cp_p_file,&
25 : cp_print_key_finished_output,&
26 : cp_print_key_should_output,&
27 : cp_print_key_unit_nr
28 : USE cp_subsys_types, ONLY: cp_subsys_get,&
29 : cp_subsys_type
30 : USE cp_units, ONLY: cp_unit_from_cp2k
31 : USE distribution_1d_types, ONLY: distribution_1d_type
32 : USE ewald_environment_types, ONLY: ewald_env_get,&
33 : ewald_environment_type
34 : USE ewald_pw_methods, ONLY: ewald_pw_grid_update
35 : USE ewald_pw_types, ONLY: ewald_pw_type
36 : USE ewalds, ONLY: ewald_evaluate,&
37 : ewald_print,&
38 : ewald_self,&
39 : ewald_self_atom
40 : USE ewalds_multipole, ONLY: ewald_multipole_evaluate
41 : USE exclusion_types, ONLY: exclusion_type
42 : USE fist_efield_methods, ONLY: fist_dipole,&
43 : fist_efield_energy_force
44 : USE fist_efield_types, ONLY: fist_efield_type
45 : USE fist_energy_types, ONLY: fist_energy_type
46 : USE fist_environment_types, ONLY: fist_env_get,&
47 : fist_environment_type
48 : USE fist_intra_force, ONLY: force_intra_control
49 : USE fist_neighbor_list_control, ONLY: list_control
50 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
51 : USE fist_nonbond_force, ONLY: bonded_correct_gaussian,&
52 : force_nonbond
53 : USE fist_pol_scf, ONLY: fist_pol_evaluate
54 : USE input_constants, ONLY: do_fist_pol_none
55 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
56 : section_vals_type
57 : USE kinds, ONLY: dp
58 : USE manybody_eam, ONLY: density_nonbond
59 : USE manybody_potential, ONLY: energy_manybody,&
60 : force_nonbond_manybody
61 : USE message_passing, ONLY: mp_para_env_type
62 : USE molecule_kind_types, ONLY: molecule_kind_type
63 : USE molecule_types, ONLY: molecule_type
64 : USE multipole_types, ONLY: multipole_type
65 : USE particle_types, ONLY: particle_type
66 : USE pme, ONLY: pme_evaluate
67 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
68 : do_ewald_none,&
69 : do_ewald_pme,&
70 : do_ewald_spme
71 : USE shell_potential_types, ONLY: shell_kind_type
72 : USE spme, ONLY: spme_evaluate
73 : USE virial_types, ONLY: virial_type
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 : PRIVATE
78 : PUBLIC :: fist_calc_energy_force
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_force'
80 :
81 : TYPE debug_variables_type
82 : REAL(KIND=dp) :: pot_manybody, pot_bend, pot_torsion
83 : REAL(KIND=dp) :: pot_nonbond, pot_g, pot_bond
84 : REAL(KIND=dp) :: pot_imptors, pot_urey_bradley
85 : REAL(KIND=dp) :: pot_opbend
86 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_nonbond, f_g, f_bond, f_bend, &
87 : f_torsion, f_imptors, f_ub, &
88 : f_opbend
89 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_nonbond, pv_g, pv_bond, pv_ub, &
90 : pv_bend, pv_torsion, pv_imptors, &
91 : pv_opbend
92 : END TYPE debug_variables_type
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief Calculates the total potential energy, total force, and the
98 : !> total pressure tensor from the potentials
99 : !> \param fist_env ...
100 : !> \param debug ...
101 : !> \par History
102 : !> Harald Forbert(Dec-2000): Changes for multiple linked lists
103 : !> cjm, 20-Feb-2001: box_ref used to initialize ewald. Now
104 : !> have consistent restarts with npt and ewald
105 : !> JGH(15-Mar-2001): box_change replaces ensemble parameter
106 : !> Call ewald_setup if box has changed
107 : !> Consistent setup for PME and SPME
108 : !> cjm, 28-Feb-2006: box_change is gone
109 : !> \author CJM & JGH
110 : ! **************************************************************************************************
111 83088 : SUBROUTINE fist_calc_energy_force(fist_env, debug)
112 : TYPE(fist_environment_type), POINTER :: fist_env
113 : TYPE(debug_variables_type), INTENT(INOUT), &
114 : OPTIONAL :: debug
115 :
116 : CHARACTER(len=*), PARAMETER :: routineN = 'fist_calc_energy_force'
117 :
118 : INTEGER :: do_ipol, ewald_type, fg_coulomb_size, handle, i, ii, ikind, iparticle_kind, &
119 : iparticle_local, iw, iw2, j, natoms, nlocal_particles, node, nparticle_kind, &
120 : nparticle_local, nshell, shell_index
121 : LOGICAL :: do_multipoles, shell_model_ad, &
122 : shell_present, use_virial
123 : REAL(KIND=dp) :: ef_ener, fc, fs, mass, pot_bend, pot_bond, pot_imptors, pot_manybody, &
124 : pot_nonbond, pot_opbend, pot_shell, pot_torsion, pot_urey_bradley, vg_coulomb, xdum1, &
125 : xdum2, xdum3
126 83088 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ef_f, f_nonbond, f_total, fcore_nonbond, &
127 83088 : fcore_shell_bonded, fcore_total, fg_coulomb, fgcore_coulomb, fgshell_coulomb, &
128 83088 : fshell_nonbond, fshell_total
129 : REAL(KIND=dp), DIMENSION(3, 3) :: ef_pv, ident, pv_bc, pv_bend, pv_bond, &
130 : pv_g, pv_imptors, pv_nonbond, &
131 : pv_opbend, pv_torsion, pv_urey_bradley
132 83088 : REAL(KIND=dp), DIMENSION(:), POINTER :: e_coulomb
133 83088 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pv_coulomb
134 83088 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
135 : TYPE(atomic_kind_type), POINTER :: atomic_kind
136 : TYPE(atprop_type), POINTER :: atprop_env
137 : TYPE(cell_type), POINTER :: cell
138 : TYPE(cp_logger_type), POINTER :: logger
139 : TYPE(cp_subsys_type), POINTER :: subsys
140 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
141 : TYPE(ewald_environment_type), POINTER :: ewald_env
142 : TYPE(ewald_pw_type), POINTER :: ewald_pw
143 83088 : TYPE(exclusion_type), DIMENSION(:), POINTER :: exclusions
144 : TYPE(fist_efield_type), POINTER :: efield
145 : TYPE(fist_energy_type), POINTER :: thermo
146 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
147 83088 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
148 83088 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
149 : TYPE(mp_para_env_type), POINTER :: para_env
150 : TYPE(multipole_type), POINTER :: multipoles
151 83088 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
152 83088 : shell_particle_set
153 : TYPE(section_vals_type), POINTER :: force_env_section, mm_section, &
154 : print_section
155 : TYPE(shell_kind_type), POINTER :: shell
156 : TYPE(virial_type), POINTER :: virial
157 :
158 83088 : CALL timeset(routineN, handle)
159 83088 : NULLIFY (logger)
160 83088 : NULLIFY (subsys, virial, atprop_env, para_env, force_env_section)
161 83088 : logger => cp_get_default_logger()
162 :
163 : CALL fist_env_get(fist_env, &
164 : subsys=subsys, &
165 : para_env=para_env, &
166 83088 : input=force_env_section)
167 :
168 : CALL cp_subsys_get(subsys, &
169 : virial=virial, &
170 83088 : atprop=atprop_env)
171 :
172 83088 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
173 83088 : NULLIFY (atomic_kind, atomic_kind_set, cell, ewald_pw, ewald_env, &
174 83088 : fist_nonbond_env, mm_section, local_molecules, local_particles, &
175 83088 : molecule_kind_set, molecule_set, particle_set, print_section, &
176 83088 : shell, shell_particle_set, core_particle_set, thermo, multipoles, &
177 83088 : e_coulomb, pv_coulomb)
178 :
179 83088 : mm_section => section_vals_get_subs_vals(force_env_section, "MM")
180 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%DERIVATIVES", &
181 83088 : extension=".mmLog")
182 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
183 83088 : extension=".mmLog")
184 :
185 : CALL fist_env_get(fist_env, ewald_pw=ewald_pw, ewald_env=ewald_env, &
186 : local_particles=local_particles, particle_set=particle_set, &
187 : atomic_kind_set=atomic_kind_set, molecule_set=molecule_set, &
188 : local_molecules=local_molecules, thermo=thermo, &
189 : molecule_kind_set=molecule_kind_set, fist_nonbond_env=fist_nonbond_env, &
190 : cell=cell, shell_model=shell_present, shell_model_ad=shell_model_ad, &
191 83088 : multipoles=multipoles, exclusions=exclusions, efield=efield)
192 :
193 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
194 83088 : do_ipol=do_ipol)
195 :
196 : ! Initialize ewald grids
197 83088 : CALL init_cell(cell)
198 83088 : CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
199 :
200 83088 : natoms = SIZE(particle_set)
201 83088 : nlocal_particles = 0
202 83088 : nparticle_kind = SIZE(atomic_kind_set)
203 316869 : DO ikind = 1, nparticle_kind
204 316869 : nlocal_particles = nlocal_particles + local_particles%n_el(ikind)
205 : END DO
206 :
207 249264 : ALLOCATE (f_nonbond(3, natoms))
208 33098888 : f_nonbond = 0.0_dp
209 :
210 83088 : nshell = 0
211 83088 : IF (shell_present) THEN
212 : CALL fist_env_get(fist_env, shell_particle_set=shell_particle_set, &
213 9452 : core_particle_set=core_particle_set)
214 9452 : CPASSERT(ASSOCIATED(shell_particle_set))
215 9452 : CPASSERT(ASSOCIATED(core_particle_set))
216 9452 : nshell = SIZE(shell_particle_set)
217 28356 : ALLOCATE (fshell_nonbond(3, nshell))
218 3058388 : fshell_nonbond = 0.0_dp
219 28356 : ALLOCATE (fcore_nonbond(3, nshell))
220 3058388 : fcore_nonbond = 0.0_dp
221 : ELSE
222 73636 : NULLIFY (shell_particle_set, core_particle_set)
223 : END IF
224 :
225 83088 : IF (fist_nonbond_env%do_nonbonded) THEN
226 : ! First check with list_control to update neighbor lists
227 82936 : IF (ASSOCIATED(exclusions)) THEN
228 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
229 : cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
230 78790 : core_particle_set, exclusions=exclusions)
231 : ELSE
232 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
233 : cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
234 4146 : core_particle_set)
235 : END IF
236 : END IF
237 :
238 : ! Initialize force, energy and pressure tensor arrays
239 8337038 : DO i = 1, natoms
240 8253950 : particle_set(i)%f(1) = 0.0_dp
241 8253950 : particle_set(i)%f(2) = 0.0_dp
242 8337038 : particle_set(i)%f(3) = 0.0_dp
243 : END DO
244 83088 : IF (nshell > 0) THEN
245 771686 : DO i = 1, nshell
246 762234 : shell_particle_set(i)%f(1) = 0.0_dp
247 762234 : shell_particle_set(i)%f(2) = 0.0_dp
248 762234 : shell_particle_set(i)%f(3) = 0.0_dp
249 762234 : core_particle_set(i)%f(1) = 0.0_dp
250 762234 : core_particle_set(i)%f(2) = 0.0_dp
251 771686 : core_particle_set(i)%f(3) = 0.0_dp
252 : END DO
253 : END IF
254 :
255 83088 : IF (use_virial) THEN
256 17102 : pv_bc = 0.0_dp
257 17102 : pv_bond = 0.0_dp
258 17102 : pv_bend = 0.0_dp
259 17102 : pv_torsion = 0.0_dp
260 17102 : pv_imptors = 0.0_dp
261 17102 : pv_opbend = 0.0_dp
262 17102 : pv_urey_bradley = 0.0_dp
263 17102 : pv_nonbond = 0.0_dp
264 17102 : pv_g = 0.0_dp
265 222326 : virial%pv_virial = 0.0_dp
266 : END IF
267 :
268 83088 : pot_nonbond = 0.0_dp
269 83088 : pot_manybody = 0.0_dp
270 83088 : pot_bond = 0.0_dp
271 83088 : pot_bend = 0.0_dp
272 83088 : pot_torsion = 0.0_dp
273 83088 : pot_opbend = 0.0_dp
274 83088 : pot_imptors = 0.0_dp
275 83088 : pot_urey_bradley = 0.0_dp
276 83088 : pot_shell = 0.0_dp
277 83088 : vg_coulomb = 0.0_dp
278 83088 : thermo%pot = 0.0_dp
279 83088 : thermo%harm_shell = 0.0_dp
280 :
281 : ! Get real-space non-bonded forces:
282 83088 : IF (iw > 0) THEN
283 285 : WRITE (iw, '(A)') " FIST:: FORCES IN INPUT..."
284 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
285 : END IF
286 :
287 83088 : IF (fist_nonbond_env%do_nonbonded) THEN
288 : ! Compute density for EAM
289 82936 : CALL density_nonbond(fist_nonbond_env, particle_set, cell, para_env)
290 :
291 : ! Compute embedding function and manybody energy
292 : CALL energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, &
293 82936 : cell, pot_manybody, para_env, mm_section)
294 :
295 : ! Nonbond contribution + Manybody Forces
296 82936 : IF (shell_present) THEN
297 : CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
298 : pot_nonbond, f_nonbond, pv_nonbond, &
299 : fshell_nonbond=fshell_nonbond, fcore_nonbond=fcore_nonbond, &
300 : atprop_env=atprop_env, &
301 9364 : atomic_kind_set=atomic_kind_set, use_virial=use_virial)
302 : ELSE
303 : CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
304 : pot_nonbond, f_nonbond, pv_nonbond, atprop_env=atprop_env, &
305 73572 : atomic_kind_set=atomic_kind_set, use_virial=use_virial)
306 : CALL force_nonbond_manybody(fist_nonbond_env, particle_set, cell, f_nonbond, pv_nonbond, &
307 73572 : use_virial=use_virial)
308 : END IF
309 : END IF
310 :
311 83088 : IF (iw > 0) THEN
312 285 : WRITE (iw, '(A)') " FIST:: NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
313 285 : WRITE (iw, '(3f15.9)') f_nonbond
314 285 : IF (shell_present .AND. shell_model_ad) THEN
315 44 : WRITE (iw, '(A)') " FIST:: SHELL NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
316 44 : WRITE (iw, '(3f15.9)') fshell_nonbond
317 44 : WRITE (iw, '(A)') " FIST:: CORE NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
318 44 : WRITE (iw, '(3f15.9)') fcore_nonbond
319 : END IF
320 : END IF
321 :
322 : ! Get g-space non-bonded forces:
323 83088 : IF (ewald_type /= do_ewald_none) THEN
324 : ! Determine size of the forces array
325 29690 : SELECT CASE (ewald_type)
326 : CASE (do_ewald_ewald)
327 29690 : fg_coulomb_size = nlocal_particles
328 : CASE DEFAULT
329 60992 : fg_coulomb_size = natoms
330 : END SELECT
331 : ! Allocate and zeroing arrays
332 182430 : ALLOCATE (fg_coulomb(3, fg_coulomb_size))
333 27656076 : fg_coulomb = 0.0_dp
334 60992 : IF (shell_present) THEN
335 28314 : ALLOCATE (fgshell_coulomb(3, nshell))
336 28314 : ALLOCATE (fgcore_coulomb(3, nshell))
337 3058294 : fgshell_coulomb = 0.0_dp
338 3058294 : fgcore_coulomb = 0.0_dp
339 : END IF
340 60992 : IF (shell_present .AND. do_multipoles) THEN
341 0 : CPABORT("Multipoles and Core-Shell model not implemented.")
342 : END IF
343 : ! If not multipole: Compute self-interaction and neutralizing background
344 : ! for multipoles is handled separately..
345 60992 : IF (.NOT. do_multipoles) THEN
346 : CALL ewald_self(ewald_env, cell, atomic_kind_set, local_particles, &
347 59834 : thermo%e_self, thermo%e_neut, fist_nonbond_env%charges)
348 59834 : IF (atprop_env%energy) THEN
349 : CALL ewald_self_atom(ewald_env, atomic_kind_set, local_particles, &
350 636 : atprop_env%atener, fist_nonbond_env%charges)
351 6054 : atprop_env%atener = atprop_env%atener + thermo%e_neut/SIZE(atprop_env%atener)
352 : END IF
353 : END IF
354 :
355 : ! Polarizable force-field
356 60992 : IF (do_ipol /= do_fist_pol_none) THEN
357 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
358 104 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
359 0 : CPABORT("Polarizable force field and array charges not implemented!")
360 : END IF
361 : ! Converge the dipoles self-consistently
362 : CALL fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
363 : ewald_pw, fist_nonbond_env, cell, particle_set, &
364 : local_particles, thermo, vg_coulomb, pot_nonbond, f_nonbond, &
365 104 : fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section, do_ipol)
366 : ELSE
367 : ! Non-Polarizable force-field
368 29586 : SELECT CASE (ewald_type)
369 : CASE (do_ewald_ewald)
370 : ! Parallelisation over atoms --> allocate local atoms
371 29586 : IF (shell_present) THEN
372 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
373 0 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
374 0 : CPABORT("Core-Shell and array charges not implemented!")
375 : END IF
376 0 : IF (do_multipoles) THEN
377 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within Ewald sum!")
378 : ELSE
379 0 : CPABORT("Core-Shell model not yet implemented within Ewald sums.")
380 : END IF
381 : ELSE
382 29586 : IF (do_multipoles) THEN
383 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
384 1054 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
385 0 : CPABORT("Multipole Ewald and array charges not implemented!")
386 : END IF
387 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
388 : particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
389 : thermo%e_self, multipoles%task, do_correction_bonded=.TRUE., do_forces=.TRUE., &
390 : do_stress=use_virial, do_efield=.FALSE., radii=multipoles%radii, &
391 : charges=multipoles%charges, dipoles=multipoles%dipoles, &
392 : quadrupoles=multipoles%quadrupoles, forces_local=fg_coulomb, &
393 : forces_glob=f_nonbond, pv_local=pv_g, pv_glob=pv_nonbond, iw=iw2, &
394 1054 : do_debug=.TRUE., atomic_kind_set=atomic_kind_set, mm_section=mm_section)
395 : ELSE
396 28532 : IF (atprop_env%energy) THEN
397 505 : ALLOCATE (e_coulomb(fg_coulomb_size))
398 : END IF
399 28532 : IF (atprop_env%stress) THEN
400 505 : ALLOCATE (pv_coulomb(3, 3, fg_coulomb_size))
401 : END IF
402 : CALL ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
403 : local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial=use_virial, &
404 : charges=fist_nonbond_env%charges, e_coulomb=e_coulomb, &
405 28532 : pv_coulomb=pv_coulomb)
406 : END IF
407 : END IF
408 : CASE (do_ewald_pme)
409 : ! Parallelisation over grids --> allocate all atoms
410 1818 : IF (shell_present) THEN
411 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
412 0 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
413 0 : CPABORT("Core-Shell and array charges not implemented!")
414 : END IF
415 0 : IF (do_multipoles) THEN
416 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a PME scheme!")
417 : ELSE
418 : CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
419 : pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
420 : fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
421 0 : atprop=atprop_env)
422 0 : CALL para_env%sum(fgshell_coulomb)
423 0 : CALL para_env%sum(fgcore_coulomb)
424 : END IF
425 : ELSE
426 1818 : IF (do_multipoles) THEN
427 0 : CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
428 : ELSE
429 : CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
430 : pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
431 1818 : atprop=atprop_env)
432 : END IF
433 : END IF
434 1818 : CALL para_env%sum(fg_coulomb)
435 : CASE (do_ewald_spme)
436 : ! Parallelisation over grids --> allocate all atoms
437 29484 : IF (shell_present) THEN
438 : ! Check if an array of charges was provided and in case abort due to lack of implementation
439 9438 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
440 0 : CPABORT("Core-Shell and array charges not implemented!")
441 : END IF
442 9438 : IF (do_multipoles) THEN
443 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a SPME scheme!")
444 : ELSE
445 : CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
446 : pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
447 : fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
448 9438 : atprop=atprop_env)
449 9438 : CALL para_env%sum(fgshell_coulomb)
450 9438 : CALL para_env%sum(fgcore_coulomb)
451 : END IF
452 : ELSE
453 20046 : IF (do_multipoles) THEN
454 0 : CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
455 : ELSE
456 : CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
457 : pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
458 20046 : atprop=atprop_env)
459 : END IF
460 : END IF
461 90372 : CALL para_env%sum(fg_coulomb)
462 : END SELECT
463 : END IF
464 :
465 : ! Subtract the interaction between screening charges. This is a
466 : ! correction in real-space and depends on the neighborlists. Therefore
467 : ! it is only executed if fist_nonbond_env%do_nonbonded is set.
468 60992 : IF (fist_nonbond_env%do_nonbonded) THEN
469 : ! Correction for core-shell model
470 60904 : IF (shell_present) THEN
471 : CALL bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
472 : local_particles, particle_set, ewald_env, thermo%e_bonded, &
473 : pv_bc, shell_particle_set=shell_particle_set, &
474 : core_particle_set=core_particle_set, atprop_env=atprop_env, cell=cell, &
475 9350 : use_virial=use_virial)
476 : ELSE
477 51554 : IF (.NOT. do_multipoles) THEN
478 : CALL bonded_correct_gaussian(fist_nonbond_env, &
479 : atomic_kind_set, local_particles, particle_set, &
480 : ewald_env, thermo%e_bonded, pv_bc=pv_bc, atprop_env=atprop_env, cell=cell, &
481 50396 : use_virial=use_virial)
482 : END IF
483 : END IF
484 : END IF
485 :
486 60992 : IF (.NOT. do_multipoles) THEN
487 : ! Multipole code has its own printing routine.
488 : CALL ewald_print(iw2, pot_nonbond, vg_coulomb, thermo%e_self, thermo%e_neut, &
489 59834 : thermo%e_bonded)
490 : END IF
491 : ELSE
492 22096 : IF (use_virial) THEN
493 3178 : pv_g = 0.0_dp
494 3178 : pv_bc = 0.0_dp
495 : END IF
496 22096 : thermo%e_neut = 0.0_dp
497 : END IF
498 :
499 83088 : IF (iw > 0) THEN
500 285 : IF (ALLOCATED(fg_coulomb)) THEN
501 206 : WRITE (iw, '(A)') " FIST:: NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
502 206 : WRITE (iw, '(3f15.9)') ((fg_coulomb(j, i), j=1, 3), i=1, SIZE(fg_coulomb, 2))
503 : END IF
504 285 : IF (shell_present .AND. shell_model_ad .AND. ALLOCATED(fgshell_coulomb)) THEN
505 44 : WRITE (iw, '(A)') " FIST:: SHELL NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
506 44 : WRITE (iw, '(3f15.9)') ((fgshell_coulomb(j, i), j=1, 3), i=1, SIZE(fgshell_coulomb, 2))
507 44 : WRITE (iw, '(A)') " FIST:: CORE NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
508 44 : WRITE (iw, '(3f15.9)') ((fgcore_coulomb(j, i), j=1, 3), i=1, SIZE(fgcore_coulomb, 2))
509 : END IF
510 : END IF
511 :
512 : ! calculate the action of an (external) electric field
513 83088 : IF (efield%apply_field) THEN
514 348 : ALLOCATE (ef_f(3, natoms))
515 : CALL fist_efield_energy_force(ef_ener, ef_f, ef_pv, atomic_kind_set, particle_set, cell, &
516 116 : efield, use_virial=use_virial, iunit=iw, charges=fist_nonbond_env%charges)
517 : END IF
518 :
519 : ! Get intramolecular forces
520 83088 : IF (PRESENT(debug)) THEN
521 : CALL force_intra_control(molecule_set, molecule_kind_set, &
522 : local_molecules, particle_set, shell_particle_set, &
523 : core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
524 : pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
525 : pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
526 : debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, &
527 0 : debug%f_imptors, debug%f_opbend, cell, use_virial, atprop_env)
528 :
529 : ELSE
530 : CALL force_intra_control(molecule_set, molecule_kind_set, &
531 : local_molecules, particle_set, shell_particle_set, &
532 : core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
533 : pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
534 : pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
535 83088 : cell=cell, use_virial=use_virial, atprop_env=atprop_env)
536 : END IF
537 :
538 : ! Perform global sum of the intra-molecular (bonded) forces for the core-shell atoms
539 83088 : IF (shell_present .AND. shell_model_ad) THEN
540 28356 : ALLOCATE (fcore_shell_bonded(3, nshell))
541 3058388 : fcore_shell_bonded = 0.0_dp
542 806318 : DO i = 1, natoms
543 796866 : shell_index = particle_set(i)%shell_index
544 806318 : IF (shell_index /= 0) THEN
545 3048936 : fcore_shell_bonded(1:3, shell_index) = particle_set(i)%f(1:3)
546 : END IF
547 : END DO
548 9452 : CALL para_env%sum(fcore_shell_bonded)
549 806318 : DO i = 1, natoms
550 796866 : shell_index = particle_set(i)%shell_index
551 806318 : IF (shell_index /= 0) THEN
552 3048936 : particle_set(i)%f(1:3) = fcore_shell_bonded(1:3, shell_index)
553 : END IF
554 : END DO
555 9452 : DEALLOCATE (fcore_shell_bonded)
556 : END IF
557 :
558 83088 : IF (iw > 0) THEN
559 285 : xdum1 = cp_unit_from_cp2k(pot_bond, "kcalmol")
560 285 : xdum2 = cp_unit_from_cp2k(pot_bend, "kcalmol")
561 285 : xdum3 = cp_unit_from_cp2k(pot_urey_bradley, "kcalmol")
562 285 : WRITE (iw, '(A)') " FIST energy contributions in kcal/mol:"
563 : WRITE (iw, '(1x,"BOND = ",f13.4,'// &
564 : '2x,"ANGLE = ",f13.4,'// &
565 285 : '2x,"UBRAD = ",f13.4)') xdum1, xdum2, xdum3
566 285 : xdum1 = cp_unit_from_cp2k(pot_torsion, "kcalmol")
567 285 : xdum2 = cp_unit_from_cp2k(pot_imptors, "kcalmol")
568 285 : xdum3 = cp_unit_from_cp2k(pot_opbend, "kcalmol")
569 : WRITE (iw, '(1x,"TORSION = ",f13.4,'// &
570 : '2x,"IMPTORS = ",f13.4,'// &
571 285 : '2x,"OPBEND = ",f13.4)') xdum1, xdum2, xdum3
572 :
573 285 : WRITE (iw, '(A)') " FIST:: CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
574 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
575 285 : IF (shell_present .AND. shell_model_ad) THEN
576 44 : WRITE (iw, '(A)') " FIST:: SHELL CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
577 16940 : WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, SIZE(shell_particle_set))
578 44 : WRITE (iw, '(A)') " FIST:: CORE CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
579 16940 : WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, SIZE(core_particle_set))
580 : END IF
581 : END IF
582 :
583 : ! add up all the potential energies
584 : thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_opbend + &
585 83088 : pot_imptors + pot_urey_bradley + pot_manybody + pot_shell
586 :
587 83088 : CALL para_env%sum(thermo%pot)
588 :
589 83088 : IF (shell_present) THEN
590 9452 : thermo%harm_shell = pot_shell
591 9452 : CALL para_env%sum(thermo%harm_shell)
592 : END IF
593 : ! add g-space contributions if needed
594 83088 : IF (ewald_type /= do_ewald_none) THEN
595 : ! e_self, e_neut, and ebonded are already summed over all processors
596 : ! vg_coulomb is not calculated in parallel
597 60992 : thermo%e_gspace = vg_coulomb
598 60992 : thermo%pot = thermo%pot + thermo%e_self + thermo%e_neut
599 60992 : thermo%pot = thermo%pot + vg_coulomb + thermo%e_bonded
600 : ! add the induction energy of the dipoles for polarizable models
601 60992 : IF (do_ipol /= do_fist_pol_none) thermo%pot = thermo%pot + thermo%e_induction
602 : END IF
603 :
604 : ! add up all the forces
605 : !
606 : ! nonbonded forces might be calculated for atoms not on this node
607 : ! ewald forces are strictly local -> sum only over pnode
608 : ! We first sum the forces in f_nonbond, this allows for a more efficient
609 : ! global sum in the parallel code and in the end copy them back to part
610 249264 : ALLOCATE (f_total(3, natoms))
611 33098888 : f_total = 0.0_dp
612 8337038 : DO i = 1, natoms
613 8253950 : f_total(1, i) = particle_set(i)%f(1) + f_nonbond(1, i)
614 8253950 : f_total(2, i) = particle_set(i)%f(2) + f_nonbond(2, i)
615 8337038 : f_total(3, i) = particle_set(i)%f(3) + f_nonbond(3, i)
616 : END DO
617 83088 : IF (shell_present) THEN
618 28356 : ALLOCATE (fshell_total(3, nshell))
619 3058388 : fshell_total = 0.0_dp
620 28356 : ALLOCATE (fcore_total(3, nshell))
621 3058388 : fcore_total = 0.0_dp
622 771686 : DO i = 1, nshell
623 762234 : fcore_total(1, i) = core_particle_set(i)%f(1) + fcore_nonbond(1, i)
624 762234 : fcore_total(2, i) = core_particle_set(i)%f(2) + fcore_nonbond(2, i)
625 762234 : fcore_total(3, i) = core_particle_set(i)%f(3) + fcore_nonbond(3, i)
626 762234 : fshell_total(1, i) = shell_particle_set(i)%f(1) + fshell_nonbond(1, i)
627 762234 : fshell_total(2, i) = shell_particle_set(i)%f(2) + fshell_nonbond(2, i)
628 771686 : fshell_total(3, i) = shell_particle_set(i)%f(3) + fshell_nonbond(3, i)
629 : END DO
630 : END IF
631 :
632 83088 : IF (iw > 0) THEN
633 285 : WRITE (iw, '(A)') " FIST::(1)INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
634 285 : WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
635 285 : IF (shell_present .AND. shell_model_ad) THEN
636 44 : WRITE (iw, '(A)') " FIST::(1)SHELL INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
637 44 : WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
638 44 : WRITE (iw, '(A)') " FIST::(1)CORE INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
639 44 : WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
640 : END IF
641 : END IF
642 :
643 : ! Adding in the reciprocal forces: EWALD is a special case because of distrubted data
644 83088 : IF (ewald_type == do_ewald_ewald) THEN
645 : node = 0
646 131012 : DO iparticle_kind = 1, nparticle_kind
647 101322 : nparticle_local = local_particles%n_el(iparticle_kind)
648 558785 : DO iparticle_local = 1, nparticle_local
649 427773 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
650 427773 : node = node + 1
651 427773 : f_total(1, ii) = f_total(1, ii) + fg_coulomb(1, node)
652 427773 : f_total(2, ii) = f_total(2, ii) + fg_coulomb(2, node)
653 427773 : f_total(3, ii) = f_total(3, ii) + fg_coulomb(3, node)
654 427773 : IF (PRESENT(debug)) THEN
655 0 : debug%f_g(1, ii) = debug%f_g(1, ii) + fg_coulomb(1, node)
656 0 : debug%f_g(2, ii) = debug%f_g(2, ii) + fg_coulomb(2, node)
657 0 : debug%f_g(3, ii) = debug%f_g(3, ii) + fg_coulomb(3, node)
658 : END IF
659 427773 : IF (atprop_env%energy) THEN
660 303 : atprop_env%atener(ii) = atprop_env%atener(ii) + e_coulomb(node)
661 : END IF
662 529095 : IF (atprop_env%stress) THEN
663 7878 : atprop_env%atstress(:, :, ii) = atprop_env%atstress(:, :, ii) + pv_coulomb(:, :, node)
664 : END IF
665 : END DO
666 : END DO
667 29690 : IF (atprop_env%energy) THEN
668 202 : DEALLOCATE (e_coulomb)
669 : END IF
670 29690 : IF (atprop_env%stress) THEN
671 202 : DEALLOCATE (pv_coulomb)
672 : END IF
673 : END IF
674 :
675 83088 : IF (iw > 0) THEN
676 285 : WRITE (iw, '(A)') " FIST::(2)TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
677 285 : WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
678 285 : IF (shell_present .AND. shell_model_ad) THEN
679 44 : WRITE (iw, '(A)') " FIST::(2)SHELL TOTAL FORCES(1)+ ELECTROSTATIC FORCES "
680 44 : WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
681 44 : WRITE (iw, '(A)') " FIST::(2)CORE TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
682 44 : WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
683 : END IF
684 : END IF
685 :
686 83088 : IF (use_virial) THEN
687 : ! Add up all the pressure tensors
688 17102 : IF (ewald_type == do_ewald_none) THEN
689 41314 : virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley
690 3178 : CALL para_env%sum(virial%pv_virial)
691 : ELSE
692 13924 : ident = 0.0_dp
693 55696 : DO i = 1, 3
694 55696 : ident(i, i) = 1.0_dp
695 : END DO
696 :
697 181012 : virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc
698 13924 : CALL para_env%sum(virial%pv_virial)
699 :
700 181012 : virial%pv_virial = virial%pv_virial + ident*thermo%e_neut
701 181012 : virial%pv_virial = virial%pv_virial + pv_g
702 : END IF
703 : END IF
704 :
705 : ! Sum total forces
706 83088 : CALL para_env%sum(f_total)
707 83088 : IF (shell_present .AND. shell_model_ad) THEN
708 9452 : CALL para_env%sum(fcore_total)
709 9452 : CALL para_env%sum(fshell_total)
710 : END IF
711 :
712 : ! contributions from fields (currently all quantities are fully replicated!)
713 83088 : IF (efield%apply_field) THEN
714 116 : thermo%pot = thermo%pot + ef_ener
715 1508 : f_total(1:3, 1:natoms) = f_total(1:3, 1:natoms) + ef_f(1:3, 1:natoms)
716 116 : IF (use_virial) THEN
717 26 : virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + ef_pv(1:3, 1:3)
718 : END IF
719 116 : DEALLOCATE (ef_f)
720 : END IF
721 :
722 : ! Assign to particle_set
723 31302 : SELECT CASE (ewald_type)
724 : CASE (do_ewald_spme, do_ewald_pme)
725 31302 : IF (shell_present .AND. shell_model_ad) THEN
726 806276 : DO i = 1, natoms
727 796838 : shell_index = particle_set(i)%shell_index
728 806276 : IF (shell_index == 0) THEN
729 138496 : particle_set(i)%f(1:3) = f_total(1:3, i) + fg_coulomb(1:3, i)
730 : ELSE
731 762214 : atomic_kind => particle_set(i)%atomic_kind
732 762214 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
733 762214 : fc = shell%mass_core/mass
734 3048856 : fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
735 762214 : fs = shell%mass_shell/mass
736 3048856 : fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
737 : END IF
738 : END DO
739 :
740 771652 : DO i = 1, nshell
741 762214 : core_particle_set(i)%f(1) = fcore_total(1, i) + fgcore_coulomb(1, i)
742 762214 : core_particle_set(i)%f(2) = fcore_total(2, i) + fgcore_coulomb(2, i)
743 762214 : core_particle_set(i)%f(3) = fcore_total(3, i) + fgcore_coulomb(3, i)
744 762214 : shell_particle_set(i)%f(1) = fshell_total(1, i) + fgshell_coulomb(1, i)
745 762214 : shell_particle_set(i)%f(2) = fshell_total(2, i) + fgshell_coulomb(2, i)
746 771652 : shell_particle_set(i)%f(3) = fshell_total(3, i) + fgshell_coulomb(3, i)
747 : END DO
748 :
749 21864 : ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
750 0 : CPABORT("Non adiabatic shell-model not implemented.")
751 : ELSE
752 5696024 : DO i = 1, natoms
753 5674160 : particle_set(i)%f(1) = f_total(1, i) + fg_coulomb(1, i)
754 5674160 : particle_set(i)%f(2) = f_total(2, i) + fg_coulomb(2, i)
755 5696024 : particle_set(i)%f(3) = f_total(3, i) + fg_coulomb(3, i)
756 : END DO
757 : END IF
758 : CASE (do_ewald_ewald, do_ewald_none)
759 83088 : IF (shell_present .AND. shell_model_ad) THEN
760 42 : DO i = 1, natoms
761 28 : shell_index = particle_set(i)%shell_index
762 42 : IF (shell_index == 0) THEN
763 32 : particle_set(i)%f(1:3) = f_total(1:3, i)
764 : ELSE
765 20 : atomic_kind => particle_set(i)%atomic_kind
766 20 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
767 20 : fc = shell%mass_core/mass
768 80 : fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
769 20 : fs = shell%mass_shell/mass
770 80 : fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
771 : END IF
772 : END DO
773 34 : DO i = 1, nshell
774 20 : core_particle_set(i)%f(1) = fcore_total(1, i)
775 20 : core_particle_set(i)%f(2) = fcore_total(2, i)
776 20 : core_particle_set(i)%f(3) = fcore_total(3, i)
777 20 : shell_particle_set(i)%f(1) = fshell_total(1, i)
778 20 : shell_particle_set(i)%f(2) = fshell_total(2, i)
779 34 : shell_particle_set(i)%f(3) = fshell_total(3, i)
780 : END DO
781 51772 : ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
782 0 : CPABORT("Non adiabatic shell-model not implemented.")
783 : ELSE
784 1834696 : DO i = 1, natoms
785 1782924 : particle_set(i)%f(1) = f_total(1, i)
786 1782924 : particle_set(i)%f(2) = f_total(2, i)
787 1834696 : particle_set(i)%f(3) = f_total(3, i)
788 : END DO
789 : END IF
790 : END SELECT
791 :
792 83088 : IF (iw > 0) THEN
793 285 : WRITE (iw, '(A)') " FIST::(3)TOTAL FORCES - THE END..."
794 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, natoms)
795 285 : IF (shell_present .AND. shell_model_ad) THEN
796 44 : WRITE (iw, '(A)') " FIST::(3)SHELL TOTAL FORCES - THE END..."
797 16940 : WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, nshell)
798 44 : WRITE (iw, '(A)') " FIST::(3)CORE TOTAL FORCES - THE END..."
799 16940 : WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, nshell)
800 : END IF
801 285 : WRITE (iw, '(A,f15.9)') "Energy after FIST calculation.. exiting now ::", thermo%pot
802 : END IF
803 : !
804 : ! If we are doing debugging, check if variables are present and assign
805 : !
806 83088 : IF (PRESENT(debug)) THEN
807 0 : CALL para_env%sum(pot_manybody)
808 0 : debug%pot_manybody = pot_manybody
809 0 : CALL para_env%sum(pot_nonbond)
810 0 : debug%pot_nonbond = pot_nonbond
811 0 : CALL para_env%sum(pot_bond)
812 0 : debug%pot_bond = pot_bond
813 0 : CALL para_env%sum(pot_bend)
814 0 : debug%pot_bend = pot_bend
815 0 : CALL para_env%sum(pot_torsion)
816 0 : debug%pot_torsion = pot_torsion
817 0 : CALL para_env%sum(pot_imptors)
818 0 : debug%pot_imptors = pot_imptors
819 0 : CALL para_env%sum(pot_opbend)
820 0 : debug%pot_opbend = pot_opbend
821 0 : CALL para_env%sum(pot_urey_bradley)
822 0 : debug%pot_urey_bradley = pot_urey_bradley
823 0 : CALL para_env%sum(f_nonbond)
824 0 : debug%f_nonbond = f_nonbond
825 0 : IF (use_virial) THEN
826 0 : CALL para_env%sum(pv_nonbond)
827 0 : debug%pv_nonbond = pv_nonbond
828 0 : CALL para_env%sum(pv_bond)
829 0 : debug%pv_bond = pv_bond
830 0 : CALL para_env%sum(pv_bend)
831 0 : debug%pv_bend = pv_bend
832 0 : CALL para_env%sum(pv_torsion)
833 0 : debug%pv_torsion = pv_torsion
834 0 : CALL para_env%sum(pv_imptors)
835 0 : debug%pv_imptors = pv_imptors
836 0 : CALL para_env%sum(pv_urey_bradley)
837 0 : debug%pv_ub = pv_urey_bradley
838 : END IF
839 0 : SELECT CASE (ewald_type)
840 : CASE (do_ewald_spme, do_ewald_pme)
841 0 : debug%pot_g = vg_coulomb
842 0 : debug%pv_g = pv_g
843 0 : debug%f_g = fg_coulomb
844 : CASE (do_ewald_ewald)
845 0 : debug%pot_g = vg_coulomb
846 0 : IF (use_virial) debug%pv_g = pv_g
847 : CASE default
848 0 : debug%pot_g = 0.0_dp
849 0 : debug%f_g = 0.0_dp
850 0 : IF (use_virial) debug%pv_g = 0.0_dp
851 : END SELECT
852 : END IF
853 :
854 : ! print properties if requested
855 83088 : print_section => section_vals_get_subs_vals(mm_section, "PRINT")
856 83088 : CALL print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
857 :
858 : ! deallocating all local variables
859 83088 : IF (ALLOCATED(fg_coulomb)) THEN
860 60992 : DEALLOCATE (fg_coulomb)
861 : END IF
862 83088 : IF (ALLOCATED(f_total)) THEN
863 83088 : DEALLOCATE (f_total)
864 : END IF
865 83088 : DEALLOCATE (f_nonbond)
866 83088 : IF (shell_present) THEN
867 9452 : DEALLOCATE (fshell_total)
868 9452 : IF (ewald_type /= do_ewald_none) THEN
869 9438 : DEALLOCATE (fgshell_coulomb)
870 : END IF
871 9452 : DEALLOCATE (fshell_nonbond)
872 : END IF
873 83088 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%DERIVATIVES")
874 83088 : CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%EWALD_INFO")
875 83088 : CALL timestop(handle)
876 :
877 249264 : END SUBROUTINE fist_calc_energy_force
878 :
879 : ! **************************************************************************************************
880 : !> \brief Print properties number according the requests in input file
881 : !> \param fist_env ...
882 : !> \param print_section ...
883 : !> \param atomic_kind_set ...
884 : !> \param particle_set ...
885 : !> \param cell ...
886 : !> \par History
887 : !> [01.2006] created
888 : !> \author Teodoro Laino
889 : ! **************************************************************************************************
890 83088 : SUBROUTINE print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
891 : TYPE(fist_environment_type), POINTER :: fist_env
892 : TYPE(section_vals_type), POINTER :: print_section
893 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
894 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
895 : TYPE(cell_type), POINTER :: cell
896 :
897 : INTEGER :: unit_nr
898 : TYPE(cp_logger_type), POINTER :: logger
899 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
900 : TYPE(section_vals_type), POINTER :: print_key
901 :
902 83088 : NULLIFY (logger, print_key, fist_nonbond_env)
903 83088 : logger => cp_get_default_logger()
904 83088 : print_key => section_vals_get_subs_vals(print_section, "dipole")
905 83088 : CALL fist_env_get(fist_env, fist_nonbond_env=fist_nonbond_env)
906 83088 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
907 : cp_p_file)) THEN
908 : unit_nr = cp_print_key_unit_nr(logger, print_section, "dipole", &
909 21102 : extension=".data", middle_name="MM_DIPOLE", log_filename=.FALSE.)
910 : CALL fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
911 21102 : cell, unit_nr, fist_nonbond_env%charges)
912 21102 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
913 : END IF
914 :
915 83088 : END SUBROUTINE print_fist
916 :
917 0 : END MODULE fist_force
|
