1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Util force_env module 8!> \author Teodoro Laino [tlaino] - 02.2011 9! ************************************************************************************************** 10MODULE force_env_utils 11 12 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 13 USE cell_types, ONLY: cell_type 14 USE constraint, ONLY: rattle_control,& 15 shake_control 16 USE constraint_util, ONLY: getold 17 USE cp_subsys_types, ONLY: cp_subsys_get,& 18 cp_subsys_type 19 USE distribution_1d_types, ONLY: distribution_1d_type 20 USE force_env_types, ONLY: force_env_get,& 21 force_env_type 22 USE input_section_types, ONLY: section_vals_get,& 23 section_vals_get_subs_vals,& 24 section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE mathlib, ONLY: det_3x3,& 28 jacobi 29 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 30 USE molecule_list_types, ONLY: molecule_list_type 31 USE molecule_types, ONLY: global_constraint_type 32 USE particle_list_types, ONLY: particle_list_type 33 USE particle_types, ONLY: update_particle_set 34 USE physcon, ONLY: pascal 35#include "./base/base_uses.f90" 36 37 IMPLICIT NONE 38 39 PRIVATE 40 41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_env_utils' 42 43 PUBLIC :: force_env_shake, & 44 force_env_rattle, & 45 rescale_forces, & 46 write_stress_tensor, & 47 write_forces 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief perform shake (enforcing of constraints) 53!> \param force_env the force env to shake 54!> \param dt the dt for shake (if you are not interested in the velocities 55!> it can be any positive number) 56!> \param shake_tol the tolerance for shake 57!> \param log_unit if >0 then some information on the shake is printed, 58!> defaults to -1 59!> \param lagrange_mult ... 60!> \param dump_lm ... 61!> \param pos ... 62!> \param vel ... 63!> \param compold ... 64!> \param reset ... 65!> \author fawzi 66! ************************************************************************************************** 67 SUBROUTINE force_env_shake(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, & 68 pos, vel, compold, reset) 69 70 TYPE(force_env_type), POINTER :: force_env 71 REAL(kind=dp), INTENT(IN), OPTIONAL :: dt 72 REAL(kind=dp), INTENT(IN) :: shake_tol 73 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult 74 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm 75 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 76 OPTIONAL, TARGET :: pos, vel 77 LOGICAL, INTENT(IN), OPTIONAL :: compold, reset 78 79 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_shake', & 80 routineP = moduleN//':'//routineN 81 82 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, & 83 my_log_unit, nparticle_kind, nparticle_local 84 LOGICAL :: has_pos, has_vel, my_dump_lm 85 REAL(KIND=dp) :: mydt 86 REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_pos, my_vel 87 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 88 TYPE(cell_type), POINTER :: cell 89 TYPE(cp_subsys_type), POINTER :: subsys 90 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 91 TYPE(global_constraint_type), POINTER :: gci 92 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 93 TYPE(molecule_list_type), POINTER :: molecules 94 TYPE(particle_list_type), POINTER :: particles 95 96 CALL timeset(routineN, handle) 97 CPASSERT(ASSOCIATED(force_env)) 98 CPASSERT(force_env%ref_count > 0) 99 my_log_unit = -1 100 IF (PRESENT(log_unit)) my_log_unit = log_unit 101 my_lagrange_mult = -1 102 IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult 103 my_dump_lm = .FALSE. 104 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm 105 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, & 106 my_pos, my_vel, gci) 107 IF (PRESENT(pos)) my_pos => pos 108 IF (PRESENT(vel)) my_vel => vel 109 mydt = 0.1_dp 110 IF (PRESENT(dt)) mydt = dt 111 CALL force_env_get(force_env, subsys=subsys, cell=cell) 112 CALL cp_subsys_get(subsys, & 113 atomic_kinds=atomic_kinds, & 114 local_molecules=local_molecules, & 115 local_particles=local_particles, & 116 molecules=molecules, & 117 molecule_kinds=molecule_kinds, & 118 particles=particles, & 119 gci=gci) 120 nparticle_kind = atomic_kinds%n_els 121 IF (PRESENT(compold)) THEN 122 IF (compold) THEN 123 CALL getold(gci, local_molecules, molecules%els, molecule_kinds%els, & 124 particles%els, cell) 125 END IF 126 END IF 127 has_pos = .FALSE. 128 IF (.NOT. ASSOCIATED(my_pos)) THEN 129 has_pos = .TRUE. 130 ALLOCATE (my_pos(3, particles%n_els)) 131 my_pos = 0.0_dp 132 DO iparticle_kind = 1, nparticle_kind 133 nparticle_local = local_particles%n_el(iparticle_kind) 134 DO iparticle_local = 1, nparticle_local 135 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 136 my_pos(:, iparticle) = particles%els(iparticle)%r(:) 137 END DO 138 END DO 139 END IF 140 has_vel = .FALSE. 141 IF (.NOT. ASSOCIATED(my_vel)) THEN 142 has_vel = .TRUE. 143 ALLOCATE (my_vel(3, particles%n_els)) 144 my_vel = 0.0_dp 145 DO iparticle_kind = 1, nparticle_kind 146 nparticle_local = local_particles%n_el(iparticle_kind) 147 DO iparticle_local = 1, nparticle_local 148 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 149 my_vel(:, iparticle) = particles%els(iparticle)%v(:) 150 END DO 151 END DO 152 END IF 153 154 CALL shake_control(gci=gci, local_molecules=local_molecules, & 155 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, & 156 particle_set=particles%els, pos=my_pos, vel=my_vel, dt=mydt, & 157 shake_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, & 158 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env%group, & 159 local_particles=local_particles) 160 161 ! Possibly reset the lagrange multipliers 162 IF (PRESENT(reset)) THEN 163 IF (reset) THEN 164 ! Reset Intramolecular constraints 165 DO i = 1, SIZE(molecules%els) 166 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN 167 DO j = 1, SIZE(molecules%els(i)%lci%lcolv) 168 ! Reset langrange multiplier 169 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp 170 END DO 171 END IF 172 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN 173 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3) 174 ! Reset langrange multiplier 175 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp 176 END DO 177 END IF 178 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN 179 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6) 180 ! Reset langrange multiplier 181 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp 182 END DO 183 END IF 184 END DO 185 ! Reset Intermolecular constraints 186 IF (ASSOCIATED(gci)) THEN 187 IF (ASSOCIATED(gci%lcolv)) THEN 188 DO j = 1, SIZE(gci%lcolv) 189 ! Reset langrange multiplier 190 gci%lcolv(j)%lambda = 0.0_dp 191 END DO 192 END IF 193 IF (ASSOCIATED(gci%lg3x3)) THEN 194 DO j = 1, SIZE(gci%lg3x3) 195 ! Reset langrange multiplier 196 gci%lg3x3(j)%lambda = 0.0_dp 197 END DO 198 END IF 199 IF (ASSOCIATED(gci%lg4x6)) THEN 200 DO j = 1, SIZE(gci%lg4x6) 201 ! Reset langrange multiplier 202 gci%lg4x6(j)%lambda = 0.0_dp 203 END DO 204 END IF 205 END IF 206 END IF 207 END IF 208 209 IF (has_pos) THEN 210 CALL update_particle_set(particles%els, force_env%para_env%group, pos=my_pos) 211 DEALLOCATE (my_pos) 212 END IF 213 IF (has_vel) THEN 214 CALL update_particle_set(particles%els, force_env%para_env%group, vel=my_vel) 215 DEALLOCATE (my_vel) 216 END IF 217 CALL timestop(handle) 218 END SUBROUTINE force_env_shake 219 220! ************************************************************************************************** 221!> \brief perform rattle (enforcing of constraints on velocities) 222!> This routine can be easily adapted to performe rattle on whatever 223!> other vector different from forces.. 224!> \param force_env the force env to shake 225!> \param dt the dt for shake (if you are not interested in the velocities 226!> it can be any positive number) 227!> \param shake_tol the tolerance for shake 228!> \param log_unit if >0 then some information on the shake is printed, 229!> defaults to -1 230!> \param lagrange_mult ... 231!> \param dump_lm ... 232!> \param vel ... 233!> \param reset ... 234!> \author tlaino 235! ************************************************************************************************** 236 SUBROUTINE force_env_rattle(force_env, dt, shake_tol, log_unit, lagrange_mult, dump_lm, & 237 vel, reset) 238 239 TYPE(force_env_type), POINTER :: force_env 240 REAL(kind=dp), INTENT(in), OPTIONAL :: dt 241 REAL(kind=dp), INTENT(in) :: shake_tol 242 INTEGER, INTENT(in), OPTIONAL :: log_unit, lagrange_mult 243 LOGICAL, INTENT(IN), OPTIONAL :: dump_lm 244 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 245 OPTIONAL, TARGET :: vel 246 LOGICAL, INTENT(IN), OPTIONAL :: reset 247 248 CHARACTER(len=*), PARAMETER :: routineN = 'force_env_rattle', & 249 routineP = moduleN//':'//routineN 250 251 INTEGER :: handle, i, iparticle, iparticle_kind, iparticle_local, j, my_lagrange_mult, & 252 my_log_unit, nparticle_kind, nparticle_local 253 LOGICAL :: has_vel, my_dump_lm 254 REAL(KIND=dp) :: mydt 255 REAL(KIND=dp), DIMENSION(:, :), POINTER :: my_vel 256 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 257 TYPE(cell_type), POINTER :: cell 258 TYPE(cp_subsys_type), POINTER :: subsys 259 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 260 TYPE(global_constraint_type), POINTER :: gci 261 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 262 TYPE(molecule_list_type), POINTER :: molecules 263 TYPE(particle_list_type), POINTER :: particles 264 265 CALL timeset(routineN, handle) 266 CPASSERT(ASSOCIATED(force_env)) 267 CPASSERT(force_env%ref_count > 0) 268 my_log_unit = -1 269 IF (PRESENT(log_unit)) my_log_unit = log_unit 270 my_lagrange_mult = -1 271 IF (PRESENT(lagrange_mult)) my_lagrange_mult = lagrange_mult 272 my_dump_lm = .FALSE. 273 IF (PRESENT(dump_lm)) my_dump_lm = dump_lm 274 NULLIFY (subsys, cell, molecules, molecule_kinds, local_molecules, particles, & 275 my_vel) 276 IF (PRESENT(vel)) my_vel => vel 277 mydt = 0.1_dp 278 IF (PRESENT(dt)) mydt = dt 279 CALL force_env_get(force_env, subsys=subsys, cell=cell) 280 CALL cp_subsys_get(subsys, & 281 atomic_kinds=atomic_kinds, & 282 local_molecules=local_molecules, & 283 local_particles=local_particles, & 284 molecules=molecules, & 285 molecule_kinds=molecule_kinds, & 286 particles=particles, & 287 gci=gci) 288 nparticle_kind = atomic_kinds%n_els 289 has_vel = .FALSE. 290 IF (.NOT. ASSOCIATED(my_vel)) THEN 291 has_vel = .TRUE. 292 ALLOCATE (my_vel(3, particles%n_els)) 293 my_vel = 0.0_dp 294 DO iparticle_kind = 1, nparticle_kind 295 nparticle_local = local_particles%n_el(iparticle_kind) 296 DO iparticle_local = 1, nparticle_local 297 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 298 my_vel(:, iparticle) = particles%els(iparticle)%v(:) 299 END DO 300 END DO 301 END IF 302 303 CALL rattle_control(gci=gci, local_molecules=local_molecules, & 304 molecule_set=molecules%els, molecule_kind_set=molecule_kinds%els, & 305 particle_set=particles%els, vel=my_vel, dt=mydt, & 306 rattle_tol=shake_tol, log_unit=my_log_unit, lagrange_mult=my_lagrange_mult, & 307 dump_lm=my_dump_lm, cell=cell, group=force_env%para_env%group, & 308 local_particles=local_particles) 309 310 ! Possibly reset the lagrange multipliers 311 IF (PRESENT(reset)) THEN 312 IF (reset) THEN 313 ! Reset Intramolecular constraints 314 DO i = 1, SIZE(molecules%els) 315 IF (ASSOCIATED(molecules%els(i)%lci%lcolv)) THEN 316 DO j = 1, SIZE(molecules%els(i)%lci%lcolv) 317 ! Reset langrange multiplier 318 molecules%els(i)%lci%lcolv(j)%lambda = 0.0_dp 319 END DO 320 END IF 321 IF (ASSOCIATED(molecules%els(i)%lci%lg3x3)) THEN 322 DO j = 1, SIZE(molecules%els(i)%lci%lg3x3) 323 ! Reset langrange multiplier 324 molecules%els(i)%lci%lg3x3(j)%lambda = 0.0_dp 325 END DO 326 END IF 327 IF (ASSOCIATED(molecules%els(i)%lci%lg4x6)) THEN 328 DO j = 1, SIZE(molecules%els(i)%lci%lg4x6) 329 ! Reset langrange multiplier 330 molecules%els(i)%lci%lg4x6(j)%lambda = 0.0_dp 331 END DO 332 END IF 333 END DO 334 ! Reset Intermolecular constraints 335 IF (ASSOCIATED(gci)) THEN 336 IF (ASSOCIATED(gci%lcolv)) THEN 337 DO j = 1, SIZE(gci%lcolv) 338 ! Reset langrange multiplier 339 gci%lcolv(j)%lambda = 0.0_dp 340 END DO 341 END IF 342 IF (ASSOCIATED(gci%lg3x3)) THEN 343 DO j = 1, SIZE(gci%lg3x3) 344 ! Reset langrange multiplier 345 gci%lg3x3(j)%lambda = 0.0_dp 346 END DO 347 END IF 348 IF (ASSOCIATED(gci%lg4x6)) THEN 349 DO j = 1, SIZE(gci%lg4x6) 350 ! Reset langrange multiplier 351 gci%lg4x6(j)%lambda = 0.0_dp 352 END DO 353 END IF 354 END IF 355 END IF 356 END IF 357 358 IF (has_vel) THEN 359 CALL update_particle_set(particles%els, force_env%para_env%group, vel=my_vel) 360 END IF 361 DEALLOCATE (my_vel) 362 CALL timestop(handle) 363 END SUBROUTINE force_env_rattle 364 365! ************************************************************************************************** 366!> \brief Rescale forces if requested 367!> \param force_env the force env to shake 368!> \author tlaino 369! ************************************************************************************************** 370 SUBROUTINE rescale_forces(force_env) 371 TYPE(force_env_type), POINTER :: force_env 372 373 CHARACTER(len=*), PARAMETER :: routineN = 'rescale_forces', routineP = moduleN//':'//routineN 374 375 INTEGER :: handle, iparticle 376 LOGICAL :: explicit 377 REAL(KIND=dp) :: force(3), max_value, mod_force 378 TYPE(cp_subsys_type), POINTER :: subsys 379 TYPE(particle_list_type), POINTER :: particles 380 TYPE(section_vals_type), POINTER :: rescale_force_section 381 382 CALL timeset(routineN, handle) 383 CPASSERT(ASSOCIATED(force_env)) 384 CPASSERT(force_env%ref_count > 0) 385 rescale_force_section => section_vals_get_subs_vals(force_env%force_env_section, "RESCALE_FORCES") 386 CALL section_vals_get(rescale_force_section, explicit=explicit) 387 IF (explicit) THEN 388 CALL section_vals_val_get(rescale_force_section, "MAX_FORCE", r_val=max_value) 389 CALL force_env_get(force_env, subsys=subsys) 390 CALL cp_subsys_get(subsys, particles=particles) 391 DO iparticle = 1, SIZE(particles%els) 392 force = particles%els(iparticle)%f(:) 393 mod_force = SQRT(DOT_PRODUCT(force, force)) 394 IF ((mod_force > max_value) .AND. (mod_force /= 0.0_dp)) THEN 395 force = force/mod_force*max_value 396 particles%els(iparticle)%f(:) = force 397 END IF 398 END DO 399 END IF 400 CALL timestop(handle) 401 END SUBROUTINE rescale_forces 402 403! ************************************************************************************************** 404!> \brief Variable precision output of the stress tensor 405!> 406!> \param pv_virial ... 407!> \param output_unit ... 408!> \param cell ... 409!> \param ndigits ... 410!> \param numerical ... 411!> \author MK (26.08.2010) 412! ************************************************************************************************** 413 SUBROUTINE write_stress_tensor(pv_virial, output_unit, cell, ndigits, numerical) 414 415 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: pv_virial 416 INTEGER, INTENT(IN) :: output_unit 417 TYPE(cell_type), POINTER :: cell 418 INTEGER, INTENT(IN) :: ndigits 419 LOGICAL, INTENT(IN) :: numerical 420 421 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_stress_tensor', & 422 routineP = moduleN//':'//routineN 423 424 CHARACTER(LEN=15) :: fmtstr3 425 CHARACTER(LEN=16) :: fmtstr4 426 CHARACTER(LEN=22) :: fmtstr2 427 CHARACTER(LEN=27) :: fmtstr5 428 CHARACTER(LEN=31) :: fmtstr1 429 INTEGER :: n 430 REAL(KIND=dp), DIMENSION(3) :: eigval 431 REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor 432 433 IF (output_unit > 0) THEN 434 CPASSERT(ASSOCIATED(cell)) 435 stress_tensor(:, :) = pv_virial(:, :)/cell%deth*pascal*1.0E-9_dp 436 n = MIN(MAX(1, ndigits), 20) 437 fmtstr1 = "(/,T2,A,/,/,T13,A1,2( X,A1))" 438 WRITE (UNIT=fmtstr1(22:23), FMT="(I2)") n + 7 439 fmtstr2 = "(T3,A,T5,3(1X,F . ))" 440 WRITE (UNIT=fmtstr2(16:17), FMT="(I2)") n + 7 441 WRITE (UNIT=fmtstr2(19:20), FMT="(I2)") n 442 fmtstr3 = "(/,T3,A,F . )" 443 WRITE (UNIT=fmtstr3(10:11), FMT="(I2)") n + 8 444 WRITE (UNIT=fmtstr3(13:14), FMT="(I2)") n 445 IF (numerical) THEN 446 WRITE (UNIT=output_unit, FMT=fmtstr1) & 447 "NUMERICAL STRESS TENSOR [GPa]", "X", "Y", "Z" 448 ELSE 449 WRITE (UNIT=output_unit, FMT=fmtstr1) & 450 "STRESS TENSOR [GPa]", "X", "Y", "Z" 451 END IF 452 WRITE (UNIT=output_unit, FMT=fmtstr2) "X", stress_tensor(1, 1:3) 453 WRITE (UNIT=output_unit, FMT=fmtstr2) "Y", stress_tensor(2, 1:3) 454 WRITE (UNIT=output_unit, FMT=fmtstr2) "Z", stress_tensor(3, 1:3) 455 fmtstr4 = "(/,T3,A,ES . )" 456 WRITE (UNIT=fmtstr4(11:12), FMT="(I2)") n + 8 457 WRITE (UNIT=fmtstr4(14:15), FMT="(I2)") n 458 WRITE (UNIT=output_unit, FMT=fmtstr4) & 459 "1/3 Trace(stress tensor): ", (stress_tensor(1, 1) + & 460 stress_tensor(2, 2) + & 461 stress_tensor(3, 3))/3.0_dp, & 462 "Det(stress tensor) : ", det_3x3(stress_tensor(:, 1), & 463 stress_tensor(:, 2), & 464 stress_tensor(:, 3)) 465 eigval(:) = 0.0_dp 466 eigvec(:, :) = 0.0_dp 467 CALL jacobi(stress_tensor, eigval, eigvec) 468 fmtstr5 = "(/,/,T2,A,/,/,T5,3F . ,/)" 469 WRITE (UNIT=fmtstr5(20:21), FMT="(I2)") n + 8 470 WRITE (UNIT=fmtstr5(23:24), FMT="(I2)") n 471 WRITE (UNIT=output_unit, FMT=fmtstr5) & 472 "EIGENVECTORS AND EIGENVALUES OF THE STRESS TENSOR", & 473 eigval(1:3) 474 WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(1, 1:3) 475 WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(2, 1:3) 476 WRITE (UNIT=output_unit, FMT=fmtstr2) " ", eigvec(3, 1:3) 477 END IF 478 479 END SUBROUTINE write_stress_tensor 480 481! ************************************************************************************************** 482!> \brief Write forces 483!> 484!> \param particles ... 485!> \param output_unit ... 486!> \param label ... 487!> \param ndigits ... 488!> \param total_force ... 489!> \param grand_total_force ... 490!> \param zero_force_core_shell_atom ... 491!> \author MK (06.09.2010) 492! ************************************************************************************************** 493 SUBROUTINE write_forces(particles, output_unit, label, ndigits, total_force, & 494 grand_total_force, zero_force_core_shell_atom) 495 496 TYPE(particle_list_type), POINTER :: particles 497 INTEGER, INTENT(IN) :: output_unit 498 CHARACTER(LEN=*), INTENT(IN) :: label 499 INTEGER, INTENT(IN) :: ndigits 500 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: total_force 501 REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), & 502 OPTIONAL :: grand_total_force 503 LOGICAL, INTENT(IN), OPTIONAL :: zero_force_core_shell_atom 504 505 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_forces', routineP = moduleN//':'//routineN 506 507 CHARACTER(LEN=23) :: fmtstr3 508 CHARACTER(LEN=36) :: fmtstr2 509 CHARACTER(LEN=46) :: fmtstr1 510 INTEGER :: i, ikind, iparticle, n 511 LOGICAL :: zero_force 512 REAL(KIND=dp), DIMENSION(3) :: f 513 514 IF (output_unit > 0) THEN 515 CPASSERT(ASSOCIATED(particles)) 516 n = MIN(MAX(1, ndigits), 20) 517 fmtstr1 = "(/,T2,A,/,/,T2,A,T11,A,T18,A,T35,A1,2( X,A1))" 518 WRITE (UNIT=fmtstr1(39:40), FMT="(I2)") n + 6 519 fmtstr2 = "(T2,I6,1X,I6,T21,A,T28,3(1X,F . ))" 520 WRITE (UNIT=fmtstr2(33:34), FMT="(I2)") n 521 WRITE (UNIT=fmtstr2(30:31), FMT="(I2)") n + 6 522 fmtstr3 = "(T2,A,T28,4(1X,F . ))" 523 WRITE (UNIT=fmtstr3(20:21), FMT="(I2)") n 524 WRITE (UNIT=fmtstr3(17:18), FMT="(I2)") n + 6 525 IF (PRESENT(zero_force_core_shell_atom)) THEN 526 zero_force = zero_force_core_shell_atom 527 ELSE 528 zero_force = .FALSE. 529 END IF 530 WRITE (UNIT=output_unit, FMT=fmtstr1) & 531 label//" FORCES in [a.u.]", "# Atom", "Kind", "Element", "X", "Y", "Z" 532 total_force(1:3) = 0.0_dp 533 DO iparticle = 1, particles%n_els 534 ikind = particles%els(iparticle)%atomic_kind%kind_number 535 IF (particles%els(iparticle)%atom_index /= 0) THEN 536 i = particles%els(iparticle)%atom_index 537 ELSE 538 i = iparticle 539 END IF 540 IF (zero_force .AND. (particles%els(iparticle)%shell_index /= 0)) THEN 541 f(1:3) = 0.0_dp 542 ELSE 543 f(1:3) = particles%els(iparticle)%f(1:3) 544 END IF 545 WRITE (UNIT=output_unit, FMT=fmtstr2) & 546 i, ikind, particles%els(iparticle)%atomic_kind%element_symbol, f(1:3) 547 total_force(1:3) = total_force(1:3) + f(1:3) 548 END DO 549 WRITE (UNIT=output_unit, FMT=fmtstr3) & 550 "SUM OF "//label//" FORCES", total_force(1:3), SQRT(SUM(total_force(:)**2)) 551 END IF 552 553 IF (PRESENT(grand_total_force)) THEN 554 grand_total_force(1:3) = grand_total_force(1:3) + total_force(1:3) 555 WRITE (UNIT=output_unit, FMT="(A)") "" 556 WRITE (UNIT=output_unit, FMT=fmtstr3) & 557 "GRAND TOTAL FORCE", grand_total_force(1:3), SQRT(SUM(grand_total_force(:)**2)) 558 END IF 559 560 END SUBROUTINE write_forces 561 562END MODULE force_env_utils 563