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!> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints 9! ************************************************************************************************** 10MODULE constraint 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind 13 USE cell_types, ONLY: cell_type 14 USE colvar_types, ONLY: colvar_counters 15 USE constraint_3x3, ONLY: rattle_3x3_ext,& 16 rattle_3x3_int,& 17 rattle_roll_3x3_ext,& 18 rattle_roll_3x3_int,& 19 shake_3x3_ext,& 20 shake_3x3_int,& 21 shake_roll_3x3_ext,& 22 shake_roll_3x3_int 23 USE constraint_4x6, ONLY: rattle_4x6_ext,& 24 rattle_4x6_int,& 25 rattle_roll_4x6_ext,& 26 rattle_roll_4x6_int,& 27 shake_4x6_ext,& 28 shake_4x6_int,& 29 shake_roll_4x6_ext,& 30 shake_roll_4x6_int 31 USE constraint_clv, ONLY: & 32 rattle_colv_ext, rattle_colv_int, rattle_roll_colv_ext, rattle_roll_colv_int, & 33 shake_colv_ext, shake_colv_int, shake_roll_colv_ext, shake_roll_colv_int, & 34 shake_update_colv_ext, shake_update_colv_int 35 USE constraint_util, ONLY: check_tol,& 36 get_roll_matrix,& 37 restore_temporary_set,& 38 update_temporary_set 39 USE constraint_vsite, ONLY: shake_vsite_ext,& 40 shake_vsite_int 41 USE cp_log_handling, ONLY: cp_to_string 42 USE cp_para_types, ONLY: cp_para_env_type 43 USE distribution_1d_types, ONLY: distribution_1d_type 44 USE input_constants, ONLY: npt_f_ensemble,& 45 npt_i_ensemble 46 USE input_section_types, ONLY: section_vals_get_subs_vals,& 47 section_vals_type 48 USE kinds, ONLY: default_string_length,& 49 dp 50 USE memory_utilities, ONLY: reallocate 51 USE message_passing, ONLY: mp_sum 52 USE molecule_kind_types, ONLY: get_molecule_kind,& 53 get_molecule_kind_set,& 54 molecule_kind_type 55 USE molecule_types, ONLY: global_constraint_type,& 56 molecule_type 57 USE particle_types, ONLY: particle_type 58 USE simpar_types, ONLY: simpar_type 59#include "./base/base_uses.f90" 60 61 IMPLICIT NONE 62 63 PRIVATE 64 PUBLIC :: shake_control, & 65 rattle_control, & 66 shake_roll_control, & 67 rattle_roll_control, & 68 shake_update_targets 69 70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint' 71 INTEGER, PARAMETER, PRIVATE :: Max_Shake_Iter = 1000 72 73CONTAINS 74 75! ************************************************************************************************** 76!> \brief ... 77!> \param gci ... 78!> \param local_molecules ... 79!> \param molecule_set ... 80!> \param molecule_kind_set ... 81!> \param particle_set ... 82!> \param pos ... 83!> \param vel ... 84!> \param dt ... 85!> \param shake_tol ... 86!> \param log_unit ... 87!> \param lagrange_mult ... 88!> \param dump_lm ... 89!> \param cell ... 90!> \param group ... 91!> \param local_particles ... 92!> \par History 93!> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints 94! ************************************************************************************************** 95 SUBROUTINE shake_control(gci, local_molecules, molecule_set, molecule_kind_set, & 96 particle_set, pos, vel, dt, shake_tol, log_unit, lagrange_mult, dump_lm, & 97 cell, group, local_particles) 98 99 TYPE(global_constraint_type), POINTER :: gci 100 TYPE(distribution_1d_type), POINTER :: local_molecules 101 TYPE(molecule_type), POINTER :: molecule_set(:) 102 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 103 TYPE(particle_type), POINTER :: particle_set(:) 104 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 105 REAL(kind=dp), INTENT(in) :: dt, shake_tol 106 INTEGER, INTENT(in) :: log_unit, lagrange_mult 107 LOGICAL, INTENT(IN) :: dump_lm 108 TYPE(cell_type), POINTER :: cell 109 INTEGER, INTENT(in) :: group 110 TYPE(distribution_1d_type), POINTER :: local_particles 111 112 CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_control', routineP = moduleN//':'//routineN 113 114 INTEGER :: handle, i, ikind, imol, ishake_ext, & 115 ishake_int, k, n3x3con, n4x6con, & 116 nconstraint, nkind, nmol_per_kind, & 117 nvsitecon 118 LOGICAL :: do_ext_constraint 119 REAL(KIND=dp) :: int_max_sigma, mass, max_sigma 120 REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass 121 TYPE(atomic_kind_type), POINTER :: atomic_kind 122 TYPE(colvar_counters) :: ncolv 123 TYPE(molecule_kind_type), POINTER :: molecule_kind 124 TYPE(molecule_type), POINTER :: molecule 125 126 CALL timeset(routineN, handle) 127 nkind = SIZE(molecule_kind_set) 128 DO k = 1, SIZE(pos, 2) 129 atomic_kind => particle_set(k)%atomic_kind 130 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 131 imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp)) 132 END DO 133 do_ext_constraint = (gci%ntot /= 0) 134 ishake_ext = 0 135 max_sigma = -1.0E+10_dp 136 Shake_Inter_Loop: DO WHILE ((ABS(max_sigma) >= shake_tol) .AND. (ishake_ext <= Max_Shake_Iter)) 137 max_sigma = 0.0_dp 138 ishake_ext = ishake_ext + 1 139 ! Intramolecular Constraints 140 MOL: DO ikind = 1, nkind 141 nmol_per_kind = local_molecules%n_el(ikind) 142 DO imol = 1, nmol_per_kind 143 i = local_molecules%list(ikind)%array(imol) 144 molecule => molecule_set(i) 145 molecule_kind => molecule%molecule_kind 146 CALL get_molecule_kind(molecule_kind, ncolv=ncolv, & 147 ng3x3=n3x3con, ng4x6=n4x6con, & 148 nconstraint=nconstraint, nvsite=nvsitecon) 149 IF (nconstraint == 0) CYCLE 150 ishake_int = 0 151 int_max_sigma = -1.0E+10_dp 152 Shake_Intra_Loop: DO WHILE ((ABS(int_max_sigma) >= shake_tol) .AND. (ishake_int <= Max_Shake_Iter)) 153 int_max_sigma = 0.0_dp 154 ishake_int = ishake_int + 1 155 ! 3x3 156 IF (n3x3con /= 0) & 157 CALL shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake_int, & 158 int_max_sigma) 159 ! 4x6 160 IF (n4x6con /= 0) & 161 CALL shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake_int, & 162 int_max_sigma) 163 ! Collective Variables 164 IF (ncolv%ntot /= 0) & 165 CALL shake_colv_int(molecule, particle_set, pos, vel, dt, ishake_int, & 166 cell, imass, int_max_sigma) 167 END DO Shake_Intra_Loop 168 max_sigma = MAX(max_sigma, int_max_sigma) 169 CALL shake_int_info(log_unit, i, ishake_int, max_sigma) 170 ! Virtual Site 171 IF (nvsitecon /= 0) & 172 CALL shake_vsite_int(molecule, pos) 173 END DO 174 END DO MOL 175 ! Intermolecular constraints 176 IF (do_ext_constraint) THEN 177 CALL update_temporary_set(group, pos=pos, vel=vel) 178 ! 3x3 179 IF (gci%ng3x3 /= 0) & 180 CALL shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake_ext, & 181 max_sigma) 182 ! 4x6 183 IF (gci%ng4x6 /= 0) & 184 CALL shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake_ext, & 185 max_sigma) 186 ! Collective Variables 187 IF (gci%ncolv%ntot /= 0) & 188 CALL shake_colv_ext(gci, particle_set, pos, vel, dt, ishake_ext, & 189 cell, imass, max_sigma) 190 ! Virtual Site 191 IF (gci%nvsite /= 0) & 192 CALL shake_vsite_ext(gci, pos) 193 CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel) 194 END IF 195 CALL shake_ext_info(log_unit, ishake_ext, max_sigma) 196 END DO Shake_Inter_Loop 197 CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, & 198 molecule_kind_set, group, "S") 199 200 CALL timestop(handle) 201 END SUBROUTINE shake_control 202 203! ************************************************************************************************** 204!> \brief ... 205!> \param gci ... 206!> \param local_molecules ... 207!> \param molecule_set ... 208!> \param molecule_kind_set ... 209!> \param particle_set ... 210!> \param vel ... 211!> \param dt ... 212!> \param rattle_tol ... 213!> \param log_unit ... 214!> \param lagrange_mult ... 215!> \param dump_lm ... 216!> \param cell ... 217!> \param group ... 218!> \param local_particles ... 219!> \par History 220!> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints 221! ************************************************************************************************** 222 SUBROUTINE rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, & 223 particle_set, vel, dt, rattle_tol, log_unit, lagrange_mult, dump_lm, cell, group, & 224 local_particles) 225 226 TYPE(global_constraint_type), POINTER :: gci 227 TYPE(distribution_1d_type), POINTER :: local_molecules 228 TYPE(molecule_type), POINTER :: molecule_set(:) 229 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 230 TYPE(particle_type), POINTER :: particle_set(:) 231 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 232 REAL(kind=dp), INTENT(in) :: dt, rattle_tol 233 INTEGER, INTENT(in) :: log_unit, lagrange_mult 234 LOGICAL, INTENT(IN) :: dump_lm 235 TYPE(cell_type), POINTER :: cell 236 INTEGER, INTENT(in) :: group 237 TYPE(distribution_1d_type), POINTER :: local_particles 238 239 CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_control', routineP = moduleN//':'//routineN 240 241 INTEGER :: handle, i, ikind, imol, irattle_ext, & 242 irattle_int, k, n3x3con, n4x6con, & 243 nconstraint, nkind, nmol_per_kind 244 LOGICAL :: do_ext_constraint 245 REAL(KIND=dp) :: int_max_sigma, mass, max_sigma 246 REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass 247 TYPE(atomic_kind_type), POINTER :: atomic_kind 248 TYPE(colvar_counters) :: ncolv 249 TYPE(molecule_kind_type), POINTER :: molecule_kind 250 TYPE(molecule_type), POINTER :: molecule 251 252 CALL timeset(routineN, handle) 253 nkind = SIZE(molecule_kind_set) 254 DO k = 1, SIZE(vel, 2) 255 atomic_kind => particle_set(k)%atomic_kind 256 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 257 imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp)) 258 END DO 259 do_ext_constraint = (gci%ntot /= 0) 260 irattle_ext = 0 261 max_sigma = -1.0E+10_dp 262 Rattle_Inter_Loop: DO WHILE (ABS(max_sigma) >= rattle_tol) 263 max_sigma = 0.0_dp 264 irattle_ext = irattle_ext + 1 265 ! Intramolecular Constraints 266 MOL: DO ikind = 1, nkind 267 nmol_per_kind = local_molecules%n_el(ikind) 268 DO imol = 1, nmol_per_kind 269 i = local_molecules%list(ikind)%array(imol) 270 molecule => molecule_set(i) 271 molecule_kind => molecule%molecule_kind 272 CALL get_molecule_kind(molecule_kind, ncolv=ncolv, ng3x3=n3x3con, & 273 ng4x6=n4x6con, nconstraint=nconstraint) 274 IF (nconstraint == 0) CYCLE 275 irattle_int = 0 276 int_max_sigma = -1.0E+10_dp 277 Rattle_Intra_Loop: DO WHILE (ABS(int_max_sigma) >= rattle_tol) 278 int_max_sigma = 0.0_dp 279 irattle_int = irattle_int + 1 280 ! 3x3 281 IF (n3x3con /= 0) & 282 CALL rattle_3x3_int(molecule, particle_set, vel, dt) 283 ! 4x6 284 IF (n4x6con /= 0) & 285 CALL rattle_4x6_int(molecule, particle_set, vel, dt) 286 ! Collective Variables 287 IF (ncolv%ntot /= 0) & 288 CALL rattle_colv_int(molecule, particle_set, vel, dt, & 289 irattle_int, cell, imass, int_max_sigma) 290 END DO Rattle_Intra_Loop 291 max_sigma = MAX(max_sigma, int_max_sigma) 292 CALL rattle_int_info(log_unit, i, irattle_int, max_sigma) 293 END DO 294 END DO MOL 295 ! Intermolecular Constraints 296 IF (do_ext_constraint) THEN 297 CALL update_temporary_set(group, vel=vel) 298 ! 3x3 299 IF (gci%ng3x3 /= 0) & 300 CALL rattle_3x3_ext(gci, particle_set, vel, dt) 301 ! 4x6 302 IF (gci%ng4x6 /= 0) & 303 CALL rattle_4x6_ext(gci, particle_set, vel, dt) 304 ! Collective Variables 305 IF (gci%ncolv%ntot /= 0) & 306 CALL rattle_colv_ext(gci, particle_set, vel, dt, & 307 irattle_ext, cell, imass, max_sigma) 308 CALL restore_temporary_set(particle_set, local_particles, vel=vel) 309 END IF 310 CALL rattle_ext_info(log_unit, irattle_ext, max_sigma) 311 END DO Rattle_Inter_Loop 312 CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, & 313 molecule_kind_set, group, "R") 314 CALL timestop(handle) 315 316 END SUBROUTINE rattle_control 317 318! ************************************************************************************************** 319!> \brief ... 320!> \param gci ... 321!> \param local_molecules ... 322!> \param molecule_set ... 323!> \param molecule_kind_set ... 324!> \param particle_set ... 325!> \param pos ... 326!> \param vel ... 327!> \param dt ... 328!> \param simpar ... 329!> \param roll_tol ... 330!> \param iroll ... 331!> \param vector_r ... 332!> \param vector_v ... 333!> \param group ... 334!> \param u ... 335!> \param cell ... 336!> \param local_particles ... 337!> \par History 338!> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints 339! ************************************************************************************************** 340 SUBROUTINE shake_roll_control(gci, local_molecules, molecule_set, & 341 molecule_kind_set, particle_set, pos, vel, dt, simpar, roll_tol, iroll, & 342 vector_r, vector_v, group, u, cell, local_particles) 343 344 TYPE(global_constraint_type), POINTER :: gci 345 TYPE(distribution_1d_type), POINTER :: local_molecules 346 TYPE(molecule_type), POINTER :: molecule_set(:) 347 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 348 TYPE(particle_type), POINTER :: particle_set(:) 349 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 350 REAL(KIND=dp), INTENT(IN) :: dt 351 TYPE(simpar_type), INTENT(IN) :: simpar 352 REAL(KIND=dp), INTENT(OUT) :: roll_tol 353 INTEGER, INTENT(INOUT) :: iroll 354 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector_r, vector_v 355 INTEGER, INTENT(IN) :: group 356 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 357 OPTIONAL :: u 358 TYPE(cell_type), POINTER :: cell 359 TYPE(distribution_1d_type), POINTER :: local_particles 360 361 CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_roll_control', & 362 routineP = moduleN//':'//routineN 363 364 INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, lagrange_mult, log_unit, & 365 n3x3con, n4x6con, nconstraint, nkind, nmol_per_kind, nvsitecon 366 LOGICAL :: do_ext_constraint, dump_lm 367 REAL(KIND=dp) :: int_max_sigma, mass, max_sigma, shake_tol 368 REAL(KIND=dp), DIMENSION(3, 3) :: r_shake, v_shake 369 REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass 370 TYPE(atomic_kind_type), POINTER :: atomic_kind 371 TYPE(colvar_counters) :: ncolv 372 TYPE(molecule_kind_type), POINTER :: molecule_kind 373 TYPE(molecule_type), POINTER :: molecule 374 375 CALL timeset(routineN, handle) 376 nkind = SIZE(molecule_kind_set) 377 shake_tol = simpar%shake_tol 378 log_unit = simpar%info_constraint 379 lagrange_mult = simpar%lagrange_multipliers 380 dump_lm = simpar%dump_lm 381 ! setting up for roll 382 IF (simpar%ensemble == npt_i_ensemble) THEN 383 CALL get_roll_matrix('SHAKE', r_shake, v_shake, vector_r, vector_v) 384 ELSE IF (simpar%ensemble == npt_f_ensemble) THEN 385 CALL get_roll_matrix('SHAKE', r_shake, v_shake, vector_r, vector_v, u) 386 END IF 387 DO k = 1, SIZE(pos, 2) 388 atomic_kind => particle_set(k)%atomic_kind 389 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 390 imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp)) 391 END DO 392 do_ext_constraint = (gci%ntot /= 0) 393 ishake_ext = 0 394 max_sigma = -1.0E+10_dp 395 Shake_Inter_Loop: DO WHILE (ABS(max_sigma) >= shake_tol) 396 max_sigma = 0.0_dp 397 ishake_ext = ishake_ext + 1 398 ! Intramolecular Constraints 399 MOL: DO ikind = 1, nkind 400 nmol_per_kind = local_molecules%n_el(ikind) 401 DO imol = 1, nmol_per_kind 402 i = local_molecules%list(ikind)%array(imol) 403 molecule => molecule_set(i) 404 molecule_kind => molecule%molecule_kind 405 CALL get_molecule_kind(molecule_kind, ncolv=ncolv, & 406 ng3x3=n3x3con, ng4x6=n4x6con, & 407 nconstraint=nconstraint, nvsite=nvsitecon) 408 IF (nconstraint == 0) CYCLE 409 ishake_int = 0 410 int_max_sigma = -1.0E+10_dp 411 Shake_Roll_Intra_Loop: DO WHILE (ABS(int_max_sigma) >= shake_tol) 412 int_max_sigma = 0.0_dp 413 ishake_int = ishake_int + 1 414 ! 3x3 415 IF (n3x3con /= 0) & 416 CALL shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, & 417 v_shake, dt, ishake_int, int_max_sigma) 418 ! 4x6 419 IF (n4x6con /= 0) & 420 CALL shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, & 421 dt, ishake_int, int_max_sigma) 422 ! Collective Variables 423 IF (ncolv%ntot /= 0) & 424 CALL shake_roll_colv_int(molecule, particle_set, pos, vel, r_shake, & 425 v_shake, dt, ishake_int, cell, imass, int_max_sigma) 426 END DO Shake_Roll_Intra_Loop 427 max_sigma = MAX(max_sigma, int_max_sigma) 428 CALL shake_int_info(log_unit, i, ishake_int, max_sigma) 429 ! Virtual Site 430 IF (nvsitecon /= 0) THEN 431 CPABORT("Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!") 432 END IF 433 END DO 434 END DO MOL 435 ! Intermolecular constraints 436 IF (do_ext_constraint) THEN 437 CALL update_temporary_set(group, pos=pos, vel=vel) 438 ! 3x3 439 IF (gci%ng3x3 /= 0) & 440 CALL shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, & 441 v_shake, dt, ishake_ext, max_sigma) 442 ! 4x6 443 IF (gci%ng4x6 /= 0) & 444 CALL shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, & 445 dt, ishake_ext, max_sigma) 446 ! Collective Variables 447 IF (gci%ncolv%ntot /= 0) & 448 CALL shake_roll_colv_ext(gci, particle_set, pos, vel, r_shake, & 449 v_shake, dt, ishake_ext, cell, imass, max_sigma) 450 ! Virtual Site 451 IF (gci%nvsite /= 0) & 452 CPABORT("Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!") 453 CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel) 454 END IF 455 CALL shake_ext_info(log_unit, ishake_ext, max_sigma) 456 END DO Shake_Inter_Loop 457 CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, & 458 molecule_kind_set, group, "S") 459 CALL check_tol(roll_tol, iroll, 'SHAKE', r_shake) 460 CALL timestop(handle) 461 462 END SUBROUTINE shake_roll_control 463 464! ************************************************************************************************** 465!> \brief ... 466!> \param gci ... 467!> \param local_molecules ... 468!> \param molecule_set ... 469!> \param molecule_kind_set ... 470!> \param particle_set ... 471!> \param vel ... 472!> \param dt ... 473!> \param simpar ... 474!> \param vector ... 475!> \param veps ... 476!> \param roll_tol ... 477!> \param iroll ... 478!> \param para_env ... 479!> \param u ... 480!> \param cell ... 481!> \param local_particles ... 482!> \par History 483!> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints 484! ************************************************************************************************** 485 SUBROUTINE rattle_roll_control(gci, local_molecules, molecule_set, & 486 molecule_kind_set, particle_set, vel, dt, simpar, vector, & 487 veps, roll_tol, iroll, para_env, u, cell, local_particles) 488 489 TYPE(global_constraint_type), POINTER :: gci 490 TYPE(distribution_1d_type), POINTER :: local_molecules 491 TYPE(molecule_type), POINTER :: molecule_set(:) 492 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 493 TYPE(particle_type), POINTER :: particle_set(:) 494 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 495 REAL(KIND=dp), INTENT(IN) :: dt 496 TYPE(simpar_type), INTENT(IN) :: simpar 497 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector 498 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: veps 499 REAL(KIND=dp), INTENT(OUT) :: roll_tol 500 INTEGER, INTENT(INOUT) :: iroll 501 TYPE(cp_para_env_type), INTENT(IN) :: para_env 502 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 503 OPTIONAL :: u 504 TYPE(cell_type), POINTER :: cell 505 TYPE(distribution_1d_type), POINTER :: local_particles 506 507 CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_roll_control', & 508 routineP = moduleN//':'//routineN 509 510 INTEGER :: handle, i, ikind, imol, irattle_ext, irattle_int, k, lagrange_mult, log_unit, & 511 n3x3con, n4x6con, nconstraint, nkind, nmol_per_kind 512 LOGICAL :: do_ext_constraint, dump_lm 513 REAL(KIND=dp) :: int_max_sigma, mass, max_sigma, & 514 rattle_tol 515 REAL(KIND=dp), DIMENSION(3, 3) :: r_rattle 516 REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass 517 TYPE(atomic_kind_type), POINTER :: atomic_kind 518 TYPE(colvar_counters) :: ncolv 519 TYPE(molecule_kind_type), POINTER :: molecule_kind 520 TYPE(molecule_type), POINTER :: molecule 521 522 CALL timeset(routineN, handle) 523 ! initialize locals 524 nkind = SIZE(molecule_kind_set) 525 rattle_tol = simpar%shake_tol 526 log_unit = simpar%info_constraint 527 lagrange_mult = simpar%lagrange_multipliers 528 dump_lm = simpar%dump_lm 529 ! setting up for roll 530 IF (simpar%ensemble == npt_i_ensemble) THEN 531 CALL get_roll_matrix('RATTLE', v_shake=r_rattle, vector_v=vector) 532 ELSE IF (simpar%ensemble == npt_f_ensemble) THEN 533 CALL get_roll_matrix('RATTLE', v_shake=r_rattle, vector_v=vector, u=u) 534 END IF 535 DO k = 1, SIZE(vel, 2) 536 atomic_kind => particle_set(k)%atomic_kind 537 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 538 imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp)) 539 END DO 540 do_ext_constraint = (gci%ntot /= 0) 541 irattle_ext = 0 542 max_sigma = -1.0E+10_dp 543 Rattle_Inter_Loop: DO WHILE (ABS(max_sigma) >= rattle_tol) 544 max_sigma = 0.0_dp 545 irattle_ext = irattle_ext + 1 546 ! Intramolecular Constraints 547 MOL: DO ikind = 1, nkind 548 nmol_per_kind = local_molecules%n_el(ikind) 549 DO imol = 1, nmol_per_kind 550 i = local_molecules%list(ikind)%array(imol) 551 molecule => molecule_set(i) 552 molecule_kind => molecule%molecule_kind 553 CALL get_molecule_kind(molecule_kind, ncolv=ncolv, & 554 ng3x3=n3x3con, ng4x6=n4x6con, & 555 nconstraint=nconstraint) 556 IF (nconstraint == 0) CYCLE 557 int_max_sigma = -1.0E+10_dp 558 irattle_int = 0 559 Rattle_Roll_Intramolecular: DO WHILE (ABS(int_max_sigma) >= rattle_tol) 560 int_max_sigma = 0.0_dp 561 irattle_int = irattle_int + 1 562 ! 3x3 563 IF (n3x3con /= 0) & 564 CALL rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, & 565 veps) 566 ! 4x6 567 IF (n4x6con /= 0) & 568 CALL rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, & 569 veps) 570 ! Collective Variables 571 IF (ncolv%ntot /= 0) & 572 CALL rattle_roll_colv_int(molecule, particle_set, vel, r_rattle, dt, & 573 irattle_int, veps, cell, imass, int_max_sigma) 574 END DO Rattle_Roll_Intramolecular 575 max_sigma = MAX(max_sigma, int_max_sigma) 576 CALL rattle_int_info(log_unit, i, irattle_int, max_sigma) 577 END DO 578 END DO MOL 579 ! Intermolecular Constraints 580 IF (do_ext_constraint) THEN 581 CALL update_temporary_set(para_env%group, vel=vel) 582 ! 3x3 583 IF (gci%ng3x3 /= 0) & 584 CALL rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, & 585 veps) 586 ! 4x6 587 IF (gci%ng4x6 /= 0) & 588 CALL rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, & 589 veps) 590 ! Collective Variables 591 IF (gci%ncolv%ntot /= 0) & 592 CALL rattle_roll_colv_ext(gci, particle_set, vel, r_rattle, dt, & 593 irattle_ext, veps, cell, imass, max_sigma) 594 CALL restore_temporary_set(particle_set, local_particles, vel=vel) 595 END IF 596 CALL rattle_ext_info(log_unit, irattle_ext, max_sigma) 597 END DO Rattle_Inter_Loop 598 CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, & 599 molecule_kind_set, para_env%group, "R") 600 CALL check_tol(roll_tol, iroll, 'RATTLE', veps=veps) 601 CALL timestop(handle) 602 END SUBROUTINE rattle_roll_control 603 604! ************************************************************************************************** 605!> \brief ... 606!> \param dump_lm ... 607!> \param log_unit ... 608!> \param local_molecules ... 609!> \param molecule_set ... 610!> \param gci ... 611!> \param molecule_kind_set ... 612!> \param group ... 613!> \param id_type ... 614!> \par History 615!> Teodoro Laino [tlaino] 2007 - Dumps lagrange multipliers 616! ************************************************************************************************** 617 SUBROUTINE dump_lagrange_mult(dump_lm, log_unit, local_molecules, molecule_set, gci, & 618 molecule_kind_set, group, id_type) 619 LOGICAL, INTENT(IN) :: dump_lm 620 INTEGER, INTENT(IN) :: log_unit 621 TYPE(distribution_1d_type), POINTER :: local_molecules 622 TYPE(molecule_type), POINTER :: molecule_set(:) 623 TYPE(global_constraint_type), POINTER :: gci 624 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 625 INTEGER, INTENT(IN) :: group 626 CHARACTER(LEN=1), INTENT(IN) :: id_type 627 628 CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_lagrange_mult', & 629 routineP = moduleN//':'//routineN 630 631 CHARACTER(LEN=default_string_length) :: label 632 INTEGER :: handle, i, ikind, imol, j, my_index, & 633 n3x3con, n4x6con, nconstraint, nkind 634 LOGICAL :: do_ext_constraint, do_int_constraint 635 REAL(KIND=dp), DIMENSION(:), POINTER :: lagr 636 TYPE(colvar_counters) :: ncolv 637 TYPE(molecule_kind_type), POINTER :: molecule_kind 638 TYPE(molecule_type), POINTER :: molecule 639 640 CALL timeset(routineN, handle) 641 ! Total number of intramolecular constraints (distributed) 642 CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, & 643 nconstraint=nconstraint) 644 do_int_constraint = (nconstraint > 0) 645 do_ext_constraint = (gci%ntot > 0) 646 IF (dump_lm .AND. (do_int_constraint .OR. do_ext_constraint)) THEN 647 nkind = SIZE(molecule_kind_set) 648 ALLOCATE (lagr(nconstraint)) 649 lagr = 0.0_dp 650 ! Dump lagrange multipliers for Intramolecular Constraints 651 my_index = 0 652 IF (do_int_constraint) THEN 653 MOL: DO ikind = 1, nkind 654 molecule_kind => molecule_kind_set(ikind) 655 CALL get_molecule_kind(molecule_kind, & 656 ncolv=ncolv, & 657 ng3x3=n3x3con, & 658 ng4x6=n4x6con) 659 DO imol = 1, molecule_kind%nmolecule 660 i = molecule_kind%molecule_list(imol) 661 IF (ANY(local_molecules%list(ikind)%array == i)) THEN 662 molecule => molecule_set(i) 663 ! Collective Variables 664 DO j = 1, ncolv%ntot 665 lagr(my_index + 1) = molecule%lci%lcolv(j)%lambda 666 my_index = my_index + 1 667 END DO 668 ! 3x3 669 DO j = 1, n3x3con 670 lagr(my_index + 1:my_index + 3) = molecule%lci%lg3x3(j)%lambda(:) 671 my_index = my_index + 3 672 END DO 673 ! 4x6 674 DO j = 1, n4x6con 675 lagr(my_index + 1:my_index + 6) = molecule%lci%lg4x6(j)%lambda(:) 676 my_index = my_index + 6 677 END DO 678 ELSE 679 my_index = my_index + ncolv%ntot + 3*n3x3con + 6*n4x6con 680 END IF 681 END DO 682 END DO MOL 683 CALL mp_sum(lagr, group) 684 END IF 685 ! Intermolecular constraints 686 IF (do_ext_constraint) THEN 687 CALL reallocate(lagr, 1, SIZE(lagr) + gci%ntot) 688 ! Collective Variables 689 DO j = 1, gci%ncolv%ntot 690 lagr(my_index + 1) = gci%lcolv(j)%lambda 691 my_index = my_index + 1 692 END DO 693 ! 3x3 694 DO j = 1, gci%ng3x3 695 lagr(my_index + 1:my_index + 3) = gci%lg3x3(j)%lambda(:) 696 my_index = my_index + 3 697 END DO 698 ! 4x6 699 DO j = 1, gci%ng4x6 700 lagr(my_index + 1:my_index + 6) = gci%lg4x6(j)%lambda(:) 701 my_index = my_index + 6 702 END DO 703 END IF 704 IF (log_unit > 0) THEN 705 IF (id_type == "S") THEN 706 label = "Shake Lagrangian Multipliers:" 707 ELSEIF (id_type == "R") THEN 708 label = "Rattle Lagrangian Multipliers:" 709 ELSE 710 CPABORT("") 711 END IF 712 WRITE (log_unit, FMT='(A,T40,4F15.9)') TRIM(label), lagr(1:MIN(4, SIZE(lagr))) 713 DO j = 5, SIZE(lagr), 4 714 WRITE (log_unit, FMT='(T40,4F15.9)') lagr(j:MIN(j + 3, SIZE(lagr))) 715 END DO 716 END IF 717 DEALLOCATE (lagr) 718 END IF 719 CALL timestop(handle) 720 721 END SUBROUTINE dump_lagrange_mult 722 723! ************************************************************************************************** 724!> \brief Dumps convergence info about shake - intramolecular constraint loop 725!> \param log_unit ... 726!> \param i ... 727!> \param ishake_int ... 728!> \param max_sigma ... 729!> \par History 730!> Teodoro Laino [tlaino] 2007 - University of Zurich 731! ************************************************************************************************** 732 SUBROUTINE shake_int_info(log_unit, i, ishake_int, max_sigma) 733 INTEGER, INTENT(IN) :: log_unit, i, ishake_int 734 REAL(KIND=dp), INTENT(IN) :: max_sigma 735 736 CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_int_info', routineP = moduleN//':'//routineN 737 738 IF (log_unit > 0) THEN 739 ! Dump info if requested 740 WRITE (log_unit, '("SHAKE_INFO|",2X,2(A,I6),A,F15.9)') & 741 "Molecule Nr.:", i, " Nr. Iterations:", ishake_int, " Max. Err.:", max_sigma 742 END IF 743 ! Notify a not converged SHAKE 744 IF (ishake_int > Max_Shake_Iter) & 745 CALL cp_warn(__LOCATION__, & 746 "Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// & 747 " intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// & 748 ". CP2K continues but results could be meaningless. ") 749 END SUBROUTINE shake_int_info 750 751! ************************************************************************************************** 752!> \brief Dumps convergence info about shake - intermolecular constraint loop 753!> \param log_unit ... 754!> \param ishake_ext ... 755!> \param max_sigma ... 756!> \par History 757!> Teodoro Laino [tlaino] 2007 - University of Zurich 758! ************************************************************************************************** 759 SUBROUTINE shake_ext_info(log_unit, ishake_ext, max_sigma) 760 INTEGER, INTENT(IN) :: log_unit, ishake_ext 761 REAL(KIND=dp), INTENT(IN) :: max_sigma 762 763 CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_ext_info', routineP = moduleN//':'//routineN 764 765 IF (log_unit > 0) THEN 766 ! Dump info if requested 767 WRITE (log_unit, '("SHAKE_INFO|",2X,A,I6,A,F15.9)') & 768 "External Shake Nr. Iterations:", ishake_ext, & 769 " Max. Err.:", max_sigma 770 END IF 771 ! Notify a not converged SHAKE 772 IF (ishake_ext > Max_Shake_Iter) & 773 CALL cp_warn(__LOCATION__, & 774 "Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// & 775 " intermolecular constraint. CP2K continues but results could be meaningless. ") 776 END SUBROUTINE shake_ext_info 777 778! ************************************************************************************************** 779!> \brief Dumps convergence info about rattle - intramolecular constraint loop 780!> \param log_unit ... 781!> \param i ... 782!> \param irattle_int ... 783!> \param max_sigma ... 784!> \par History 785!> Teodoro Laino [tlaino] 2007 - University of Zurich 786! ************************************************************************************************** 787 SUBROUTINE rattle_int_info(log_unit, i, irattle_int, max_sigma) 788 INTEGER, INTENT(IN) :: log_unit, i, irattle_int 789 REAL(KIND=dp), INTENT(IN) :: max_sigma 790 791 CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_int_info', & 792 routineP = moduleN//':'//routineN 793 794 IF (log_unit > 0) THEN 795 ! Dump info if requested 796 WRITE (log_unit, '("RATTLE_INFO|",1X,2(A,I6),A,F15.9)') & 797 "Molecule Nr.:", i, " Nr. Iterations:", irattle_int, " Max. Err.:", max_sigma 798 END IF 799 ! Notify a not converged RATTLE 800 IF (irattle_int > Max_shake_Iter) & 801 CALL cp_warn(__LOCATION__, & 802 "Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// & 803 " intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// & 804 ". CP2K continues but results could be meaningless. ") 805 END SUBROUTINE rattle_int_info 806 807! ************************************************************************************************** 808!> \brief Dumps convergence info about rattle - intermolecular constraint loop 809!> \param log_unit ... 810!> \param irattle_ext ... 811!> \param max_sigma ... 812!> \par History 813!> Teodoro Laino [tlaino] 2007 - University of Zurich 814! ************************************************************************************************** 815 SUBROUTINE rattle_ext_info(log_unit, irattle_ext, max_sigma) 816 INTEGER, INTENT(IN) :: log_unit, irattle_ext 817 REAL(KIND=dp), INTENT(IN) :: max_sigma 818 819 CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_ext_info', & 820 routineP = moduleN//':'//routineN 821 822 IF (log_unit > 0) THEN 823 ! Dump info if requested 824 WRITE (log_unit, '("RATTLE_INFO|",1X,A,I6,A,F15.9)') & 825 "External Rattle Nr. Iterations:", irattle_ext, & 826 " Max. Err.:", max_sigma 827 END IF 828 ! Notify a not converged RATTLE 829 IF (irattle_ext > Max_shake_Iter) & 830 CALL cp_warn(__LOCATION__, & 831 "Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// & 832 " intermolecular constraint. CP2K continues but results could be meaningless. ") 833 END SUBROUTINE rattle_ext_info 834 835! ************************************************************************************************** 836!> \brief Updates the TARGET of the COLLECTIVE constraints if the growth speed 837!> is different from zero. 838!> \param gci ... 839!> \param local_molecules ... 840!> \param molecule_set ... 841!> \param molecule_kind_set ... 842!> \param dt ... 843!> \param root_section ... 844!> \date 02.2008 845!> \author Teodoro Laino [tlaino] - University of Zurich 846! ************************************************************************************************** 847 SUBROUTINE shake_update_targets(gci, local_molecules, molecule_set, & 848 molecule_kind_set, dt, root_section) 849 850 TYPE(global_constraint_type), POINTER :: gci 851 TYPE(distribution_1d_type), POINTER :: local_molecules 852 TYPE(molecule_type), POINTER :: molecule_set(:) 853 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 854 REAL(kind=dp), INTENT(in) :: dt 855 TYPE(section_vals_type), POINTER :: root_section 856 857 CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_update_targets', & 858 routineP = moduleN//':'//routineN 859 860 INTEGER :: handle, i, ikind, imol, nkind, & 861 nmol_per_kind 862 LOGICAL :: do_ext_constraint 863 TYPE(colvar_counters) :: ncolv 864 TYPE(molecule_kind_type), POINTER :: molecule_kind 865 TYPE(molecule_type), POINTER :: molecule 866 TYPE(section_vals_type), POINTER :: motion_section 867 868 CALL timeset(routineN, handle) 869 motion_section => section_vals_get_subs_vals(root_section, "MOTION") 870 nkind = SIZE(molecule_kind_set) 871 do_ext_constraint = (gci%ntot /= 0) 872 ! Intramolecular Constraints 873 MOL: DO ikind = 1, nkind 874 nmol_per_kind = local_molecules%n_el(ikind) 875 DO imol = 1, nmol_per_kind 876 i = local_molecules%list(ikind)%array(imol) 877 molecule => molecule_set(i) 878 molecule_kind => molecule%molecule_kind 879 CALL get_molecule_kind(molecule_kind, ncolv=ncolv) 880 881 ! Updating TARGETS for Collective Variables only 882 IF (ncolv%ntot /= 0) CALL shake_update_colv_int(molecule, dt, motion_section) 883 END DO 884 END DO MOL 885 ! Intermolecular constraints 886 IF (do_ext_constraint) THEN 887 ! Collective Variables 888 IF (gci%ncolv%ntot /= 0) CALL shake_update_colv_ext(gci, dt, motion_section) 889 END IF 890 CALL timestop(handle) 891 END SUBROUTINE shake_update_targets 892 893END MODULE constraint 894