1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> none 9! ************************************************************************************************** 10MODULE constraint_3x3 11 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE constraint_fxd, ONLY: check_fixed_atom_cns_g3x3 15 USE kinds, ONLY: dp 16 USE linear_systems, ONLY: solve_system 17 USE mathlib, ONLY: dotprod_3d,& 18 matvec_3x3 19 USE molecule_kind_types, ONLY: fixd_constraint_type,& 20 g3x3_constraint_type,& 21 get_molecule_kind,& 22 molecule_kind_type 23 USE molecule_types, ONLY: get_molecule,& 24 global_constraint_type,& 25 local_g3x3_constraint_type,& 26 molecule_type 27 USE particle_types, ONLY: particle_type 28#include "./base/base_uses.f90" 29 30 IMPLICIT NONE 31 32 PRIVATE 33 PUBLIC :: shake_3x3_int, & 34 rattle_3x3_int, & 35 shake_roll_3x3_int, & 36 rattle_roll_3x3_int, & 37 shake_3x3_ext, & 38 rattle_3x3_ext, & 39 shake_roll_3x3_ext, & 40 rattle_roll_3x3_ext 41 42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_3x3' 43 44CONTAINS 45 46! ************************************************************************************************** 47!> \brief ... 48!> \param molecule ... 49!> \param particle_set ... 50!> \param pos ... 51!> \param vel ... 52!> \param dt ... 53!> \param ishake ... 54!> \param max_sigma ... 55!> \par History 56!> none 57! ************************************************************************************************** 58 SUBROUTINE shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake, & 59 max_sigma) 60 61 TYPE(molecule_type), POINTER :: molecule 62 TYPE(particle_type), POINTER :: particle_set(:) 63 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 64 REAL(kind=dp), INTENT(in) :: dt 65 INTEGER, INTENT(IN) :: ishake 66 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 67 68 INTEGER :: first_atom, ng3x3 69 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 70 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 71 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 72 TYPE(molecule_kind_type), POINTER :: molecule_kind 73 74 NULLIFY (fixd_list) 75 molecule_kind => molecule%molecule_kind 76 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, & 77 g3x3_list=g3x3_list, fixd_list=fixd_list) 78 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3) 79 ! Real Shake 80 CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 81 particle_set, pos, vel, dt, ishake, max_sigma) 82 83 END SUBROUTINE shake_3x3_int 84 85! ************************************************************************************************** 86!> \brief ... 87!> \param molecule ... 88!> \param particle_set ... 89!> \param pos ... 90!> \param vel ... 91!> \param r_shake ... 92!> \param v_shake ... 93!> \param dt ... 94!> \param ishake ... 95!> \param max_sigma ... 96!> \par History 97!> none 98! ************************************************************************************************** 99 SUBROUTINE shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, & 100 v_shake, dt, ishake, max_sigma) 101 102 TYPE(molecule_type), POINTER :: molecule 103 TYPE(particle_type), POINTER :: particle_set(:) 104 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 105 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake 106 REAL(kind=dp), INTENT(in) :: dt 107 INTEGER, INTENT(IN) :: ishake 108 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 109 110 INTEGER :: first_atom, ng3x3 111 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 112 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 113 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 114 TYPE(molecule_kind_type), POINTER :: molecule_kind 115 116 NULLIFY (fixd_list) 117 molecule_kind => molecule%molecule_kind 118 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, & 119 g3x3_list=g3x3_list, fixd_list=fixd_list) 120 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3) 121 ! Real Shake 122 CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 123 particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma) 124 125 END SUBROUTINE shake_roll_3x3_int 126 127! ************************************************************************************************** 128!> \brief ... 129!> \param molecule ... 130!> \param particle_set ... 131!> \param vel ... 132!> \param r_rattle ... 133!> \param dt ... 134!> \param veps ... 135!> \par History 136!> none 137! ************************************************************************************************** 138 SUBROUTINE rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, veps) 139 140 TYPE(molecule_type), POINTER :: molecule 141 TYPE(particle_type), POINTER :: particle_set(:) 142 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 143 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 144 REAL(kind=dp), INTENT(in) :: dt 145 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 146 147 INTEGER :: first_atom 148 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 149 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 150 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 151 TYPE(molecule_kind_type), POINTER :: molecule_kind 152 153 NULLIFY (fixd_list) 154 molecule_kind => molecule%molecule_kind 155 CALL get_molecule_kind(molecule_kind, & 156 g3x3_list=g3x3_list, fixd_list=fixd_list) 157 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3) 158 ! Real Rattle 159 CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 160 particle_set, vel, r_rattle, dt, veps) 161 162 END SUBROUTINE rattle_roll_3x3_int 163 164! ************************************************************************************************** 165!> \brief ... 166!> \param molecule ... 167!> \param particle_set ... 168!> \param vel ... 169!> \param dt ... 170!> \par History 171!> none 172! ************************************************************************************************** 173 SUBROUTINE rattle_3x3_int(molecule, particle_set, vel, dt) 174 175 TYPE(molecule_type), POINTER :: molecule 176 TYPE(particle_type), POINTER :: particle_set(:) 177 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 178 REAL(kind=dp), INTENT(in) :: dt 179 180 INTEGER :: first_atom 181 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 182 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 183 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 184 TYPE(molecule_kind_type), POINTER :: molecule_kind 185 186 NULLIFY (fixd_list) 187 molecule_kind => molecule%molecule_kind 188 CALL get_molecule_kind(molecule_kind, & 189 g3x3_list=g3x3_list, fixd_list=fixd_list) 190 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3) 191 ! Real Rattle 192 CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 193 particle_set, vel, dt) 194 195 END SUBROUTINE rattle_3x3_int 196 197! ************************************************************************************************** 198!> \brief ... 199!> \param gci ... 200!> \param particle_set ... 201!> \param pos ... 202!> \param vel ... 203!> \param dt ... 204!> \param ishake ... 205!> \param max_sigma ... 206!> \par History 207!> none 208! ************************************************************************************************** 209 SUBROUTINE shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake, & 210 max_sigma) 211 212 TYPE(global_constraint_type), POINTER :: gci 213 TYPE(particle_type), POINTER :: particle_set(:) 214 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 215 REAL(kind=dp), INTENT(in) :: dt 216 INTEGER, INTENT(IN) :: ishake 217 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 218 219 INTEGER :: first_atom, ng3x3 220 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 221 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 222 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 223 224 first_atom = 1 225 ng3x3 = gci%ng3x3 226 g3x3_list => gci%g3x3_list 227 fixd_list => gci%fixd_list 228 lg3x3 => gci%lg3x3 229 ! Real Shake 230 CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 231 particle_set, pos, vel, dt, ishake, max_sigma) 232 233 END SUBROUTINE shake_3x3_ext 234 235! ************************************************************************************************** 236!> \brief ... 237!> \param gci ... 238!> \param particle_set ... 239!> \param pos ... 240!> \param vel ... 241!> \param r_shake ... 242!> \param v_shake ... 243!> \param dt ... 244!> \param ishake ... 245!> \param max_sigma ... 246!> \par History 247!> none 248! ************************************************************************************************** 249 SUBROUTINE shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, & 250 v_shake, dt, ishake, max_sigma) 251 252 TYPE(global_constraint_type), POINTER :: gci 253 TYPE(particle_type), POINTER :: particle_set(:) 254 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 255 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake 256 REAL(kind=dp), INTENT(in) :: dt 257 INTEGER, INTENT(IN) :: ishake 258 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 259 260 INTEGER :: first_atom, ng3x3 261 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 262 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 263 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 264 265 first_atom = 1 266 ng3x3 = gci%ng3x3 267 g3x3_list => gci%g3x3_list 268 fixd_list => gci%fixd_list 269 lg3x3 => gci%lg3x3 270 ! Real Shake 271 CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 272 particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma) 273 274 END SUBROUTINE shake_roll_3x3_ext 275 276! ************************************************************************************************** 277!> \brief ... 278!> \param gci ... 279!> \param particle_set ... 280!> \param vel ... 281!> \param r_rattle ... 282!> \param dt ... 283!> \param veps ... 284!> \par History 285!> none 286! ************************************************************************************************** 287 SUBROUTINE rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, veps) 288 289 TYPE(global_constraint_type), POINTER :: gci 290 TYPE(particle_type), POINTER :: particle_set(:) 291 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 292 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 293 REAL(kind=dp), INTENT(in) :: dt 294 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 295 296 INTEGER :: first_atom 297 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 298 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 299 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 300 301 first_atom = 1 302 g3x3_list => gci%g3x3_list 303 fixd_list => gci%fixd_list 304 lg3x3 => gci%lg3x3 305 ! Real Rattle 306 CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 307 particle_set, vel, r_rattle, dt, veps) 308 309 END SUBROUTINE rattle_roll_3x3_ext 310 311! ************************************************************************************************** 312!> \brief ... 313!> \param gci ... 314!> \param particle_set ... 315!> \param vel ... 316!> \param dt ... 317!> \par History 318!> none 319! ************************************************************************************************** 320 SUBROUTINE rattle_3x3_ext(gci, particle_set, vel, dt) 321 322 TYPE(global_constraint_type), POINTER :: gci 323 TYPE(particle_type), POINTER :: particle_set(:) 324 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 325 REAL(kind=dp), INTENT(in) :: dt 326 327 INTEGER :: first_atom 328 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 329 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 330 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 331 332 first_atom = 1 333 g3x3_list => gci%g3x3_list 334 fixd_list => gci%fixd_list 335 lg3x3 => gci%lg3x3 336 ! Real Rattle 337 CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 338 particle_set, vel, dt) 339 340 END SUBROUTINE rattle_3x3_ext 341 342! ************************************************************************************************** 343!> \brief ... 344!> \param fixd_list ... 345!> \param g3x3_list ... 346!> \param lg3x3 ... 347!> \param first_atom ... 348!> \param ng3x3 ... 349!> \param particle_set ... 350!> \param pos ... 351!> \param vel ... 352!> \param dt ... 353!> \param ishake ... 354!> \param max_sigma ... 355!> \par History 356!> none 357! ************************************************************************************************** 358 SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 359 particle_set, pos, vel, dt, ishake, max_sigma) 360 361 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 362 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 363 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 364 INTEGER, INTENT(IN) :: first_atom, ng3x3 365 TYPE(particle_type), POINTER :: particle_set(:) 366 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 367 REAL(kind=dp), INTENT(in) :: dt 368 INTEGER, INTENT(IN) :: ishake 369 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 370 371 INTEGER :: iconst, index_a, index_b, index_c 372 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, & 373 imass3, sigma 374 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, r0_12, r0_13, r0_23, r1, & 375 r12, r13, r2, r23, r3, v1, v2, v3, vec 376 REAL(KIND=dp), DIMENSION(3, 1) :: bvec 377 REAL(KIND=dp), DIMENSION(3, 3) :: amat, atemp 378 TYPE(atomic_kind_type), POINTER :: atomic_kind 379 380 dtsqby2 = dt*dt*.5_dp 381 idtsq = 1.0_dp/dt/dt 382 dtby2 = dt*.5_dp 383 DO iconst = 1, ng3x3 384 IF (g3x3_list(iconst)%restraint%active) CYCLE 385 index_a = g3x3_list(iconst)%a + first_atom - 1 386 index_b = g3x3_list(iconst)%b + first_atom - 1 387 index_c = g3x3_list(iconst)%c + first_atom - 1 388 IF (ishake == 1) THEN 389 r0_12(:) = pos(:, index_a) - pos(:, index_b) 390 r0_13(:) = pos(:, index_a) - pos(:, index_c) 391 r0_23(:) = pos(:, index_b) - pos(:, index_c) 392 atomic_kind => particle_set(index_a)%atomic_kind 393 imass1 = 1.0_dp/atomic_kind%mass 394 atomic_kind => particle_set(index_b)%atomic_kind 395 imass2 = 1.0_dp/atomic_kind%mass 396 atomic_kind => particle_set(index_c)%atomic_kind 397 imass3 = 1.0_dp/atomic_kind%mass 398 lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - & 399 lg3x3(iconst)%rb_old) 400 lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - & 401 lg3x3(iconst)%rc_old) 402 lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - & 403 lg3x3(iconst)%rc_old) 404 ! Check for fixed atom constraints 405 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, & 406 index_a, index_b, index_c, fixd_list, lg3x3(iconst)) 407 ! construct matrix 408 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, lg3x3(iconst)%fa) 409 amat(1, 2) = imass1*DOTPROD_3D(r0_12, lg3x3(iconst)%fb) 410 amat(1, 3) = -imass2*DOTPROD_3D(r0_12, lg3x3(iconst)%fc) 411 amat(2, 1) = imass1*DOTPROD_3D(r0_13, lg3x3(iconst)%fa) 412 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, lg3x3(iconst)%fb) 413 amat(2, 3) = imass3*DOTPROD_3D(r0_13, lg3x3(iconst)%fc) 414 amat(3, 1) = -imass2*DOTPROD_3D(r0_23, lg3x3(iconst)%fa) 415 amat(3, 2) = imass3*DOTPROD_3D(r0_23, lg3x3(iconst)%fb) 416 amat(3, 3) = (imass3 + imass2)*DOTPROD_3D(r0_23, lg3x3(iconst)%fc) 417 ! Store values 418 lg3x3(iconst)%r0_12 = r0_12 419 lg3x3(iconst)%r0_13 = r0_13 420 lg3x3(iconst)%r0_23 = r0_23 421 lg3x3(iconst)%amat = amat 422 lg3x3(iconst)%lambda_old = 0.0_dp 423 lg3x3(iconst)%del_lambda = 0.0_dp 424 ELSE 425 ! Retrieve values 426 r0_12 = lg3x3(iconst)%r0_12 427 r0_13 = lg3x3(iconst)%r0_13 428 r0_23 = lg3x3(iconst)%r0_23 429 amat = lg3x3(iconst)%amat 430 imass1 = lg3x3(iconst)%imass1 431 imass2 = lg3x3(iconst)%imass2 432 imass3 = lg3x3(iconst)%imass3 433 END IF 434 ! Iterate until convergence 435 vec = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*(imass1 + imass2) + & 436 lg3x3(iconst)%lambda(2)*imass1*lg3x3(iconst)%fb - & 437 lg3x3(iconst)%lambda(3)*imass2*lg3x3(iconst)%fc 438 bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab & 439 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12) 440 vec = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass1 + & 441 lg3x3(iconst)%lambda(2)*(imass1 + imass3)*lg3x3(iconst)%fb + & 442 lg3x3(iconst)%lambda(3)*imass3*lg3x3(iconst)%fc 443 bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac & 444 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13) 445 vec = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass2 + & 446 lg3x3(iconst)%lambda(2)*imass3*lg3x3(iconst)%fb + & 447 lg3x3(iconst)%lambda(3)*(imass2 + imass3)*lg3x3(iconst)%fc 448 bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc & 449 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23) 450 bvec = bvec*idtsq 451 ! get lambda 452 atemp = amat 453 CALL solve_system(atemp, 3, bvec) 454 lg3x3(iconst)%lambda(:) = bvec(:, 1) 455 lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - & 456 lg3x3(iconst)%lambda_old(:) 457 lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:) 458 fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + & 459 lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb 460 fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + & 461 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc 462 fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - & 463 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc 464 r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:) 465 r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:) 466 r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:) 467 v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:) 468 v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:) 469 v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:) 470 r12 = r1 - r2 471 r13 = r1 - r3 472 r23 = r2 - r3 473 474 ! compute the tolerance: 475 sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* & 476 g3x3_list(iconst)%dab 477 max_sigma = MAX(max_sigma, ABS(sigma)) 478 sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* & 479 g3x3_list(iconst)%dac 480 max_sigma = MAX(max_sigma, ABS(sigma)) 481 sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* & 482 g3x3_list(iconst)%dbc 483 max_sigma = MAX(max_sigma, ABS(sigma)) 484 ! update positions with full multiplier 485 pos(:, index_a) = r1(:) 486 pos(:, index_b) = r2(:) 487 pos(:, index_c) = r3(:) 488 489 ! update velocites with full multiplier 490 vel(:, index_a) = v1(:) 491 vel(:, index_b) = v2(:) 492 vel(:, index_c) = v3(:) 493 END DO 494 495 END SUBROUTINE shake_3x3_low 496 497! ************************************************************************************************** 498!> \brief ... 499!> \param fixd_list ... 500!> \param g3x3_list ... 501!> \param lg3x3 ... 502!> \param first_atom ... 503!> \param ng3x3 ... 504!> \param particle_set ... 505!> \param pos ... 506!> \param vel ... 507!> \param r_shake ... 508!> \param v_shake ... 509!> \param dt ... 510!> \param ishake ... 511!> \param max_sigma ... 512!> \par History 513!> none 514! ************************************************************************************************** 515 SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, & 516 particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma) 517 518 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 519 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 520 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 521 INTEGER, INTENT(IN) :: first_atom, ng3x3 522 TYPE(particle_type), POINTER :: particle_set(:) 523 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 524 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake 525 REAL(kind=dp), INTENT(in) :: dt 526 INTEGER, INTENT(IN) :: ishake 527 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 528 529 INTEGER :: iconst, index_a, index_b, index_c 530 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, & 531 imass3, sigma 532 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, fc1, fc2, & 533 fc3, r0_12, r0_13, r0_23, r1, r12, & 534 r13, r2, r23, r3, v1, v2, v3, vec 535 REAL(KIND=dp), DIMENSION(3, 1) :: bvec 536 REAL(KIND=dp), DIMENSION(3, 3) :: amat, atemp 537 TYPE(atomic_kind_type), POINTER :: atomic_kind 538 539 dtsqby2 = dt*dt*.5_dp 540 idtsq = 1.0_dp/dt/dt 541 dtby2 = dt*.5_dp 542 DO iconst = 1, ng3x3 543 IF (g3x3_list(iconst)%restraint%active) CYCLE 544 index_a = g3x3_list(iconst)%a + first_atom - 1 545 index_b = g3x3_list(iconst)%b + first_atom - 1 546 index_c = g3x3_list(iconst)%c + first_atom - 1 547 IF (ishake == 1) THEN 548 r0_12(:) = pos(:, index_a) - pos(:, index_b) 549 r0_13(:) = pos(:, index_a) - pos(:, index_c) 550 r0_23(:) = pos(:, index_b) - pos(:, index_c) 551 atomic_kind => particle_set(index_a)%atomic_kind 552 imass1 = 1.0_dp/atomic_kind%mass 553 atomic_kind => particle_set(index_b)%atomic_kind 554 imass2 = 1.0_dp/atomic_kind%mass 555 atomic_kind => particle_set(index_c)%atomic_kind 556 imass3 = 1.0_dp/atomic_kind%mass 557 lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - & 558 lg3x3(iconst)%rb_old) 559 lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - & 560 lg3x3(iconst)%rc_old) 561 lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - & 562 lg3x3(iconst)%rc_old) 563 ! Check for fixed atom constraints 564 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, & 565 index_a, index_b, index_c, fixd_list, lg3x3(iconst)) 566 ! rotate fconst: 567 CALL matvec_3x3(f_roll1, r_shake, lg3x3(iconst)%fa) 568 CALL matvec_3x3(f_roll2, r_shake, lg3x3(iconst)%fb) 569 CALL matvec_3x3(f_roll3, r_shake, lg3x3(iconst)%fc) 570 ! construct matrix 571 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, f_roll1) 572 amat(1, 2) = imass1*DOTPROD_3D(r0_12, f_roll2) 573 amat(1, 3) = -imass2*DOTPROD_3D(r0_12, f_roll3) 574 amat(2, 1) = imass1*DOTPROD_3D(r0_13, f_roll1) 575 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, f_roll2) 576 amat(2, 3) = imass3*DOTPROD_3D(r0_13, f_roll3) 577 amat(3, 1) = -imass2*DOTPROD_3D(r0_23, f_roll1) 578 amat(3, 2) = imass3*DOTPROD_3D(r0_23, f_roll2) 579 amat(3, 3) = (imass3 + imass2)*DOTPROD_3D(r0_23, f_roll3) 580 ! Store values 581 lg3x3(iconst)%r0_12 = r0_12 582 lg3x3(iconst)%r0_13 = r0_13 583 lg3x3(iconst)%r0_23 = r0_23 584 lg3x3(iconst)%f_roll1 = f_roll1 585 lg3x3(iconst)%f_roll2 = f_roll2 586 lg3x3(iconst)%f_roll3 = f_roll3 587 lg3x3(iconst)%amat = amat 588 lg3x3(iconst)%lambda_old = 0.0_dp 589 lg3x3(iconst)%del_lambda = 0.0_dp 590 ELSE 591 ! Retrieve values 592 r0_12 = lg3x3(iconst)%r0_12 593 r0_13 = lg3x3(iconst)%r0_13 594 r0_23 = lg3x3(iconst)%r0_23 595 f_roll1 = lg3x3(iconst)%f_roll1 596 f_roll2 = lg3x3(iconst)%f_roll2 597 f_roll3 = lg3x3(iconst)%f_roll3 598 amat = lg3x3(iconst)%amat 599 imass1 = lg3x3(iconst)%imass1 600 imass2 = lg3x3(iconst)%imass2 601 imass3 = lg3x3(iconst)%imass3 602 END IF 603 ! Iterate until convergence 604 vec = lg3x3(iconst)%lambda(1)*f_roll1*(imass1 + imass2) + & 605 lg3x3(iconst)%lambda(2)*imass1*f_roll2 - & 606 lg3x3(iconst)%lambda(3)*imass2*f_roll3 607 bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab & 608 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12) 609 vec = lg3x3(iconst)%lambda(1)*f_roll1*imass1 + & 610 lg3x3(iconst)%lambda(2)*(imass1 + imass3)*f_roll2 + & 611 lg3x3(iconst)%lambda(3)*imass3*f_roll3 612 bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac & 613 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13) 614 vec = -lg3x3(iconst)%lambda(1)*f_roll1*imass2 + & 615 lg3x3(iconst)%lambda(2)*imass3*f_roll2 + & 616 lg3x3(iconst)%lambda(3)*(imass2 + imass3)*f_roll3 617 bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc & 618 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23) 619 bvec = bvec*idtsq 620 621 ! get lambda 622 atemp = amat 623 CALL solve_system(atemp, 3, bvec) 624 lg3x3(iconst)%lambda(:) = bvec(:, 1) 625 lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - & 626 lg3x3(iconst)%lambda_old(:) 627 lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:) 628 629 fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + & 630 lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb 631 fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + & 632 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc 633 fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - & 634 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc 635 CALL MATVEC_3x3(vec, r_shake, fc1) 636 r1(:) = pos(:, index_a) + imass1*dtsqby2*vec 637 CALL MATVEC_3x3(vec, r_shake, fc2) 638 r2(:) = pos(:, index_b) + imass2*dtsqby2*vec 639 CALL MATVEC_3x3(vec, r_shake, fc3) 640 r3(:) = pos(:, index_c) + imass3*dtsqby2*vec 641 CALL MATVEC_3x3(vec, v_shake, fc1) 642 v1(:) = vel(:, index_a) + imass1*dtby2*vec 643 CALL MATVEC_3x3(vec, v_shake, fc2) 644 v2(:) = vel(:, index_b) + imass2*dtby2*vec 645 CALL MATVEC_3x3(vec, v_shake, fc3) 646 v3(:) = vel(:, index_c) + imass3*dtby2*vec 647 r12 = r1 - r2 648 r13 = r1 - r3 649 r23 = r2 - r3 650 651 ! compute the tolerance: 652 sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* & 653 g3x3_list(iconst)%dab 654 max_sigma = MAX(max_sigma, ABS(sigma)) 655 sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* & 656 g3x3_list(iconst)%dac 657 max_sigma = MAX(max_sigma, ABS(sigma)) 658 sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* & 659 g3x3_list(iconst)%dbc 660 max_sigma = MAX(max_sigma, ABS(sigma)) 661 662 ! update positions with full multiplier 663 pos(:, index_a) = r1(:) 664 pos(:, index_b) = r2(:) 665 pos(:, index_c) = r3(:) 666 667 ! update velocites with full multiplier 668 vel(:, index_a) = v1(:) 669 vel(:, index_b) = v2(:) 670 vel(:, index_c) = v3(:) 671 END DO 672 END SUBROUTINE shake_roll_3x3_low 673 674! ************************************************************************************************** 675!> \brief ... 676!> \param fixd_list ... 677!> \param g3x3_list ... 678!> \param lg3x3 ... 679!> \param first_atom ... 680!> \param particle_set ... 681!> \param vel ... 682!> \param r_rattle ... 683!> \param dt ... 684!> \param veps ... 685!> \par History 686!> none 687! ************************************************************************************************** 688 SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 689 particle_set, vel, r_rattle, dt, veps) 690 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 691 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 692 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 693 INTEGER, INTENT(IN) :: first_atom 694 TYPE(particle_type), POINTER :: particle_set(:) 695 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 696 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 697 REAL(kind=dp), INTENT(in) :: dt 698 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 699 700 INTEGER :: iconst, index_a, index_b, index_c 701 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, mass 702 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, fc1, fc2, & 703 fc3, lambda, r12, r13, r23, v12, v13, & 704 v23, vec 705 REAL(KIND=dp), DIMENSION(3, 1) :: bvec 706 REAL(KIND=dp), DIMENSION(3, 3) :: amat 707 TYPE(atomic_kind_type), POINTER :: atomic_kind 708 709 idt = 1.0_dp/dt 710 dtby2 = dt*.5_dp 711 DO iconst = 1, SIZE(g3x3_list) 712 IF (g3x3_list(iconst)%restraint%active) CYCLE 713 index_a = g3x3_list(iconst)%a + first_atom - 1 714 index_b = g3x3_list(iconst)%b + first_atom - 1 715 index_c = g3x3_list(iconst)%c + first_atom - 1 716 v12(:) = vel(:, index_a) - vel(:, index_b) 717 v13(:) = vel(:, index_a) - vel(:, index_c) 718 v23(:) = vel(:, index_b) - vel(:, index_c) 719 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:) 720 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:) 721 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:) 722 atomic_kind => particle_set(index_a)%atomic_kind 723 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 724 imass1 = 1.0_dp/mass 725 atomic_kind => particle_set(index_b)%atomic_kind 726 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 727 imass2 = 1.0_dp/mass 728 atomic_kind => particle_set(index_c)%atomic_kind 729 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 730 imass3 = 1.0_dp/mass 731 lg3x3(iconst)%fa = -2.0_dp*r12 732 lg3x3(iconst)%fb = -2.0_dp*r13 733 lg3x3(iconst)%fc = -2.0_dp*r23 734 ! Check for fixed atom constraints 735 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, & 736 index_a, index_b, index_c, fixd_list, lg3x3(iconst)) 737 ! roll the fc 738 CALL MATVEC_3x3(f_roll1, r_rattle, lg3x3(iconst)%fa) 739 CALL MATVEC_3x3(f_roll2, r_rattle, lg3x3(iconst)%fb) 740 CALL MATVEC_3x3(f_roll3, r_rattle, lg3x3(iconst)%fc) 741 742 ! construct matrix 743 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, f_roll1) 744 amat(1, 2) = imass1*DOTPROD_3D(r12, f_roll2) 745 amat(1, 3) = -imass2*DOTPROD_3D(r12, f_roll3) 746 amat(2, 1) = imass1*DOTPROD_3D(r13, f_roll1) 747 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, f_roll2) 748 amat(2, 3) = imass3*DOTPROD_3D(r13, f_roll3) 749 amat(3, 1) = -imass2*DOTPROD_3D(r23, f_roll1) 750 amat(3, 2) = imass3*DOTPROD_3D(r23, f_roll2) 751 amat(3, 3) = (imass2 + imass3)*DOTPROD_3D(r23, f_roll3) 752 753 ! construct solution vector 754 CALL matvec_3x3(vec, veps, r12) 755 bvec(1, 1) = DOTPROD_3D(r12, v12 + vec) 756 CALL matvec_3x3(vec, veps, r13) 757 bvec(2, 1) = DOTPROD_3D(r13, v13 + vec) 758 CALL matvec_3x3(vec, veps, r23) 759 bvec(3, 1) = DOTPROD_3D(r23, v23 + vec) 760 bvec = -bvec*2.0_dp*idt 761 762 ! get lambda 763 CALL solve_system(amat, 3, bvec) 764 lambda(:) = bvec(:, 1) 765 lg3x3(iconst)%lambda(:) = lambda 766 767 fc1 = lambda(1)*f_roll1 + & 768 lambda(2)*f_roll2 769 fc2 = -lambda(1)*f_roll1 + & 770 lambda(3)*f_roll3 771 fc3 = -lambda(2)*f_roll2 - & 772 lambda(3)*f_roll3 773 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:) 774 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:) 775 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:) 776 END DO 777 END SUBROUTINE rattle_roll_3x3_low 778 779! ************************************************************************************************** 780!> \brief ... 781!> \param fixd_list ... 782!> \param g3x3_list ... 783!> \param lg3x3 ... 784!> \param first_atom ... 785!> \param particle_set ... 786!> \param vel ... 787!> \param dt ... 788!> \par History 789!> none 790! ************************************************************************************************** 791 SUBROUTINE rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, & 792 particle_set, vel, dt) 793 794 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 795 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 796 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:) 797 INTEGER, INTENT(IN) :: first_atom 798 TYPE(particle_type), POINTER :: particle_set(:) 799 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 800 REAL(kind=dp), INTENT(in) :: dt 801 802 INTEGER :: iconst, index_a, index_b, index_c 803 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, mass 804 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, r12, r13, r23, v12, v13, & 805 v23 806 REAL(KIND=dp), DIMENSION(3, 1) :: bvec 807 REAL(KIND=dp), DIMENSION(3, 3) :: amat 808 TYPE(atomic_kind_type), POINTER :: atomic_kind 809 810 idt = 1.0_dp/dt 811 dtby2 = dt*.5_dp 812 DO iconst = 1, SIZE(g3x3_list) 813 IF (g3x3_list(iconst)%restraint%active) CYCLE 814 index_a = g3x3_list(iconst)%a + first_atom - 1 815 index_b = g3x3_list(iconst)%b + first_atom - 1 816 index_c = g3x3_list(iconst)%c + first_atom - 1 817 v12(:) = vel(:, index_a) - vel(:, index_b) 818 v13(:) = vel(:, index_a) - vel(:, index_c) 819 v23(:) = vel(:, index_b) - vel(:, index_c) 820 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:) 821 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:) 822 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:) 823 atomic_kind => particle_set(index_a)%atomic_kind 824 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 825 imass1 = 1.0_dp/mass 826 atomic_kind => particle_set(index_b)%atomic_kind 827 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 828 imass2 = 1.0_dp/mass 829 atomic_kind => particle_set(index_c)%atomic_kind 830 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 831 imass3 = 1.0_dp/mass 832 lg3x3(iconst)%fa = -2.0_dp*r12 833 lg3x3(iconst)%fb = -2.0_dp*r13 834 lg3x3(iconst)%fc = -2.0_dp*r23 835 ! Check for fixed atom constraints 836 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, & 837 index_a, index_b, index_c, fixd_list, lg3x3(iconst)) 838 ! construct matrix 839 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, lg3x3(iconst)%fa) 840 amat(1, 2) = imass1*DOTPROD_3D(r12, lg3x3(iconst)%fb) 841 amat(1, 3) = -imass2*DOTPROD_3D(r12, lg3x3(iconst)%fc) 842 amat(2, 1) = imass1*DOTPROD_3D(r13, lg3x3(iconst)%fa) 843 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, lg3x3(iconst)%fb) 844 amat(2, 3) = imass3*DOTPROD_3D(r13, lg3x3(iconst)%fc) 845 amat(3, 1) = -imass2*DOTPROD_3D(r23, lg3x3(iconst)%fa) 846 amat(3, 2) = imass3*DOTPROD_3D(r23, lg3x3(iconst)%fb) 847 amat(3, 3) = (imass2 + imass3)*DOTPROD_3D(r23, lg3x3(iconst)%fc) 848 849 ! construct solution vector 850 bvec(1, 1) = DOTPROD_3D(r12, v12) 851 bvec(2, 1) = DOTPROD_3D(r13, v13) 852 bvec(3, 1) = DOTPROD_3D(r23, v23) 853 bvec = -bvec*2.0_dp*idt 854 855 ! get lambda 856 CALL solve_system(amat, 3, bvec) 857 lg3x3(iconst)%lambda(:) = bvec(:, 1) 858 859 fc1 = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + & 860 lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb 861 fc2 = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + & 862 lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc 863 fc3 = -lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb - & 864 lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc 865 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:) 866 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:) 867 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:) 868 END DO 869 END SUBROUTINE rattle_3x3_low 870 871END MODULE constraint_3x3 872