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