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 = 0.0_dp, pot_bend = 0.0_dp, pot_torsion = 0.0_dp
83 : REAL(KIND=dp) :: pot_nonbond = 0.0_dp, pot_g = 0.0_dp, pot_bond = 0.0_dp
84 : REAL(KIND=dp) :: pot_imptors = 0.0_dp, pot_urey_bradley = 0.0_dp
85 : REAL(KIND=dp) :: pot_opbend = 0.0_dp
86 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_nonbond => NULL(), f_g => NULL(), f_bond => NULL(), f_bend => NULL(), &
87 : f_torsion => NULL(), f_imptors => NULL(), f_ub => NULL(), &
88 : f_opbend => NULL()
89 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_nonbond = 0.0_dp, pv_g = 0.0_dp, pv_bond = 0.0_dp, pv_ub = 0.0_dp, &
90 : pv_bend = 0.0_dp, pv_torsion = 0.0_dp, pv_imptors = 0.0_dp, &
91 : pv_opbend = 0.0_dp
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 80943 : 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 80943 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ef_f, f_nonbond, f_total, fcore_nonbond, &
127 80943 : fcore_shell_bonded, fcore_total, fg_coulomb, fgcore_coulomb, fgshell_coulomb, &
128 80943 : 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 80943 : REAL(KIND=dp), DIMENSION(:), POINTER :: e_coulomb
133 80943 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
134 : TYPE(atomic_kind_type), POINTER :: atomic_kind
135 : TYPE(atprop_type), POINTER :: atprop_env
136 : TYPE(cell_type), POINTER :: cell
137 : TYPE(cp_logger_type), POINTER :: logger
138 : TYPE(cp_subsys_type), POINTER :: subsys
139 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
140 : TYPE(ewald_environment_type), POINTER :: ewald_env
141 : TYPE(ewald_pw_type), POINTER :: ewald_pw
142 80943 : TYPE(exclusion_type), DIMENSION(:), POINTER :: exclusions
143 : TYPE(fist_efield_type), POINTER :: efield
144 : TYPE(fist_energy_type), POINTER :: thermo
145 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
146 80943 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
147 80943 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
148 : TYPE(mp_para_env_type), POINTER :: para_env
149 : TYPE(multipole_type), POINTER :: multipoles
150 80943 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, &
151 80943 : shell_particle_set
152 : TYPE(section_vals_type), POINTER :: force_env_section, mm_section, &
153 : print_section
154 : TYPE(shell_kind_type), POINTER :: shell
155 : TYPE(virial_type), POINTER :: virial
156 :
157 80943 : CALL timeset(routineN, handle)
158 80943 : NULLIFY (logger)
159 80943 : NULLIFY (subsys, virial, atprop_env, para_env, force_env_section)
160 80943 : logger => cp_get_default_logger()
161 :
162 : CALL fist_env_get(fist_env, &
163 : subsys=subsys, &
164 : para_env=para_env, &
165 80943 : input=force_env_section)
166 :
167 : CALL cp_subsys_get(subsys, &
168 : virial=virial, &
169 80943 : atprop=atprop_env)
170 :
171 80943 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
172 80943 : NULLIFY (atomic_kind, atomic_kind_set, cell, ewald_pw, ewald_env, &
173 80943 : fist_nonbond_env, mm_section, local_molecules, local_particles, &
174 80943 : molecule_kind_set, molecule_set, particle_set, print_section, &
175 80943 : shell, shell_particle_set, core_particle_set, thermo, multipoles, &
176 80943 : e_coulomb)
177 :
178 80943 : mm_section => section_vals_get_subs_vals(force_env_section, "MM")
179 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%DERIVATIVES", &
180 80943 : extension=".mmLog")
181 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
182 80943 : extension=".mmLog")
183 :
184 : CALL fist_env_get(fist_env, ewald_pw=ewald_pw, ewald_env=ewald_env, &
185 : local_particles=local_particles, particle_set=particle_set, &
186 : atomic_kind_set=atomic_kind_set, molecule_set=molecule_set, &
187 : local_molecules=local_molecules, thermo=thermo, &
188 : molecule_kind_set=molecule_kind_set, fist_nonbond_env=fist_nonbond_env, &
189 : cell=cell, shell_model=shell_present, shell_model_ad=shell_model_ad, &
190 80943 : multipoles=multipoles, exclusions=exclusions, efield=efield)
191 :
192 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
193 80943 : do_ipol=do_ipol)
194 :
195 : ! Initialize ewald grids
196 80943 : CALL init_cell(cell)
197 80943 : CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
198 :
199 80943 : natoms = SIZE(particle_set)
200 80943 : nlocal_particles = 0
201 80943 : nparticle_kind = SIZE(atomic_kind_set)
202 312620 : DO ikind = 1, nparticle_kind
203 312620 : nlocal_particles = nlocal_particles + local_particles%n_el(ikind)
204 : END DO
205 :
206 242829 : ALLOCATE (f_nonbond(3, natoms))
207 32985667 : f_nonbond = 0.0_dp
208 :
209 80943 : nshell = 0
210 80943 : IF (shell_present) THEN
211 : CALL fist_env_get(fist_env, shell_particle_set=shell_particle_set, &
212 9432 : core_particle_set=core_particle_set)
213 9432 : CPASSERT(ASSOCIATED(shell_particle_set))
214 9432 : CPASSERT(ASSOCIATED(core_particle_set))
215 9432 : nshell = SIZE(shell_particle_set)
216 28296 : ALLOCATE (fshell_nonbond(3, nshell))
217 3054432 : fshell_nonbond = 0.0_dp
218 28296 : ALLOCATE (fcore_nonbond(3, nshell))
219 3054432 : fcore_nonbond = 0.0_dp
220 : ELSE
221 71511 : NULLIFY (shell_particle_set, core_particle_set)
222 : END IF
223 :
224 80943 : IF (fist_nonbond_env%do_nonbonded) THEN
225 : ! First check with list_control to update neighbor lists
226 80791 : IF (ASSOCIATED(exclusions)) THEN
227 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
228 : cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
229 76657 : core_particle_set, exclusions=exclusions)
230 : ELSE
231 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
232 : cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
233 4134 : core_particle_set)
234 : END IF
235 : END IF
236 :
237 : ! Initialize force, energy and pressure tensor arrays
238 8307124 : DO i = 1, natoms
239 8226181 : particle_set(i)%f(1) = 0.0_dp
240 8226181 : particle_set(i)%f(2) = 0.0_dp
241 8307124 : particle_set(i)%f(3) = 0.0_dp
242 : END DO
243 80943 : IF (nshell > 0) THEN
244 770682 : DO i = 1, nshell
245 761250 : shell_particle_set(i)%f(1) = 0.0_dp
246 761250 : shell_particle_set(i)%f(2) = 0.0_dp
247 761250 : shell_particle_set(i)%f(3) = 0.0_dp
248 761250 : core_particle_set(i)%f(1) = 0.0_dp
249 761250 : core_particle_set(i)%f(2) = 0.0_dp
250 770682 : core_particle_set(i)%f(3) = 0.0_dp
251 : END DO
252 : END IF
253 :
254 80943 : IF (use_virial) THEN
255 17120 : pv_bc = 0.0_dp
256 17120 : pv_bond = 0.0_dp
257 17120 : pv_bend = 0.0_dp
258 17120 : pv_torsion = 0.0_dp
259 17120 : pv_imptors = 0.0_dp
260 17120 : pv_opbend = 0.0_dp
261 17120 : pv_urey_bradley = 0.0_dp
262 17120 : pv_nonbond = 0.0_dp
263 17120 : pv_g = 0.0_dp
264 222560 : virial%pv_virial = 0.0_dp
265 : END IF
266 :
267 80943 : pot_nonbond = 0.0_dp
268 80943 : pot_manybody = 0.0_dp
269 80943 : pot_bond = 0.0_dp
270 80943 : pot_bend = 0.0_dp
271 80943 : pot_torsion = 0.0_dp
272 80943 : pot_opbend = 0.0_dp
273 80943 : pot_imptors = 0.0_dp
274 80943 : pot_urey_bradley = 0.0_dp
275 80943 : pot_shell = 0.0_dp
276 80943 : vg_coulomb = 0.0_dp
277 80943 : thermo%pot = 0.0_dp
278 80943 : thermo%harm_shell = 0.0_dp
279 :
280 : ! Get real-space non-bonded forces:
281 80943 : IF (iw > 0) THEN
282 285 : WRITE (iw, '(A)') " FIST:: FORCES IN INPUT..."
283 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
284 : END IF
285 :
286 80943 : IF (fist_nonbond_env%do_nonbonded) THEN
287 : ! Compute density for EAM
288 80791 : CALL density_nonbond(fist_nonbond_env, particle_set, cell, para_env)
289 :
290 : ! Compute embedding function and manybody energy
291 : CALL energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, &
292 80791 : cell, pot_manybody, para_env, mm_section, use_virial)
293 :
294 : ! Nonbond contribution + Manybody Forces
295 80791 : IF (shell_present) THEN
296 : CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
297 : pot_nonbond, f_nonbond, pv_nonbond, &
298 : fshell_nonbond=fshell_nonbond, fcore_nonbond=fcore_nonbond, &
299 : atprop_env=atprop_env, &
300 9344 : atomic_kind_set=atomic_kind_set, use_virial=use_virial)
301 : ELSE
302 : CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
303 : pot_nonbond, f_nonbond, pv_nonbond, atprop_env=atprop_env, &
304 71447 : atomic_kind_set=atomic_kind_set, use_virial=use_virial)
305 : CALL force_nonbond_manybody(fist_nonbond_env, particle_set, cell, f_nonbond, pv_nonbond, &
306 71447 : use_virial=use_virial)
307 : END IF
308 : END IF
309 :
310 80943 : IF (iw > 0) THEN
311 285 : WRITE (iw, '(A)') " FIST:: NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
312 285 : WRITE (iw, '(3f15.9)') f_nonbond
313 285 : IF (shell_present .AND. shell_model_ad) THEN
314 44 : WRITE (iw, '(A)') " FIST:: SHELL NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
315 44 : WRITE (iw, '(3f15.9)') fshell_nonbond
316 44 : WRITE (iw, '(A)') " FIST:: CORE NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
317 44 : WRITE (iw, '(3f15.9)') fcore_nonbond
318 : END IF
319 : END IF
320 :
321 : ! Get g-space non-bonded forces:
322 80943 : IF (ewald_type /= do_ewald_none) THEN
323 : ! Determine size of the forces array
324 : SELECT CASE (ewald_type)
325 : CASE (do_ewald_ewald)
326 31318 : fg_coulomb_size = nlocal_particles
327 : CASE DEFAULT
328 61035 : fg_coulomb_size = natoms
329 : END SELECT
330 : ! Allocate and zeroing arrays
331 182559 : ALLOCATE (fg_coulomb(3, fg_coulomb_size))
332 27634675 : fg_coulomb = 0.0_dp
333 61035 : IF (shell_present) THEN
334 28254 : ALLOCATE (fgshell_coulomb(3, nshell))
335 28254 : ALLOCATE (fgcore_coulomb(3, nshell))
336 3054338 : fgshell_coulomb = 0.0_dp
337 3054338 : fgcore_coulomb = 0.0_dp
338 : END IF
339 61035 : IF (shell_present .AND. do_multipoles) THEN
340 0 : CPABORT("Multipoles and Core-Shell model not implemented.")
341 : END IF
342 : ! If not multipole: Compute self-interaction and neutralizing background
343 : ! for multipoles is handled separately..
344 61035 : IF (.NOT. do_multipoles) THEN
345 : CALL ewald_self(ewald_env, cell, atomic_kind_set, local_particles, &
346 59877 : thermo%e_self, thermo%e_neut, fist_nonbond_env%charges)
347 59877 : IF (atprop_env%energy) THEN
348 : CALL ewald_self_atom(ewald_env, atomic_kind_set, local_particles, &
349 636 : atprop_env%atener, fist_nonbond_env%charges)
350 6054 : atprop_env%atener = atprop_env%atener + thermo%e_neut/SIZE(atprop_env%atener)
351 : END IF
352 : END IF
353 :
354 : ! Polarizable force-field
355 61035 : IF (do_ipol /= do_fist_pol_none) THEN
356 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
357 104 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
358 0 : CPABORT("Polarizable force field and array charges not implemented!")
359 : END IF
360 : ! Converge the dipoles self-consistently
361 : CALL fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
362 : ewald_pw, fist_nonbond_env, cell, particle_set, &
363 : local_particles, thermo, vg_coulomb, pot_nonbond, f_nonbond, &
364 104 : fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section, do_ipol)
365 : ELSE
366 : ! Non-Polarizable force-field
367 29613 : SELECT CASE (ewald_type)
368 : CASE (do_ewald_ewald)
369 : ! Parallelisation over atoms --> allocate local atoms
370 29613 : IF (shell_present) THEN
371 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
372 0 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
373 0 : CPABORT("Core-Shell and array charges not implemented!")
374 : END IF
375 0 : IF (do_multipoles) THEN
376 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within Ewald sum!")
377 : ELSE
378 0 : CPABORT("Core-Shell model not yet implemented within Ewald sums.")
379 : END IF
380 : ELSE
381 29613 : IF (do_multipoles) THEN
382 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
383 1054 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
384 0 : CPABORT("Multipole Ewald and array charges not implemented!")
385 : END IF
386 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
387 : particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
388 : thermo%e_self, multipoles%task, do_correction_bonded=.TRUE., do_forces=.TRUE., &
389 : do_stress=use_virial, do_efield=.FALSE., radii=multipoles%radii, &
390 : charges=multipoles%charges, dipoles=multipoles%dipoles, &
391 : quadrupoles=multipoles%quadrupoles, forces_local=fg_coulomb, &
392 : forces_glob=f_nonbond, pv_local=pv_g, pv_glob=pv_nonbond, iw=iw2, &
393 1054 : do_debug=.TRUE., atomic_kind_set=atomic_kind_set, mm_section=mm_section)
394 : ELSE
395 28559 : IF (atprop_env%energy) THEN
396 505 : ALLOCATE (e_coulomb(fg_coulomb_size))
397 : END IF
398 : CALL ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
399 : local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial=use_virial, &
400 28559 : charges=fist_nonbond_env%charges, e_coulomb=e_coulomb)
401 : END IF
402 : END IF
403 : CASE (do_ewald_pme)
404 : ! Parallelisation over grids --> allocate all atoms
405 1818 : IF (shell_present) THEN
406 : ! Check if an array of chagres was provided and in case abort due to lack of implementation
407 0 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
408 0 : CPABORT("Core-Shell and array charges not implemented!")
409 : END IF
410 0 : IF (do_multipoles) THEN
411 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a PME scheme!")
412 : ELSE
413 : CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
414 : pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
415 : fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
416 0 : atprop=atprop_env)
417 0 : CALL para_env%sum(fgshell_coulomb)
418 0 : CALL para_env%sum(fgcore_coulomb)
419 : END IF
420 : ELSE
421 1818 : IF (do_multipoles) THEN
422 0 : CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
423 : ELSE
424 : CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
425 : pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
426 1818 : atprop=atprop_env)
427 : END IF
428 : END IF
429 1818 : CALL para_env%sum(fg_coulomb)
430 : CASE (do_ewald_spme)
431 : ! Parallelisation over grids --> allocate all atoms
432 29500 : IF (shell_present) THEN
433 : ! Check if an array of charges was provided and in case abort due to lack of implementation
434 9418 : IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
435 0 : CPABORT("Core-Shell and array charges not implemented!")
436 : END IF
437 9418 : IF (do_multipoles) THEN
438 0 : CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a SPME scheme!")
439 : ELSE
440 : CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
441 : pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
442 : fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
443 9418 : atprop=atprop_env)
444 9418 : CALL para_env%sum(fgshell_coulomb)
445 9418 : CALL para_env%sum(fgcore_coulomb)
446 : END IF
447 : ELSE
448 20082 : IF (do_multipoles) THEN
449 0 : CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
450 : ELSE
451 : CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
452 : pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
453 20082 : atprop=atprop_env)
454 : END IF
455 : END IF
456 90431 : CALL para_env%sum(fg_coulomb)
457 : END SELECT
458 : END IF
459 :
460 : ! Subtract the interaction between screening charges. This is a
461 : ! correction in real-space and depends on the neighborlists. Therefore
462 : ! it is only executed if fist_nonbond_env%do_nonbonded is set.
463 61035 : IF (fist_nonbond_env%do_nonbonded) THEN
464 : ! Correction for core-shell model
465 60947 : IF (shell_present) THEN
466 : CALL bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
467 : local_particles, particle_set, ewald_env, thermo%e_bonded, &
468 : pv_bc, shell_particle_set=shell_particle_set, &
469 : core_particle_set=core_particle_set, atprop_env=atprop_env, cell=cell, &
470 9330 : use_virial=use_virial)
471 : ELSE
472 51617 : IF (.NOT. do_multipoles) THEN
473 : CALL bonded_correct_gaussian(fist_nonbond_env, &
474 : atomic_kind_set, local_particles, particle_set, &
475 : ewald_env, thermo%e_bonded, pv_bc=pv_bc, atprop_env=atprop_env, cell=cell, &
476 50459 : use_virial=use_virial)
477 : END IF
478 : END IF
479 : END IF
480 :
481 61035 : IF (.NOT. do_multipoles) THEN
482 : ! Multipole code has its own printing routine.
483 : CALL ewald_print(iw2, pot_nonbond, vg_coulomb, thermo%e_self, thermo%e_neut, &
484 59877 : thermo%e_bonded)
485 : END IF
486 : ELSE
487 19908 : IF (use_virial) THEN
488 3180 : pv_g = 0.0_dp
489 3180 : pv_bc = 0.0_dp
490 : END IF
491 19908 : thermo%e_neut = 0.0_dp
492 : END IF
493 :
494 80943 : IF (iw > 0) THEN
495 285 : IF (ALLOCATED(fg_coulomb)) THEN
496 206 : WRITE (iw, '(A)') " FIST:: NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
497 206 : WRITE (iw, '(3f15.9)') ((fg_coulomb(j, i), j=1, 3), i=1, SIZE(fg_coulomb, 2))
498 : END IF
499 285 : IF (shell_present .AND. shell_model_ad .AND. ALLOCATED(fgshell_coulomb)) THEN
500 44 : WRITE (iw, '(A)') " FIST:: SHELL NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
501 44 : WRITE (iw, '(3f15.9)') ((fgshell_coulomb(j, i), j=1, 3), i=1, SIZE(fgshell_coulomb, 2))
502 44 : WRITE (iw, '(A)') " FIST:: CORE NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
503 44 : WRITE (iw, '(3f15.9)') ((fgcore_coulomb(j, i), j=1, 3), i=1, SIZE(fgcore_coulomb, 2))
504 : END IF
505 : END IF
506 :
507 : ! calculate the action of an (external) electric field
508 80943 : IF (efield%apply_field) THEN
509 348 : ALLOCATE (ef_f(3, natoms))
510 : CALL fist_efield_energy_force(ef_ener, ef_f, ef_pv, atomic_kind_set, particle_set, cell, &
511 116 : efield, use_virial=use_virial, iunit=iw, charges=fist_nonbond_env%charges)
512 : END IF
513 :
514 : ! Get intramolecular forces
515 80943 : IF (PRESENT(debug)) THEN
516 : CALL force_intra_control(molecule_set, molecule_kind_set, &
517 : local_molecules, particle_set, shell_particle_set, &
518 : core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
519 : pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
520 : pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
521 : debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, &
522 0 : debug%f_imptors, debug%f_opbend, cell, use_virial, atprop_env)
523 :
524 : ELSE
525 : CALL force_intra_control(molecule_set, molecule_kind_set, &
526 : local_molecules, particle_set, shell_particle_set, &
527 : core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
528 : pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
529 : pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
530 80943 : cell=cell, use_virial=use_virial, atprop_env=atprop_env)
531 : END IF
532 :
533 : ! Perform global sum of the intra-molecular (bonded) forces for the core-shell atoms
534 80943 : IF (shell_present .AND. shell_model_ad) THEN
535 28296 : ALLOCATE (fcore_shell_bonded(3, nshell))
536 3054432 : fcore_shell_bonded = 0.0_dp
537 805170 : DO i = 1, natoms
538 795738 : shell_index = particle_set(i)%shell_index
539 805170 : IF (shell_index /= 0) THEN
540 3045000 : fcore_shell_bonded(1:3, shell_index) = particle_set(i)%f(1:3)
541 : END IF
542 : END DO
543 9432 : CALL para_env%sum(fcore_shell_bonded)
544 805170 : DO i = 1, natoms
545 795738 : shell_index = particle_set(i)%shell_index
546 805170 : IF (shell_index /= 0) THEN
547 3045000 : particle_set(i)%f(1:3) = fcore_shell_bonded(1:3, shell_index)
548 : END IF
549 : END DO
550 9432 : DEALLOCATE (fcore_shell_bonded)
551 : END IF
552 :
553 80943 : IF (iw > 0) THEN
554 285 : xdum1 = cp_unit_from_cp2k(pot_bond, "kcalmol")
555 285 : xdum2 = cp_unit_from_cp2k(pot_bend, "kcalmol")
556 285 : xdum3 = cp_unit_from_cp2k(pot_urey_bradley, "kcalmol")
557 285 : WRITE (iw, '(A)') " FIST energy contributions in kcal/mol:"
558 : WRITE (iw, '(1x,"BOND = ",f13.4,'// &
559 : '2x,"ANGLE = ",f13.4,'// &
560 285 : '2x,"UBRAD = ",f13.4)') xdum1, xdum2, xdum3
561 285 : xdum1 = cp_unit_from_cp2k(pot_torsion, "kcalmol")
562 285 : xdum2 = cp_unit_from_cp2k(pot_imptors, "kcalmol")
563 285 : xdum3 = cp_unit_from_cp2k(pot_opbend, "kcalmol")
564 : WRITE (iw, '(1x,"TORSION = ",f13.4,'// &
565 : '2x,"IMPTORS = ",f13.4,'// &
566 285 : '2x,"OPBEND = ",f13.4)') xdum1, xdum2, xdum3
567 :
568 285 : WRITE (iw, '(A)') " FIST:: CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
569 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
570 285 : IF (shell_present .AND. shell_model_ad) THEN
571 44 : WRITE (iw, '(A)') " FIST:: SHELL CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
572 16940 : WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, SIZE(shell_particle_set))
573 44 : WRITE (iw, '(A)') " FIST:: CORE CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
574 16940 : WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, SIZE(core_particle_set))
575 : END IF
576 : END IF
577 :
578 : ! add up all the potential energies
579 : thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_opbend + &
580 80943 : pot_imptors + pot_urey_bradley + pot_manybody + pot_shell
581 :
582 80943 : CALL para_env%sum(thermo%pot)
583 :
584 80943 : IF (shell_present) THEN
585 9432 : thermo%harm_shell = pot_shell
586 9432 : CALL para_env%sum(thermo%harm_shell)
587 : END IF
588 : ! add g-space contributions if needed
589 80943 : IF (ewald_type /= do_ewald_none) THEN
590 : ! e_self, e_neut, and ebonded are already summed over all processors
591 : ! vg_coulomb is not calculated in parallel
592 61035 : thermo%e_gspace = vg_coulomb
593 61035 : thermo%pot = thermo%pot + thermo%e_self + thermo%e_neut
594 61035 : thermo%pot = thermo%pot + vg_coulomb + thermo%e_bonded
595 : ! add the induction energy of the dipoles for polarizable models
596 61035 : IF (do_ipol /= do_fist_pol_none) thermo%pot = thermo%pot + thermo%e_induction
597 : END IF
598 :
599 : ! add up all the forces
600 : !
601 : ! nonbonded forces might be calculated for atoms not on this node
602 : ! ewald forces are strictly local -> sum only over pnode
603 : ! We first sum the forces in f_nonbond, this allows for a more efficient
604 : ! global sum in the parallel code and in the end copy them back to part
605 242829 : ALLOCATE (f_total(3, natoms))
606 32985667 : f_total = 0.0_dp
607 8307124 : DO i = 1, natoms
608 8226181 : f_total(1, i) = particle_set(i)%f(1) + f_nonbond(1, i)
609 8226181 : f_total(2, i) = particle_set(i)%f(2) + f_nonbond(2, i)
610 8307124 : f_total(3, i) = particle_set(i)%f(3) + f_nonbond(3, i)
611 : END DO
612 80943 : IF (shell_present) THEN
613 28296 : ALLOCATE (fshell_total(3, nshell))
614 3054432 : fshell_total = 0.0_dp
615 28296 : ALLOCATE (fcore_total(3, nshell))
616 3054432 : fcore_total = 0.0_dp
617 770682 : DO i = 1, nshell
618 761250 : fcore_total(1, i) = core_particle_set(i)%f(1) + fcore_nonbond(1, i)
619 761250 : fcore_total(2, i) = core_particle_set(i)%f(2) + fcore_nonbond(2, i)
620 761250 : fcore_total(3, i) = core_particle_set(i)%f(3) + fcore_nonbond(3, i)
621 761250 : fshell_total(1, i) = shell_particle_set(i)%f(1) + fshell_nonbond(1, i)
622 761250 : fshell_total(2, i) = shell_particle_set(i)%f(2) + fshell_nonbond(2, i)
623 770682 : fshell_total(3, i) = shell_particle_set(i)%f(3) + fshell_nonbond(3, i)
624 : END DO
625 : END IF
626 :
627 80943 : IF (iw > 0) THEN
628 285 : WRITE (iw, '(A)') " FIST::(1)INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
629 285 : WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
630 285 : IF (shell_present .AND. shell_model_ad) THEN
631 44 : WRITE (iw, '(A)') " FIST::(1)SHELL INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
632 44 : WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
633 44 : WRITE (iw, '(A)') " FIST::(1)CORE INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
634 44 : WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
635 : END IF
636 : END IF
637 :
638 : ! Adding in the reciprocal forces: EWALD is a special case because of distrubted data
639 80943 : IF (ewald_type == do_ewald_ewald) THEN
640 : node = 0
641 131093 : DO iparticle_kind = 1, nparticle_kind
642 101376 : nparticle_local = local_particles%n_el(iparticle_kind)
643 559433 : DO iparticle_local = 1, nparticle_local
644 428340 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
645 428340 : node = node + 1
646 428340 : f_total(1, ii) = f_total(1, ii) + fg_coulomb(1, node)
647 428340 : f_total(2, ii) = f_total(2, ii) + fg_coulomb(2, node)
648 428340 : f_total(3, ii) = f_total(3, ii) + fg_coulomb(3, node)
649 428340 : IF (PRESENT(debug)) THEN
650 0 : debug%f_g(1, ii) = debug%f_g(1, ii) + fg_coulomb(1, node)
651 0 : debug%f_g(2, ii) = debug%f_g(2, ii) + fg_coulomb(2, node)
652 0 : debug%f_g(3, ii) = debug%f_g(3, ii) + fg_coulomb(3, node)
653 : END IF
654 529716 : IF (atprop_env%energy) THEN
655 303 : atprop_env%atener(ii) = atprop_env%atener(ii) + e_coulomb(node)
656 : END IF
657 : END DO
658 : END DO
659 29717 : IF (atprop_env%energy) THEN
660 202 : DEALLOCATE (e_coulomb)
661 : END IF
662 : END IF
663 :
664 80943 : IF (iw > 0) THEN
665 285 : WRITE (iw, '(A)') " FIST::(2)TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
666 285 : WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
667 285 : IF (shell_present .AND. shell_model_ad) THEN
668 44 : WRITE (iw, '(A)') " FIST::(2)SHELL TOTAL FORCES(1)+ ELECTROSTATIC FORCES "
669 44 : WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
670 44 : WRITE (iw, '(A)') " FIST::(2)CORE TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
671 44 : WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
672 : END IF
673 : END IF
674 :
675 80943 : IF (use_virial) THEN
676 : ! Add up all the pressure tensors
677 17120 : IF (ewald_type == do_ewald_none) THEN
678 41340 : virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley
679 3180 : CALL para_env%sum(virial%pv_virial)
680 : ELSE
681 13940 : ident = 0.0_dp
682 55760 : DO i = 1, 3
683 55760 : ident(i, i) = 1.0_dp
684 : END DO
685 :
686 181220 : virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc
687 13940 : CALL para_env%sum(virial%pv_virial)
688 :
689 181220 : virial%pv_virial = virial%pv_virial + ident*thermo%e_neut
690 181220 : virial%pv_virial = virial%pv_virial + pv_g
691 : END IF
692 : END IF
693 :
694 : ! Sum total forces
695 80943 : CALL para_env%sum(f_total)
696 80943 : IF (shell_present .AND. shell_model_ad) THEN
697 9432 : CALL para_env%sum(fcore_total)
698 9432 : CALL para_env%sum(fshell_total)
699 : END IF
700 :
701 : ! contributions from fields (currently all quantities are fully replicated!)
702 80943 : IF (efield%apply_field) THEN
703 116 : thermo%pot = thermo%pot + ef_ener
704 1508 : f_total(1:3, 1:natoms) = f_total(1:3, 1:natoms) + ef_f(1:3, 1:natoms)
705 116 : IF (use_virial) THEN
706 26 : virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + ef_pv(1:3, 1:3)
707 : END IF
708 116 : DEALLOCATE (ef_f)
709 : END IF
710 :
711 : ! Assign to particle_set
712 31318 : SELECT CASE (ewald_type)
713 : CASE (do_ewald_spme, do_ewald_pme)
714 31318 : IF (shell_present .AND. shell_model_ad) THEN
715 805128 : DO i = 1, natoms
716 795710 : shell_index = particle_set(i)%shell_index
717 805128 : IF (shell_index == 0) THEN
718 137920 : particle_set(i)%f(1:3) = f_total(1:3, i) + fg_coulomb(1:3, i)
719 : ELSE
720 761230 : atomic_kind => particle_set(i)%atomic_kind
721 761230 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
722 761230 : fc = shell%mass_core/mass
723 3044920 : fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
724 761230 : fs = shell%mass_shell/mass
725 3044920 : fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
726 : END IF
727 : END DO
728 :
729 770648 : DO i = 1, nshell
730 761230 : core_particle_set(i)%f(1) = fcore_total(1, i) + fgcore_coulomb(1, i)
731 761230 : core_particle_set(i)%f(2) = fcore_total(2, i) + fgcore_coulomb(2, i)
732 761230 : core_particle_set(i)%f(3) = fcore_total(3, i) + fgcore_coulomb(3, i)
733 761230 : shell_particle_set(i)%f(1) = fshell_total(1, i) + fgshell_coulomb(1, i)
734 761230 : shell_particle_set(i)%f(2) = fshell_total(2, i) + fgshell_coulomb(2, i)
735 770648 : shell_particle_set(i)%f(3) = fshell_total(3, i) + fgshell_coulomb(3, i)
736 : END DO
737 :
738 21900 : ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
739 0 : CPABORT("Non adiabatic shell-model not implemented.")
740 : ELSE
741 5691260 : DO i = 1, natoms
742 5669360 : particle_set(i)%f(1) = f_total(1, i) + fg_coulomb(1, i)
743 5669360 : particle_set(i)%f(2) = f_total(2, i) + fg_coulomb(2, i)
744 5691260 : particle_set(i)%f(3) = f_total(3, i) + fg_coulomb(3, i)
745 : END DO
746 : END IF
747 : CASE (do_ewald_ewald, do_ewald_none)
748 80943 : IF (shell_present .AND. shell_model_ad) THEN
749 42 : DO i = 1, natoms
750 28 : shell_index = particle_set(i)%shell_index
751 42 : IF (shell_index == 0) THEN
752 32 : particle_set(i)%f(1:3) = f_total(1:3, i)
753 : ELSE
754 20 : atomic_kind => particle_set(i)%atomic_kind
755 20 : CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
756 20 : fc = shell%mass_core/mass
757 80 : fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
758 20 : fs = shell%mass_shell/mass
759 80 : fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
760 : END IF
761 : END DO
762 34 : DO i = 1, nshell
763 20 : core_particle_set(i)%f(1) = fcore_total(1, i)
764 20 : core_particle_set(i)%f(2) = fcore_total(2, i)
765 20 : core_particle_set(i)%f(3) = fcore_total(3, i)
766 20 : shell_particle_set(i)%f(1) = fshell_total(1, i)
767 20 : shell_particle_set(i)%f(2) = fshell_total(2, i)
768 34 : shell_particle_set(i)%f(3) = fshell_total(3, i)
769 : END DO
770 49611 : ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
771 0 : CPABORT("Non adiabatic shell-model not implemented.")
772 : ELSE
773 1810694 : DO i = 1, natoms
774 1761083 : particle_set(i)%f(1) = f_total(1, i)
775 1761083 : particle_set(i)%f(2) = f_total(2, i)
776 1810694 : particle_set(i)%f(3) = f_total(3, i)
777 : END DO
778 : END IF
779 : END SELECT
780 :
781 80943 : IF (iw > 0) THEN
782 285 : WRITE (iw, '(A)') " FIST::(3)TOTAL FORCES - THE END..."
783 48061 : WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, natoms)
784 285 : IF (shell_present .AND. shell_model_ad) THEN
785 44 : WRITE (iw, '(A)') " FIST::(3)SHELL TOTAL FORCES - THE END..."
786 16940 : WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, nshell)
787 44 : WRITE (iw, '(A)') " FIST::(3)CORE TOTAL FORCES - THE END..."
788 16940 : WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, nshell)
789 : END IF
790 285 : WRITE (iw, '(A,f15.9)') "Energy after FIST calculation.. exiting now ::", thermo%pot
791 : END IF
792 : !
793 : ! If we are doing debugging, check if variables are present and assign
794 : !
795 80943 : IF (PRESENT(debug)) THEN
796 0 : CALL para_env%sum(pot_manybody)
797 0 : debug%pot_manybody = pot_manybody
798 0 : CALL para_env%sum(pot_nonbond)
799 0 : debug%pot_nonbond = pot_nonbond
800 0 : CALL para_env%sum(pot_bond)
801 0 : debug%pot_bond = pot_bond
802 0 : CALL para_env%sum(pot_bend)
803 0 : debug%pot_bend = pot_bend
804 0 : CALL para_env%sum(pot_torsion)
805 0 : debug%pot_torsion = pot_torsion
806 0 : CALL para_env%sum(pot_imptors)
807 0 : debug%pot_imptors = pot_imptors
808 0 : CALL para_env%sum(pot_opbend)
809 0 : debug%pot_opbend = pot_opbend
810 0 : CALL para_env%sum(pot_urey_bradley)
811 0 : debug%pot_urey_bradley = pot_urey_bradley
812 0 : CALL para_env%sum(f_nonbond)
813 0 : debug%f_nonbond = f_nonbond
814 0 : IF (use_virial) THEN
815 0 : CALL para_env%sum(pv_nonbond)
816 0 : debug%pv_nonbond = pv_nonbond
817 0 : CALL para_env%sum(pv_bond)
818 0 : debug%pv_bond = pv_bond
819 0 : CALL para_env%sum(pv_bend)
820 0 : debug%pv_bend = pv_bend
821 0 : CALL para_env%sum(pv_torsion)
822 0 : debug%pv_torsion = pv_torsion
823 0 : CALL para_env%sum(pv_imptors)
824 0 : debug%pv_imptors = pv_imptors
825 0 : CALL para_env%sum(pv_urey_bradley)
826 0 : debug%pv_ub = pv_urey_bradley
827 : END IF
828 0 : SELECT CASE (ewald_type)
829 : CASE (do_ewald_spme, do_ewald_pme)
830 0 : debug%pot_g = vg_coulomb
831 0 : debug%pv_g = pv_g
832 0 : debug%f_g = fg_coulomb
833 : CASE (do_ewald_ewald)
834 0 : debug%pot_g = vg_coulomb
835 0 : IF (use_virial) debug%pv_g = pv_g
836 : CASE default
837 0 : debug%pot_g = 0.0_dp
838 0 : debug%f_g = 0.0_dp
839 0 : IF (use_virial) debug%pv_g = 0.0_dp
840 : END SELECT
841 : END IF
842 :
843 : ! print properties if requested
844 80943 : print_section => section_vals_get_subs_vals(mm_section, "PRINT")
845 80943 : CALL print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
846 :
847 : ! deallocating all local variables
848 80943 : IF (ALLOCATED(fg_coulomb)) THEN
849 61035 : DEALLOCATE (fg_coulomb)
850 : END IF
851 80943 : IF (ALLOCATED(f_total)) THEN
852 80943 : DEALLOCATE (f_total)
853 : END IF
854 80943 : DEALLOCATE (f_nonbond)
855 80943 : IF (shell_present) THEN
856 9432 : DEALLOCATE (fshell_total)
857 9432 : IF (ewald_type /= do_ewald_none) THEN
858 9418 : DEALLOCATE (fgshell_coulomb)
859 : END IF
860 9432 : DEALLOCATE (fshell_nonbond)
861 : END IF
862 80943 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%DERIVATIVES")
863 80943 : CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%EWALD_INFO")
864 80943 : CALL timestop(handle)
865 :
866 242829 : END SUBROUTINE fist_calc_energy_force
867 :
868 : ! **************************************************************************************************
869 : !> \brief Print properties number according the requests in input file
870 : !> \param fist_env ...
871 : !> \param print_section ...
872 : !> \param atomic_kind_set ...
873 : !> \param particle_set ...
874 : !> \param cell ...
875 : !> \par History
876 : !> [01.2006] created
877 : !> \author Teodoro Laino
878 : ! **************************************************************************************************
879 80943 : SUBROUTINE print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
880 : TYPE(fist_environment_type), POINTER :: fist_env
881 : TYPE(section_vals_type), POINTER :: print_section
882 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
883 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
884 : TYPE(cell_type), POINTER :: cell
885 :
886 : INTEGER :: unit_nr
887 : TYPE(cp_logger_type), POINTER :: logger
888 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
889 : TYPE(section_vals_type), POINTER :: print_key
890 :
891 80943 : NULLIFY (logger, print_key, fist_nonbond_env)
892 80943 : logger => cp_get_default_logger()
893 80943 : print_key => section_vals_get_subs_vals(print_section, "dipole")
894 80943 : CALL fist_env_get(fist_env, fist_nonbond_env=fist_nonbond_env)
895 80943 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
896 : cp_p_file)) THEN
897 : unit_nr = cp_print_key_unit_nr(logger, print_section, "dipole", &
898 21102 : extension=".data", middle_name="MM_DIPOLE", log_filename=.FALSE.)
899 : CALL fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
900 21102 : cell, unit_nr, fist_nonbond_env%charges)
901 21102 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
902 : END IF
903 :
904 80943 : END SUBROUTINE print_fist
905 :
906 0 : END MODULE fist_force
|
