1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Contains routines useful for the application of constraints during MD 8!> \par History 9!> none 10! ************************************************************************************************** 11MODULE constraint_util 12 USE cell_types, ONLY: cell_type 13 USE colvar_methods, ONLY: colvar_eval_mol_f 14 USE distribution_1d_types, ONLY: distribution_1d_type 15 USE kinds, ONLY: dp 16 USE mathlib, ONLY: matmul_3x3,& 17 transpose_3d 18 USE message_passing, ONLY: mp_sum 19 USE molecule_kind_types, ONLY: colvar_constraint_type,& 20 fixd_constraint_type,& 21 g3x3_constraint_type,& 22 g4x6_constraint_type,& 23 get_molecule_kind,& 24 molecule_kind_type 25 USE molecule_types, ONLY: get_molecule,& 26 global_constraint_type,& 27 local_colvar_constraint_type,& 28 local_constraint_type,& 29 local_g3x3_constraint_type,& 30 local_g4x6_constraint_type,& 31 molecule_type 32 USE particle_types, ONLY: particle_type 33 USE virial_types, ONLY: virial_type 34#include "./base/base_uses.f90" 35 36 IMPLICIT NONE 37 38 PRIVATE 39 PUBLIC :: getold, & 40 pv_constraint, & 41 check_tol, & 42 get_roll_matrix, & 43 restore_temporary_set, & 44 update_temporary_set 45 46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_util' 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief saves all of the old variables 52!> \param gci ... 53!> \param local_molecules ... 54!> \param molecule_set ... 55!> \param molecule_kind_set ... 56!> \param particle_set ... 57!> \param cell ... 58!> \par History 59!> none 60! ************************************************************************************************** 61 SUBROUTINE getold(gci, local_molecules, molecule_set, molecule_kind_set, & 62 particle_set, cell) 63 64 TYPE(global_constraint_type), POINTER :: gci 65 TYPE(distribution_1d_type), POINTER :: local_molecules 66 TYPE(molecule_type), POINTER :: molecule_set(:) 67 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 68 TYPE(particle_type), POINTER :: particle_set(:) 69 TYPE(cell_type), POINTER :: cell 70 71 CHARACTER(LEN=*), PARAMETER :: routineN = 'getold', routineP = moduleN//':'//routineN 72 73 INTEGER :: first_atom, i, ikind, imol, n3x3con, & 74 n4x6con, nkind, nmol_per_kind 75 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 76 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 77 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 78 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 79 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 80 TYPE(local_constraint_type), POINTER :: lci 81 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 82 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 83 TYPE(molecule_kind_type), POINTER :: molecule_kind 84 TYPE(molecule_type), POINTER :: molecule 85 86 NULLIFY (fixd_list) 87 nkind = SIZE(molecule_kind_set) 88 ! Intramolecular constraints 89 MOL: DO ikind = 1, nkind 90 nmol_per_kind = local_molecules%n_el(ikind) 91 DO imol = 1, nmol_per_kind 92 i = local_molecules%list(ikind)%array(imol) 93 molecule => molecule_set(i) 94 molecule_kind => molecule%molecule_kind 95 CALL get_molecule_kind(molecule_kind, ng3x3=n3x3con, ng4x6=n4x6con, & 96 colv_list=colv_list, g3x3_list=g3x3_list, g4x6_list=g4x6_list, & 97 fixd_list=fixd_list) 98 CALL get_molecule(molecule, lci=lci) 99 IF (.NOT. ASSOCIATED(lci)) CYCLE 100 CALL get_molecule(molecule, first_atom=first_atom, & 101 lcolv=lcolv, lg3x3=lg3x3, lg4x6=lg4x6) 102 CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, & 103 lcolv, lg3x3, lg4x6, first_atom, particle_set, cell) 104 END DO 105 END DO MOL 106 ! Intermolecular constraints 107 IF (gci%ntot > 0) THEN 108 n3x3con = gci%ng3x3 109 n4x6con = gci%ng4x6 110 colv_list => gci%colv_list 111 g3x3_list => gci%g3x3_list 112 g4x6_list => gci%g4x6_list 113 fixd_list => gci%fixd_list 114 lcolv => gci%lcolv 115 lg3x3 => gci%lg3x3 116 lg4x6 => gci%lg4x6 117 CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, & 118 lcolv, lg3x3, lg4x6, 1, particle_set, cell) 119 END IF 120 END SUBROUTINE getold 121 122! ************************************************************************************************** 123!> \brief saves all of the old variables - Low Level 124!> \param n3x3con ... 125!> \param n4x6con ... 126!> \param colv_list ... 127!> \param g3x3_list ... 128!> \param g4x6_list ... 129!> \param fixd_list ... 130!> \param lcolv ... 131!> \param lg3x3 ... 132!> \param lg4x6 ... 133!> \param first_atom ... 134!> \param particle_set ... 135!> \param cell ... 136!> \par History 137!> none 138! ************************************************************************************************** 139 SUBROUTINE getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, & 140 lcolv, lg3x3, lg4x6, first_atom, particle_set, cell) 141 142 INTEGER, INTENT(IN) :: n3x3con, n4x6con 143 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 144 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 145 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 146 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 147 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 148 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 149 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 150 INTEGER, INTENT(IN) :: first_atom 151 TYPE(particle_type), POINTER :: particle_set(:) 152 TYPE(cell_type), POINTER :: cell 153 154 CHARACTER(LEN=*), PARAMETER :: routineN = 'getold_low', routineP = moduleN//':'//routineN 155 156 INTEGER :: iconst, index 157 158 IF (ASSOCIATED(colv_list)) THEN 159 ! Collective constraints 160 DO iconst = 1, SIZE(colv_list) 161 CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, & 162 particles=particle_set, fixd_list=fixd_list) 163 ENDDO 164 END IF 165 ! 3x3 constraints 166 DO iconst = 1, n3x3con 167 index = g3x3_list(iconst)%a + first_atom - 1 168 lg3x3(iconst)%ra_old = particle_set(index)%r 169 index = g3x3_list(iconst)%b + first_atom - 1 170 lg3x3(iconst)%rb_old = particle_set(index)%r 171 index = g3x3_list(iconst)%c + first_atom - 1 172 lg3x3(iconst)%rc_old = particle_set(index)%r 173 ENDDO 174 ! 4x6 constraints 175 DO iconst = 1, n4x6con 176 index = g4x6_list(iconst)%a + first_atom - 1 177 lg4x6(iconst)%ra_old = particle_set(index)%r 178 index = g4x6_list(iconst)%b + first_atom - 1 179 lg4x6(iconst)%rb_old = particle_set(index)%r 180 index = g4x6_list(iconst)%c + first_atom - 1 181 lg4x6(iconst)%rc_old = particle_set(index)%r 182 index = g4x6_list(iconst)%d + first_atom - 1 183 lg4x6(iconst)%rd_old = particle_set(index)%r 184 ENDDO 185 186 END SUBROUTINE getold_low 187 188! ************************************************************************************************** 189!> \brief ... 190!> \param gci ... 191!> \param local_molecules ... 192!> \param molecule_set ... 193!> \param molecule_kind_set ... 194!> \param particle_set ... 195!> \param virial ... 196!> \param group ... 197!> \par History 198!> none 199! ************************************************************************************************** 200 SUBROUTINE pv_constraint(gci, local_molecules, molecule_set, molecule_kind_set, & 201 particle_set, virial, group) 202 203 TYPE(global_constraint_type), POINTER :: gci 204 TYPE(distribution_1d_type), POINTER :: local_molecules 205 TYPE(molecule_type), POINTER :: molecule_set(:) 206 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 207 TYPE(particle_type), POINTER :: particle_set(:) 208 TYPE(virial_type), INTENT(INOUT) :: virial 209 INTEGER, INTENT(IN) :: group 210 211 CHARACTER(LEN=*), PARAMETER :: routineN = 'pv_constraint', routineP = moduleN//':'//routineN 212 213 INTEGER :: first_atom, i, ikind, imol, ng3x3, & 214 ng4x6, nkind, nmol_per_kind 215 REAL(KIND=dp) :: pv(3, 3) 216 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 217 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 218 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 219 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 220 TYPE(local_constraint_type), POINTER :: lci 221 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 222 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 223 TYPE(molecule_kind_type), POINTER :: molecule_kind 224 TYPE(molecule_type), POINTER :: molecule 225 226 pv = 0.0_dp 227 nkind = SIZE(molecule_kind_set) 228 ! Intramolecular Constraints 229 MOL: DO ikind = 1, nkind 230 nmol_per_kind = local_molecules%n_el(ikind) 231 DO imol = 1, nmol_per_kind 232 i = local_molecules%list(ikind)%array(imol) 233 molecule => molecule_set(i) 234 molecule_kind => molecule%molecule_kind 235 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, & 236 ng4x6=ng4x6, g3x3_list=g3x3_list, g4x6_list=g4x6_list, & 237 colv_list=colv_list) 238 CALL get_molecule(molecule, lci=lci) 239 IF (.NOT. ASSOCIATED(lci)) CYCLE 240 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3, & 241 lg4x6=lg4x6, lcolv=lcolv) 242 CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, & 243 first_atom, lg3x3, lg4x6, lcolv, particle_set, pv) 244 END DO 245 END DO MOL 246 ! Intermolecular constraints 247 IF (gci%ntot > 0) THEN 248 ng3x3 = gci%ng3x3 249 ng4x6 = gci%ng4x6 250 colv_list => gci%colv_list 251 g3x3_list => gci%g3x3_list 252 g4x6_list => gci%g4x6_list 253 lcolv => gci%lcolv 254 lg3x3 => gci%lg3x3 255 lg4x6 => gci%lg4x6 256 CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, & 257 1, lg3x3, lg4x6, lcolv, particle_set, pv) 258 END IF 259 CALL mp_sum(pv, group) 260 ! Symmetrize PV 261 pv(1, 2) = 0.5_dp*(pv(1, 2) + pv(2, 1)) 262 pv(2, 1) = pv(1, 2) 263 pv(1, 3) = 0.5_dp*(pv(1, 3) + pv(3, 1)) 264 pv(3, 1) = pv(1, 3) 265 pv(3, 2) = 0.5_dp*(pv(3, 2) + pv(2, 3)) 266 pv(2, 3) = pv(3, 2) 267 ! Store in virial type 268 virial%pv_constraint = pv 269 270 END SUBROUTINE pv_constraint 271 272! ************************************************************************************************** 273!> \brief ... 274!> \param ng3x3 ... 275!> \param ng4x6 ... 276!> \param g3x3_list ... 277!> \param g4x6_list ... 278!> \param colv_list ... 279!> \param first_atom ... 280!> \param lg3x3 ... 281!> \param lg4x6 ... 282!> \param lcolv ... 283!> \param particle_set ... 284!> \param pv ... 285!> \par History 286!> none 287! ************************************************************************************************** 288 SUBROUTINE pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, & 289 first_atom, lg3x3, lg4x6, lcolv, particle_set, pv) 290 291 INTEGER, INTENT(IN) :: ng3x3, ng4x6 292 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 293 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 294 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 295 INTEGER, INTENT(IN) :: first_atom 296 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 297 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 298 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 299 TYPE(particle_type), POINTER :: particle_set(:) 300 REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv 301 302 CHARACTER(LEN=*), PARAMETER :: routineN = 'pv_constraint_low', & 303 routineP = moduleN//':'//routineN 304 305 INTEGER :: iconst, index_a, index_b, index_c, & 306 index_d 307 REAL(KIND=dp) :: fc1(3), fc2(3), fc3(3), fc4(3), & 308 lambda_3x3(3), lambda_4x6(6) 309 310 IF (ASSOCIATED(colv_list)) THEN 311 ! Colvar Constraints 312 DO iconst = 1, SIZE(colv_list) 313 CALL pv_colv_eval(pv, lcolv(iconst), particle_set) 314 END DO 315 END IF 316 ! 3x3 317 DO iconst = 1, ng3x3 318 ! pv gets updated with FULL multiplier 319 lambda_3x3 = lg3x3(iconst)%lambda 320 321 fc1 = lambda_3x3(1)*lg3x3(iconst)%fa + & 322 lambda_3x3(2)*lg3x3(iconst)%fb 323 fc2 = -lambda_3x3(1)*lg3x3(iconst)%fa + & 324 lambda_3x3(3)*lg3x3(iconst)%fc 325 fc3 = -lambda_3x3(2)*lg3x3(iconst)%fb - & 326 lambda_3x3(3)*lg3x3(iconst)%fc 327 index_a = g3x3_list(iconst)%a + first_atom - 1 328 index_b = g3x3_list(iconst)%b + first_atom - 1 329 index_c = g3x3_list(iconst)%c + first_atom - 1 330 331 !pv(1,1) 332 pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1) 333 pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1) 334 pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1) 335 !pv(1,2) 336 pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2) 337 pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2) 338 pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2) 339 !pv(1,3) 340 pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3) 341 pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3) 342 pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3) 343 !pv(2,1) 344 pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1) 345 pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1) 346 pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1) 347 !pv(2,2) 348 pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2) 349 pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2) 350 pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2) 351 !pv(2,3) 352 pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3) 353 pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3) 354 pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3) 355 !pv(3,1) 356 pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1) 357 pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1) 358 pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1) 359 !pv(3,2) 360 pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2) 361 pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2) 362 pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2) 363 !pv(3,3) 364 pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3) 365 pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3) 366 pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3) 367 END DO 368 369 ! 4x6 370 DO iconst = 1, ng4x6 371 ! pv gets updated with FULL multiplier 372 lambda_4x6 = lg4x6(iconst)%lambda 373 374 fc1 = lambda_4x6(1)*lg4x6(iconst)%fa + & 375 lambda_4x6(2)*lg4x6(iconst)%fb + & 376 lambda_4x6(3)*lg4x6(iconst)%fc 377 fc2 = -lambda_4x6(1)*lg4x6(iconst)%fa + & 378 lambda_4x6(4)*lg4x6(iconst)%fd + & 379 lambda_4x6(5)*lg4x6(iconst)%fe 380 fc3 = -lambda_4x6(2)*lg4x6(iconst)%fb - & 381 lambda_4x6(4)*lg4x6(iconst)%fd + & 382 lambda_4x6(6)*lg4x6(iconst)%ff 383 fc4 = -lambda_4x6(3)*lg4x6(iconst)%fc - & 384 lambda_4x6(5)*lg4x6(iconst)%fe - & 385 lambda_4x6(6)*lg4x6(iconst)%ff 386 index_a = g4x6_list(iconst)%a + first_atom - 1 387 index_b = g4x6_list(iconst)%b + first_atom - 1 388 index_c = g4x6_list(iconst)%c + first_atom - 1 389 index_d = g4x6_list(iconst)%d + first_atom - 1 390 391 !pv(1,1) 392 pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1) 393 pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1) 394 pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1) 395 pv(1, 1) = pv(1, 1) + fc4(1)*particle_set(index_d)%r(1) 396 !pv(1,2) 397 pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2) 398 pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2) 399 pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2) 400 pv(1, 2) = pv(1, 2) + fc4(1)*particle_set(index_d)%r(2) 401 !pv(1,3) 402 pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3) 403 pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3) 404 pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3) 405 pv(1, 3) = pv(1, 3) + fc4(1)*particle_set(index_d)%r(3) 406 !pv(2,1) 407 pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1) 408 pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1) 409 pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1) 410 pv(2, 1) = pv(2, 1) + fc4(2)*particle_set(index_d)%r(1) 411 !pv(2,2) 412 pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2) 413 pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2) 414 pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2) 415 pv(2, 2) = pv(2, 2) + fc4(2)*particle_set(index_d)%r(2) 416 !pv(2,3) 417 pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3) 418 pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3) 419 pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3) 420 pv(2, 3) = pv(2, 3) + fc4(2)*particle_set(index_d)%r(3) 421 !pv(3,1) 422 pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1) 423 pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1) 424 pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1) 425 pv(3, 1) = pv(3, 1) + fc4(3)*particle_set(index_d)%r(1) 426 !pv(3,2) 427 pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2) 428 pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2) 429 pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2) 430 pv(3, 2) = pv(3, 2) + fc4(3)*particle_set(index_d)%r(2) 431 !pv(3,3) 432 pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3) 433 pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3) 434 pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3) 435 pv(3, 3) = pv(3, 3) + fc4(3)*particle_set(index_d)%r(3) 436 END DO 437 438 END SUBROUTINE pv_constraint_low 439 440! ************************************************************************************************** 441!> \brief ... 442!> \param pv ... 443!> \param lcolv ... 444!> \param particle_set ... 445!> \par History 446!> none 447! ************************************************************************************************** 448 SUBROUTINE pv_colv_eval(pv, lcolv, particle_set) 449 REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv 450 TYPE(local_colvar_constraint_type), INTENT(IN) :: lcolv 451 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 452 453 INTEGER :: i, iatm, ind, j 454 REAL(KIND=dp) :: lambda, tmp 455 REAL(KIND=dp), DIMENSION(3) :: f 456 457 DO iatm = 1, SIZE(lcolv%colvar_old%i_atom) 458 ind = lcolv%colvar_old%i_atom(iatm) 459 f = -lcolv%colvar_old%dsdr(:, iatm) 460 ! pv gets updated with FULL multiplier 461 lambda = lcolv%lambda 462 DO i = 1, 3 463 tmp = lambda*particle_set(ind)%r(i) 464 DO j = 1, 3 465 pv(j, i) = pv(j, i) + f(j)*tmp 466 END DO 467 END DO 468 END DO 469 470 END SUBROUTINE pv_colv_eval 471 472! ************************************************************************************************** 473!> \brief ... 474!> \param roll_tol ... 475!> \param iroll ... 476!> \param char ... 477!> \param matrix ... 478!> \param veps ... 479!> \par History 480!> none 481! ************************************************************************************************** 482 SUBROUTINE check_tol(roll_tol, iroll, char, matrix, veps) 483 484 REAL(KIND=dp), INTENT(OUT) :: roll_tol 485 INTEGER, INTENT(INOUT) :: iroll 486 CHARACTER(LEN=*), INTENT(IN) :: char 487 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 488 OPTIONAL :: matrix, veps 489 490 CHARACTER(LEN=*), PARAMETER :: routineN = 'check_tol', routineP = moduleN//':'//routineN 491 492 REAL(KIND=dp) :: local_tol 493 REAL(KIND=dp), DIMENSION(3, 3) :: diff_rattle, diff_shake 494 REAL(KIND=dp), DIMENSION(3, 3), SAVE :: matrix_old, veps_old 495 496 SELECT CASE (char) 497 CASE ('SHAKE') 498 IF (iroll == 1) THEN 499 matrix_old = matrix 500 roll_tol = -1.E10_dp 501 ELSE 502 roll_tol = 0.0_dp 503 diff_shake = ABS(matrix_old - matrix) 504 local_tol = MAXVAL(diff_shake) 505 roll_tol = MAX(roll_tol, local_tol) 506 matrix_old = matrix 507 END IF 508 iroll = iroll + 1 509 CASE ('RATTLE') 510 IF (iroll == 1) THEN 511 veps_old = veps 512 roll_tol = -1.E+10_dp 513 ELSE 514 roll_tol = 0.0_dp 515 ! compute tolerance on veps 516 diff_rattle = ABS(veps - veps_old) 517 local_tol = MAXVAL(diff_rattle) 518 roll_tol = MAX(roll_tol, local_tol) 519 veps_old = veps 520 END IF 521 iroll = iroll + 1 522 END SELECT 523 524 END SUBROUTINE check_tol 525 526! ************************************************************************************************** 527!> \brief ... 528!> \param char ... 529!> \param r_shake ... 530!> \param v_shake ... 531!> \param vector_r ... 532!> \param vector_v ... 533!> \param u ... 534!> \par History 535!> none 536! ************************************************************************************************** 537 SUBROUTINE get_roll_matrix(char, r_shake, v_shake, vector_r, vector_v, u) 538 539 CHARACTER(len=*), INTENT(IN) :: char 540 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), & 541 OPTIONAL :: r_shake, v_shake 542 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: vector_r, vector_v 543 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 544 OPTIONAL :: u 545 546 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_roll_matrix', & 547 routineP = moduleN//':'//routineN 548 549 INTEGER :: i 550 REAL(KIND=dp), DIMENSION(3, 3) :: diag 551 552 IF (PRESENT(r_shake)) r_shake = 0.0_dp 553 IF (PRESENT(v_shake)) v_shake = 0.0_dp 554 diag = 0.0_dp 555 556 SELECT CASE (char) 557 CASE ('SHAKE') 558 IF (PRESENT(u) .AND. PRESENT(vector_v) .AND. & 559 PRESENT(vector_r)) THEN 560 diag(1, 1) = vector_r(1) 561 diag(2, 2) = vector_r(2) 562 diag(3, 3) = vector_r(3) 563 r_shake = MATMUL_3X3(MATMUL_3X3(u, diag), TRANSPOSE_3D(u)) 564 diag(1, 1) = vector_v(1) 565 diag(2, 2) = vector_v(2) 566 diag(3, 3) = vector_v(3) 567 v_shake = MATMUL_3X3(MATMUL_3X3(u, diag), TRANSPOSE_3D(u)) 568 diag = MATMUL_3x3(r_shake, v_shake) 569 r_shake = diag 570 ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. & 571 PRESENT(vector_r)) THEN 572 DO i = 1, 3 573 r_shake(i, i) = vector_r(i)*vector_v(i) 574 v_shake(i, i) = vector_v(i) 575 END DO 576 ELSE 577 CPABORT("Not sufficient parameters") 578 END IF 579 CASE ('RATTLE') 580 IF (PRESENT(u) .AND. PRESENT(vector_v)) THEN 581 diag(1, 1) = vector_v(1) 582 diag(2, 2) = vector_v(2) 583 diag(3, 3) = vector_v(3) 584 v_shake = MATMUL_3x3(MATMUL_3X3(u, diag), TRANSPOSE_3D(u)) 585 ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN 586 DO i = 1, 3 587 v_shake(i, i) = vector_v(i) 588 END DO 589 ELSE 590 CPABORT("Not sufficient parameters") 591 END IF 592 END SELECT 593 594 END SUBROUTINE get_roll_matrix 595 596! ************************************************************************************************** 597!> \brief ... 598!> \param particle_set ... 599!> \param local_particles ... 600!> \param pos ... 601!> \param vel ... 602!> \par History 603!> Teodoro Laino [tlaino] 2007 604! ************************************************************************************************** 605 SUBROUTINE restore_temporary_set(particle_set, local_particles, pos, vel) 606 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 607 TYPE(distribution_1d_type), POINTER :: local_particles 608 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 609 OPTIONAL :: pos, vel 610 611 CHARACTER(LEN=*), PARAMETER :: routineN = 'restore_temporary_set', & 612 routineP = moduleN//':'//routineN 613 614 INTEGER :: iparticle, iparticle_kind, & 615 iparticle_local, nparticle_local 616 LOGICAL, ALLOCATABLE, DIMENSION(:) :: wrk 617 618 ALLOCATE (wrk(SIZE(particle_set))) 619 wrk = .TRUE. 620 DO iparticle_kind = 1, SIZE(local_particles%n_el) 621 nparticle_local = local_particles%n_el(iparticle_kind) 622 DO iparticle_local = 1, nparticle_local 623 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local) 624 wrk(iparticle) = .FALSE. 625 END DO 626 END DO 627 IF (PRESENT(vel)) THEN 628 DO iparticle = 1, SIZE(particle_set) 629 IF (wrk(iparticle)) THEN 630 vel(:, iparticle) = 0.0_dp 631 END IF 632 END DO 633 END IF 634 IF (PRESENT(pos)) THEN 635 DO iparticle = 1, SIZE(particle_set) 636 IF (wrk(iparticle)) THEN 637 pos(:, iparticle) = 0.0_dp 638 END IF 639 END DO 640 END IF 641 DEALLOCATE (wrk) 642 END SUBROUTINE restore_temporary_set 643 644! ************************************************************************************************** 645!> \brief ... 646!> \param group ... 647!> \param pos ... 648!> \param vel ... 649!> \par History 650!> Teodoro Laino [tlaino] 2007 651! ************************************************************************************************** 652 SUBROUTINE update_temporary_set(group, pos, vel) 653 INTEGER, INTENT(IN) :: group 654 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 655 OPTIONAL :: pos, vel 656 657 IF (PRESENT(pos)) THEN 658 CALL mp_sum(pos, group) 659 END IF 660 IF (PRESENT(vel)) THEN 661 CALL mp_sum(vel, group) 662 END IF 663 END SUBROUTINE update_temporary_set 664 665END MODULE constraint_util 666