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_4x6 11 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE constraint_fxd, ONLY: check_fixed_atom_cns_g4x6 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 g4x6_constraint_type,& 21 get_molecule_kind,& 22 molecule_kind_type 23 USE molecule_types, ONLY: get_molecule,& 24 global_constraint_type,& 25 local_g4x6_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_4x6_int, & 34 rattle_4x6_int, & 35 shake_roll_4x6_int, & 36 rattle_roll_4x6_int, & 37 shake_4x6_ext, & 38 rattle_4x6_ext, & 39 shake_roll_4x6_ext, & 40 rattle_roll_4x6_ext 41 42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_4x6' 43 44CONTAINS 45 46! ************************************************************************************************** 47!> \brief Intramolecular shake_4x6 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_4x6_int(molecule, particle_set, pos, vel, dt, ishake, & 59 max_sigma) 60 TYPE(molecule_type), POINTER :: molecule 61 TYPE(particle_type), POINTER :: particle_set(:) 62 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 63 REAL(kind=dp), INTENT(in) :: dt 64 INTEGER, INTENT(IN) :: ishake 65 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 66 67 INTEGER :: first_atom, ng4x6 68 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 69 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 70 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 71 TYPE(molecule_kind_type), POINTER :: molecule_kind 72 73 NULLIFY (fixd_list) 74 molecule_kind => molecule%molecule_kind 75 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, & 76 fixd_list=fixd_list, g4x6_list=g4x6_list) 77 CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6) 78 ! Real Shake 79 CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 80 particle_set, pos, vel, dt, ishake, max_sigma) 81 82 END SUBROUTINE shake_4x6_int 83 84! ************************************************************************************************** 85!> \brief Intramolecular shake_4x6_roll 86!> \param molecule ... 87!> \param particle_set ... 88!> \param pos ... 89!> \param vel ... 90!> \param r_shake ... 91!> \param dt ... 92!> \param ishake ... 93!> \param max_sigma ... 94!> \par History 95!> none 96! ************************************************************************************************** 97 SUBROUTINE shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, & 98 dt, ishake, max_sigma) 99 TYPE(molecule_type), POINTER :: molecule 100 TYPE(particle_type), POINTER :: particle_set(:) 101 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 102 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake 103 REAL(kind=dp), INTENT(in) :: dt 104 INTEGER, INTENT(IN) :: ishake 105 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 106 107 INTEGER :: first_atom, ng4x6 108 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 109 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 110 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 111 TYPE(molecule_kind_type), POINTER :: molecule_kind 112 113 NULLIFY (fixd_list) 114 molecule_kind => molecule%molecule_kind 115 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, & 116 fixd_list=fixd_list, g4x6_list=g4x6_list) 117 CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6) 118 ! Real Shake 119 CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 120 particle_set, pos, vel, r_shake, dt, ishake, max_sigma) 121 122 END SUBROUTINE shake_roll_4x6_int 123 124! ************************************************************************************************** 125!> \brief Intramolecular rattle_4x6 126!> \param molecule ... 127!> \param particle_set ... 128!> \param vel ... 129!> \param dt ... 130!> \par History 131!> none 132! ************************************************************************************************** 133 SUBROUTINE rattle_4x6_int(molecule, particle_set, vel, dt) 134 TYPE(molecule_type), POINTER :: molecule 135 TYPE(particle_type), POINTER :: particle_set(:) 136 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 137 REAL(kind=dp), INTENT(in) :: dt 138 139 INTEGER :: first_atom, ng4x6 140 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 141 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 142 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 143 TYPE(molecule_kind_type), POINTER :: molecule_kind 144 145 NULLIFY (fixd_list) 146 molecule_kind => molecule%molecule_kind 147 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, & 148 fixd_list=fixd_list, g4x6_list=g4x6_list) 149 CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6) 150 ! Real Rattle 151 CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 152 particle_set, vel, dt) 153 154 END SUBROUTINE rattle_4x6_int 155 156! ************************************************************************************************** 157!> \brief Intramolecular rattle_4x6_roll 158!> \param molecule ... 159!> \param particle_set ... 160!> \param vel ... 161!> \param r_rattle ... 162!> \param dt ... 163!> \param veps ... 164!> \par History 165!> none 166! ************************************************************************************************** 167 SUBROUTINE rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, veps) 168 TYPE(molecule_type), POINTER :: molecule 169 TYPE(particle_type), POINTER :: particle_set(:) 170 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 171 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 172 REAL(kind=dp), INTENT(in) :: dt 173 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 174 175 INTEGER :: first_atom, ng4x6 176 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 177 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 178 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 179 TYPE(molecule_kind_type), POINTER :: molecule_kind 180 181 NULLIFY (fixd_list) 182 molecule_kind => molecule%molecule_kind 183 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, & 184 fixd_list=fixd_list, g4x6_list=g4x6_list) 185 CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6) 186 ! Real Rattle 187 CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 188 particle_set, vel, r_rattle, dt, veps) 189 190 END SUBROUTINE rattle_roll_4x6_int 191 192! ************************************************************************************************** 193!> \brief Intramolecular shake_4x6 194!> \param gci ... 195!> \param particle_set ... 196!> \param pos ... 197!> \param vel ... 198!> \param dt ... 199!> \param ishake ... 200!> \param max_sigma ... 201!> \par History 202!> none 203! ************************************************************************************************** 204 SUBROUTINE shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake, & 205 max_sigma) 206 207 TYPE(global_constraint_type), POINTER :: gci 208 TYPE(particle_type), POINTER :: particle_set(:) 209 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 210 REAL(kind=dp), INTENT(in) :: dt 211 INTEGER, INTENT(IN) :: ishake 212 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 213 214 INTEGER :: first_atom, ng4x6 215 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 216 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 217 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 218 219 first_atom = 1 220 ng4x6 = gci%ng4x6 221 fixd_list => gci%fixd_list 222 g4x6_list => gci%g4x6_list 223 lg4x6 => gci%lg4x6 224 ! Real Shake 225 CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 226 particle_set, pos, vel, dt, ishake, max_sigma) 227 228 END SUBROUTINE shake_4x6_ext 229 230! ************************************************************************************************** 231!> \brief Intramolecular shake_4x6_roll 232!> \param gci ... 233!> \param particle_set ... 234!> \param pos ... 235!> \param vel ... 236!> \param r_shake ... 237!> \param dt ... 238!> \param ishake ... 239!> \param max_sigma ... 240!> \par History 241!> none 242! ************************************************************************************************** 243 SUBROUTINE shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, & 244 dt, ishake, max_sigma) 245 246 TYPE(global_constraint_type), POINTER :: gci 247 TYPE(particle_type), POINTER :: particle_set(:) 248 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 249 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake 250 REAL(kind=dp), INTENT(in) :: dt 251 INTEGER, INTENT(IN) :: ishake 252 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 253 254 INTEGER :: first_atom, ng4x6 255 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 256 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 257 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 258 259 first_atom = 1 260 ng4x6 = gci%ng4x6 261 fixd_list => gci%fixd_list 262 g4x6_list => gci%g4x6_list 263 lg4x6 => gci%lg4x6 264 ! Real Shake 265 CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 266 particle_set, pos, vel, r_shake, dt, ishake, max_sigma) 267 268 END SUBROUTINE shake_roll_4x6_ext 269 270! ************************************************************************************************** 271!> \brief Intramolecular rattle_4x6 272!> \param gci ... 273!> \param particle_set ... 274!> \param vel ... 275!> \param dt ... 276!> \par History 277!> none 278! ************************************************************************************************** 279 SUBROUTINE rattle_4x6_ext(gci, particle_set, vel, dt) 280 281 TYPE(global_constraint_type), POINTER :: gci 282 TYPE(particle_type), POINTER :: particle_set(:) 283 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 284 REAL(kind=dp), INTENT(in) :: dt 285 286 INTEGER :: first_atom, ng4x6 287 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 288 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 289 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 290 291 first_atom = 1 292 ng4x6 = gci%ng4x6 293 fixd_list => gci%fixd_list 294 g4x6_list => gci%g4x6_list 295 lg4x6 => gci%lg4x6 296 ! Real Rattle 297 CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 298 particle_set, vel, dt) 299 300 END SUBROUTINE rattle_4x6_ext 301 302! ************************************************************************************************** 303!> \brief Intramolecular rattle_4x6_roll 304!> \param gci ... 305!> \param particle_set ... 306!> \param vel ... 307!> \param r_rattle ... 308!> \param dt ... 309!> \param veps ... 310!> \par History 311!> none 312! ************************************************************************************************** 313 SUBROUTINE rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, veps) 314 315 TYPE(global_constraint_type), POINTER :: gci 316 TYPE(particle_type), POINTER :: particle_set(:) 317 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 318 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 319 REAL(kind=dp), INTENT(in) :: dt 320 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 321 322 INTEGER :: first_atom, ng4x6 323 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 324 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 325 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 326 327 first_atom = 1 328 ng4x6 = gci%ng4x6 329 fixd_list => gci%fixd_list 330 g4x6_list => gci%g4x6_list 331 lg4x6 => gci%lg4x6 332 ! Real Rattle 333 CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 334 particle_set, vel, r_rattle, dt, veps) 335 336 END SUBROUTINE rattle_roll_4x6_ext 337 338! ************************************************************************************************** 339!> \brief ... 340!> \param fixd_list ... 341!> \param g4x6_list ... 342!> \param lg4x6 ... 343!> \param ng4x6 ... 344!> \param first_atom ... 345!> \param particle_set ... 346!> \param pos ... 347!> \param vel ... 348!> \param dt ... 349!> \param ishake ... 350!> \param max_sigma ... 351!> \par History 352!> none 353! ************************************************************************************************** 354 SUBROUTINE shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 355 particle_set, pos, vel, dt, ishake, max_sigma) 356 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 357 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 358 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 359 INTEGER, INTENT(IN) :: ng4x6, first_atom 360 TYPE(particle_type), POINTER :: particle_set(:) 361 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 362 REAL(kind=dp), INTENT(in) :: dt 363 INTEGER, INTENT(IN) :: ishake 364 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 365 366 INTEGER :: iconst, index_a, index_b, index_c, & 367 index_d 368 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, & 369 imass3, imass4, sigma 370 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r0_12, r0_13, r0_14, & 371 r0_23, r0_24, r0_34, r1, r12, r13, & 372 r14, r2, r23, r24, r3, r34, r4, v1, & 373 v2, v3, v4, vec 374 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 375 REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp 376 TYPE(atomic_kind_type), POINTER :: atomic_kind 377 378 dtsqby2 = dt*dt*.5_dp 379 idtsq = 1.0_dp/dt/dt 380 dtby2 = dt*.5_dp 381 DO iconst = 1, ng4x6 382 IF (g4x6_list(iconst)%restraint%active) CYCLE 383 index_a = g4x6_list(iconst)%a + first_atom - 1 384 index_b = g4x6_list(iconst)%b + first_atom - 1 385 index_c = g4x6_list(iconst)%c + first_atom - 1 386 index_d = g4x6_list(iconst)%d + first_atom - 1 387 IF (ishake == 1) THEN 388 r0_12(:) = pos(:, index_a) - pos(:, index_b) 389 r0_13(:) = pos(:, index_a) - pos(:, index_c) 390 r0_14(:) = pos(:, index_a) - pos(:, index_d) 391 r0_23(:) = pos(:, index_b) - pos(:, index_c) 392 r0_24(:) = pos(:, index_b) - pos(:, index_d) 393 r0_34(:) = pos(:, index_c) - pos(:, index_d) 394 atomic_kind => particle_set(index_a)%atomic_kind 395 imass1 = 1.0_dp/atomic_kind%mass 396 atomic_kind => particle_set(index_b)%atomic_kind 397 imass2 = 1.0_dp/atomic_kind%mass 398 atomic_kind => particle_set(index_c)%atomic_kind 399 imass3 = 1.0_dp/atomic_kind%mass 400 atomic_kind => particle_set(index_d)%atomic_kind 401 imass4 = 1.0_dp/atomic_kind%mass 402 lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - & 403 lg4x6(iconst)%rb_old) 404 lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - & 405 lg4x6(iconst)%rc_old) 406 lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - & 407 lg4x6(iconst)%rd_old) 408 lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - & 409 lg4x6(iconst)%rc_old) 410 lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - & 411 lg4x6(iconst)%rd_old) 412 lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - & 413 lg4x6(iconst)%rd_old) 414 ! Check for fixed atom constraints 415 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, & 416 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 417 ! construct matrix 418 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, lg4x6(iconst)%fa) 419 amat(1, 2) = imass1*DOTPROD_3D(r0_12, lg4x6(iconst)%fb) 420 amat(1, 3) = imass1*DOTPROD_3D(r0_12, lg4x6(iconst)%fc) 421 amat(1, 4) = -imass2*DOTPROD_3D(r0_12, lg4x6(iconst)%fd) 422 amat(1, 5) = -imass2*DOTPROD_3D(r0_12, lg4x6(iconst)%fe) 423 amat(1, 6) = 0.0_dp 424 425 amat(2, 1) = imass1*DOTPROD_3D(r0_13, lg4x6(iconst)%fa) 426 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, lg4x6(iconst)%fb) 427 amat(2, 3) = imass1*DOTPROD_3D(r0_13, lg4x6(iconst)%fc) 428 amat(2, 4) = imass3*DOTPROD_3D(r0_13, lg4x6(iconst)%fd) 429 amat(2, 5) = 0.0_dp 430 amat(2, 6) = -imass3*DOTPROD_3D(r0_13, lg4x6(iconst)%ff) 431 432 amat(3, 1) = imass1*DOTPROD_3D(r0_14, lg4x6(iconst)%fa) 433 amat(3, 2) = imass1*DOTPROD_3D(r0_14, lg4x6(iconst)%fb) 434 amat(3, 3) = (imass1 + imass4)*DOTPROD_3D(r0_14, lg4x6(iconst)%fc) 435 amat(3, 4) = 0.0_dp 436 amat(3, 5) = imass4*DOTPROD_3D(r0_14, lg4x6(iconst)%fe) 437 amat(3, 6) = imass4*DOTPROD_3D(r0_14, lg4x6(iconst)%ff) 438 439 amat(4, 1) = -imass2*DOTPROD_3D(r0_23, lg4x6(iconst)%fa) 440 amat(4, 2) = imass3*DOTPROD_3D(r0_23, lg4x6(iconst)%fb) 441 amat(4, 3) = 0.0_dp 442 amat(4, 4) = (imass3 + imass2)*DOTPROD_3D(r0_23, lg4x6(iconst)%fd) 443 amat(4, 5) = imass2*DOTPROD_3D(r0_23, lg4x6(iconst)%fe) 444 amat(4, 6) = -imass3*DOTPROD_3D(r0_23, lg4x6(iconst)%ff) 445 446 amat(5, 1) = -imass2*DOTPROD_3D(r0_24, lg4x6(iconst)%fa) 447 amat(5, 2) = 0.0_dp 448 amat(5, 3) = imass4*DOTPROD_3D(r0_24, lg4x6(iconst)%fc) 449 amat(5, 4) = imass2*DOTPROD_3D(r0_24, lg4x6(iconst)%fd) 450 amat(5, 5) = (imass4 + imass2)*DOTPROD_3D(r0_24, lg4x6(iconst)%fe) 451 amat(5, 6) = imass4*DOTPROD_3D(r0_24, lg4x6(iconst)%ff) 452 453 amat(6, 1) = 0.0_dp 454 amat(6, 2) = -imass3*DOTPROD_3D(r0_34, lg4x6(iconst)%fb) 455 amat(6, 3) = imass4*DOTPROD_3D(r0_34, lg4x6(iconst)%fc) 456 amat(6, 4) = -imass3*DOTPROD_3D(r0_34, lg4x6(iconst)%fd) 457 amat(6, 5) = imass4*DOTPROD_3D(r0_34, lg4x6(iconst)%fe) 458 amat(6, 6) = (imass3 + imass4)*DOTPROD_3D(r0_34, lg4x6(iconst)%ff) 459 ! Store values 460 lg4x6(iconst)%r0_12 = r0_12 461 lg4x6(iconst)%r0_13 = r0_13 462 lg4x6(iconst)%r0_14 = r0_14 463 lg4x6(iconst)%r0_23 = r0_23 464 lg4x6(iconst)%r0_24 = r0_24 465 lg4x6(iconst)%r0_34 = r0_34 466 lg4x6(iconst)%amat = amat 467 lg4x6(iconst)%imass1 = imass1 468 lg4x6(iconst)%imass2 = imass2 469 lg4x6(iconst)%imass3 = imass3 470 lg4x6(iconst)%imass4 = imass4 471 lg4x6(iconst)%lambda_old = 0.0_dp 472 lg4x6(iconst)%del_lambda = 0.0_dp 473 ELSE 474 ! Retrieve values 475 r0_12 = lg4x6(iconst)%r0_12 476 r0_13 = lg4x6(iconst)%r0_13 477 r0_14 = lg4x6(iconst)%r0_14 478 r0_23 = lg4x6(iconst)%r0_23 479 r0_24 = lg4x6(iconst)%r0_24 480 r0_34 = lg4x6(iconst)%r0_34 481 amat = lg4x6(iconst)%amat 482 imass1 = lg4x6(iconst)%imass1 483 imass2 = lg4x6(iconst)%imass2 484 imass3 = lg4x6(iconst)%imass3 485 imass4 = lg4x6(iconst)%imass4 486 END IF 487 488 ! Iterate until convergence: 489 vec = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa*(imass1 + imass2) + & 490 lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + & 491 lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc - & 492 lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd - & 493 lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe 494 bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab & 495 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12) 496 vec = lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb*(imass1 + imass3) + & 497 lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + & 498 lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc + & 499 lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd - & 500 lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff 501 bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac & 502 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13) 503 vec = lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc*(imass1 + imass4) + & 504 lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + & 505 lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + & 506 lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe + & 507 lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff 508 bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad & 509 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_14, r0_14) 510 vec = lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd*(imass2 + imass3) - & 511 lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + & 512 lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + & 513 lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe - & 514 lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff 515 bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc & 516 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23) 517 vec = lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe*(imass2 + imass4) - & 518 lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + & 519 lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc + & 520 lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd + & 521 lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff 522 bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd & 523 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_24, r0_24) 524 vec = lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff*(imass3 + imass4) - & 525 lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + & 526 lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc - & 527 lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd + & 528 lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe 529 bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd & 530 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_34, r0_34) 531 bvec = bvec*idtsq 532 533 ! get lambda 534 atemp = amat 535 CALL solve_system(atemp, 6, bvec) 536 lg4x6(iconst)%lambda(:) = bvec(:, 1) 537 lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - & 538 lg4x6(iconst)%lambda_old(:) 539 lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:) 540 541 fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + & 542 lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + & 543 lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc 544 fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + & 545 lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + & 546 lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe 547 fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - & 548 lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + & 549 lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff 550 fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - & 551 lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - & 552 lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff 553 r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:) 554 r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:) 555 r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:) 556 r4(:) = pos(:, index_d) + imass4*dtsqby2*fc4(:) 557 v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:) 558 v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:) 559 v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:) 560 v4(:) = vel(:, index_d) + imass4*dtby2*fc4(:) 561 r12 = r1 - r2 562 r13 = r1 - r3 563 r14 = r1 - r4 564 r23 = r2 - r3 565 r24 = r2 - r4 566 r34 = r3 - r4 567 ! compute the tolerance: 568 sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* & 569 g4x6_list(iconst)%dab 570 max_sigma = MAX(max_sigma, ABS(sigma)) 571 sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* & 572 g4x6_list(iconst)%dac 573 max_sigma = MAX(max_sigma, ABS(sigma)) 574 sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* & 575 g4x6_list(iconst)%dad 576 max_sigma = MAX(max_sigma, ABS(sigma)) 577 sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* & 578 g4x6_list(iconst)%dbc 579 max_sigma = MAX(max_sigma, ABS(sigma)) 580 sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* & 581 g4x6_list(iconst)%dbd 582 max_sigma = MAX(max_sigma, ABS(sigma)) 583 sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* & 584 g4x6_list(iconst)%dcd 585 max_sigma = MAX(max_sigma, ABS(sigma)) 586 587 ! update positions with full multiplier 588 pos(:, index_a) = r1(:) 589 pos(:, index_b) = r2(:) 590 pos(:, index_c) = r3(:) 591 pos(:, index_d) = r4(:) 592 593 ! update velocites with full multiplier 594 vel(:, index_a) = v1(:) 595 vel(:, index_b) = v2(:) 596 vel(:, index_c) = v3(:) 597 vel(:, index_d) = v4(:) 598 END DO 599 END SUBROUTINE shake_4x6_low 600 601! ************************************************************************************************** 602!> \brief ... 603!> \param fixd_list ... 604!> \param g4x6_list ... 605!> \param lg4x6 ... 606!> \param ng4x6 ... 607!> \param first_atom ... 608!> \param particle_set ... 609!> \param pos ... 610!> \param vel ... 611!> \param r_shake ... 612!> \param dt ... 613!> \param ishake ... 614!> \param max_sigma ... 615!> \par History 616!> none 617! ************************************************************************************************** 618 SUBROUTINE shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, & 619 particle_set, pos, vel, r_shake, dt, ishake, max_sigma) 620 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 621 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 622 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 623 INTEGER, INTENT(IN) :: ng4x6, first_atom 624 TYPE(particle_type), POINTER :: particle_set(:) 625 REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :) 626 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake 627 REAL(kind=dp), INTENT(in) :: dt 628 INTEGER, INTENT(IN) :: ishake 629 REAL(KIND=dp), INTENT(INOUT) :: max_sigma 630 631 INTEGER :: iconst, index_a, index_b, index_c, & 632 index_d 633 REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, & 634 imass3, imass4, sigma 635 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, f_roll5, f_roll6, fc1, & 636 fc2, fc3, fc4, r0_12, r0_13, r0_14, r0_23, r0_24, r0_34, r1, r12, r13, r14, r2, r23, r24, & 637 r3, r34, r4, v1, v2, v3, v4, vec 638 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 639 REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp 640 TYPE(atomic_kind_type), POINTER :: atomic_kind 641 642 dtsqby2 = dt*dt*.5_dp 643 idtsq = 1.0_dp/dt/dt 644 dtby2 = dt*.5_dp 645 DO iconst = 1, ng4x6 646 IF (g4x6_list(iconst)%restraint%active) CYCLE 647 index_a = g4x6_list(iconst)%a + first_atom - 1 648 index_b = g4x6_list(iconst)%b + first_atom - 1 649 index_c = g4x6_list(iconst)%c + first_atom - 1 650 index_d = g4x6_list(iconst)%d + first_atom - 1 651 IF (ishake == 1) THEN 652 r0_12(:) = pos(:, index_a) - pos(:, index_b) 653 r0_13(:) = pos(:, index_a) - pos(:, index_c) 654 r0_23(:) = pos(:, index_b) - pos(:, index_c) 655 r0_14(:) = pos(:, index_a) - pos(:, index_d) 656 r0_24(:) = pos(:, index_b) - pos(:, index_d) 657 r0_34(:) = pos(:, index_c) - pos(:, index_d) 658 atomic_kind => particle_set(index_a)%atomic_kind 659 imass1 = 1.0_dp/atomic_kind%mass 660 atomic_kind => particle_set(index_b)%atomic_kind 661 imass2 = 1.0_dp/atomic_kind%mass 662 atomic_kind => particle_set(index_c)%atomic_kind 663 imass3 = 1.0_dp/atomic_kind%mass 664 atomic_kind => particle_set(index_d)%atomic_kind 665 imass4 = 1.0_dp/atomic_kind%mass 666 lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - & 667 lg4x6(iconst)%rb_old) 668 lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - & 669 lg4x6(iconst)%rc_old) 670 lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - & 671 lg4x6(iconst)%rd_old) 672 lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - & 673 lg4x6(iconst)%rc_old) 674 lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - & 675 lg4x6(iconst)%rd_old) 676 lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - & 677 lg4x6(iconst)%rd_old) 678 ! Check for fixed atom constraints 679 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, & 680 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 681 ! rotate fconst: 682 CALL matvec_3x3(f_roll1, r_shake, lg4x6(iconst)%fa) 683 CALL matvec_3x3(f_roll2, r_shake, lg4x6(iconst)%fb) 684 CALL matvec_3x3(f_roll3, r_shake, lg4x6(iconst)%fc) 685 CALL matvec_3x3(f_roll4, r_shake, lg4x6(iconst)%fd) 686 CALL matvec_3x3(f_roll5, r_shake, lg4x6(iconst)%fe) 687 CALL matvec_3x3(f_roll6, r_shake, lg4x6(iconst)%ff) 688 689 ! construct matrix 690 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, f_roll1) 691 amat(1, 2) = imass1*DOTPROD_3D(r0_12, f_roll2) 692 amat(1, 3) = imass1*DOTPROD_3D(r0_12, f_roll3) 693 amat(1, 4) = -imass2*DOTPROD_3D(r0_12, f_roll4) 694 amat(1, 5) = -imass2*DOTPROD_3D(r0_12, f_roll5) 695 amat(1, 6) = 0.0_dp 696 697 amat(2, 1) = imass1*DOTPROD_3D(r0_13, f_roll1) 698 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, f_roll2) 699 amat(2, 3) = imass1*DOTPROD_3D(r0_13, f_roll3) 700 amat(2, 4) = imass3*DOTPROD_3D(r0_13, f_roll4) 701 amat(2, 5) = 0.0_dp 702 amat(2, 6) = -imass3*DOTPROD_3D(r0_13, f_roll6) 703 704 amat(3, 1) = imass1*DOTPROD_3D(r0_14, f_roll1) 705 amat(3, 2) = imass1*DOTPROD_3D(r0_14, f_roll2) 706 amat(3, 3) = (imass1 + imass4)*DOTPROD_3D(r0_14, f_roll3) 707 amat(3, 4) = 0.0_dp 708 amat(3, 5) = imass4*DOTPROD_3D(r0_14, f_roll5) 709 amat(3, 6) = imass4*DOTPROD_3D(r0_14, f_roll6) 710 711 amat(4, 1) = -imass2*DOTPROD_3D(r0_23, f_roll1) 712 amat(4, 2) = imass3*DOTPROD_3D(r0_23, f_roll2) 713 amat(4, 3) = 0.0_dp 714 amat(4, 4) = (imass3 + imass2)*DOTPROD_3D(r0_23, f_roll4) 715 amat(4, 5) = imass2*DOTPROD_3D(r0_23, f_roll5) 716 amat(4, 6) = -imass3*DOTPROD_3D(r0_23, f_roll6) 717 718 amat(5, 1) = -imass2*DOTPROD_3D(r0_24, f_roll1) 719 amat(5, 2) = 0.0_dp 720 amat(5, 3) = imass4*DOTPROD_3D(r0_24, f_roll3) 721 amat(5, 4) = imass2*DOTPROD_3D(r0_24, f_roll4) 722 amat(5, 5) = (imass4 + imass2)*DOTPROD_3D(r0_24, f_roll5) 723 amat(5, 6) = imass4*DOTPROD_3D(r0_24, f_roll6) 724 725 amat(6, 1) = 0.0_dp 726 amat(6, 2) = -imass3*DOTPROD_3D(r0_34, f_roll2) 727 amat(6, 3) = imass4*DOTPROD_3D(r0_34, f_roll3) 728 amat(6, 4) = -imass3*DOTPROD_3D(r0_34, f_roll4) 729 amat(6, 5) = imass4*DOTPROD_3D(r0_34, f_roll5) 730 amat(6, 6) = (imass3 + imass4)*DOTPROD_3D(r0_34, f_roll6) 731 ! Store values 732 lg4x6(iconst)%r0_12 = r0_12 733 lg4x6(iconst)%r0_13 = r0_13 734 lg4x6(iconst)%r0_14 = r0_14 735 lg4x6(iconst)%r0_23 = r0_23 736 lg4x6(iconst)%r0_24 = r0_24 737 lg4x6(iconst)%r0_34 = r0_34 738 lg4x6(iconst)%f_roll1 = f_roll1 739 lg4x6(iconst)%f_roll2 = f_roll2 740 lg4x6(iconst)%f_roll3 = f_roll3 741 lg4x6(iconst)%f_roll4 = f_roll4 742 lg4x6(iconst)%f_roll5 = f_roll5 743 lg4x6(iconst)%f_roll6 = f_roll6 744 lg4x6(iconst)%amat = amat 745 lg4x6(iconst)%imass1 = imass1 746 lg4x6(iconst)%imass2 = imass2 747 lg4x6(iconst)%imass3 = imass3 748 lg4x6(iconst)%imass4 = imass4 749 lg4x6(iconst)%lambda_old = 0.0_dp 750 lg4x6(iconst)%del_lambda = 0.0_dp 751 ELSE 752 ! Retrieve values 753 r0_12 = lg4x6(iconst)%r0_12 754 r0_13 = lg4x6(iconst)%r0_13 755 r0_14 = lg4x6(iconst)%r0_14 756 r0_23 = lg4x6(iconst)%r0_23 757 r0_24 = lg4x6(iconst)%r0_24 758 r0_34 = lg4x6(iconst)%r0_34 759 f_roll1 = lg4x6(iconst)%f_roll1 760 f_roll2 = lg4x6(iconst)%f_roll2 761 f_roll3 = lg4x6(iconst)%f_roll3 762 f_roll4 = lg4x6(iconst)%f_roll4 763 f_roll5 = lg4x6(iconst)%f_roll5 764 f_roll6 = lg4x6(iconst)%f_roll6 765 amat = lg4x6(iconst)%amat 766 imass1 = lg4x6(iconst)%imass1 767 imass2 = lg4x6(iconst)%imass2 768 imass3 = lg4x6(iconst)%imass3 769 imass4 = lg4x6(iconst)%imass4 770 END IF 771 772 ! Iterate until convergence: 773 vec = lg4x6(iconst)%lambda(1)*f_roll1*(imass1 + imass2) + & 774 lg4x6(iconst)%lambda(2)*imass1*f_roll2 + & 775 lg4x6(iconst)%lambda(3)*imass1*f_roll3 - & 776 lg4x6(iconst)%lambda(4)*imass2*f_roll4 - & 777 lg4x6(iconst)%lambda(5)*imass2*f_roll5 778 bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab & 779 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12) 780 vec = lg4x6(iconst)%lambda(2)*f_roll2*(imass1 + imass3) + & 781 lg4x6(iconst)%lambda(1)*imass1*f_roll1 + & 782 lg4x6(iconst)%lambda(3)*imass1*f_roll3 + & 783 lg4x6(iconst)%lambda(4)*imass3*f_roll4 - & 784 lg4x6(iconst)%lambda(6)*imass3*f_roll6 785 bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac & 786 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13) 787 vec = lg4x6(iconst)%lambda(3)*f_roll3*(imass1 + imass4) + & 788 lg4x6(iconst)%lambda(1)*imass1*f_roll1 + & 789 lg4x6(iconst)%lambda(2)*imass1*f_roll2 + & 790 lg4x6(iconst)%lambda(5)*imass4*f_roll5 + & 791 lg4x6(iconst)%lambda(6)*imass4*f_roll6 792 bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad & 793 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_14, r0_14) 794 vec = lg4x6(iconst)%lambda(4)*f_roll4*(imass2 + imass3) - & 795 lg4x6(iconst)%lambda(1)*imass2*f_roll1 + & 796 lg4x6(iconst)%lambda(2)*imass3*f_roll2 + & 797 lg4x6(iconst)%lambda(5)*imass2*f_roll5 - & 798 lg4x6(iconst)%lambda(6)*imass3*f_roll6 799 bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc & 800 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23) 801 vec = lg4x6(iconst)%lambda(5)*f_roll5*(imass2 + imass4) - & 802 lg4x6(iconst)%lambda(1)*imass2*f_roll1 + & 803 lg4x6(iconst)%lambda(3)*imass4*f_roll3 + & 804 lg4x6(iconst)%lambda(4)*imass2*f_roll4 + & 805 lg4x6(iconst)%lambda(6)*imass4*f_roll6 806 bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd & 807 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_24, r0_24) 808 vec = lg4x6(iconst)%lambda(6)*f_roll6*(imass3 + imass4) - & 809 lg4x6(iconst)%lambda(2)*imass3*f_roll2 + & 810 lg4x6(iconst)%lambda(3)*imass4*f_roll3 - & 811 lg4x6(iconst)%lambda(4)*imass3*f_roll4 + & 812 lg4x6(iconst)%lambda(5)*imass4*f_roll5 813 bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd & 814 - dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_34, r0_34) 815 bvec = bvec*idtsq 816 817 ! get lambda 818 atemp = amat 819 CALL solve_system(atemp, 6, bvec) 820 lg4x6(iconst)%lambda(:) = bvec(:, 1) 821 lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - & 822 lg4x6(iconst)%lambda_old(:) 823 lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:) 824 825 fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + & 826 lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + & 827 lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc 828 fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + & 829 lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + & 830 lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe 831 fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - & 832 lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + & 833 lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff 834 fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - & 835 lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - & 836 lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff 837 CALL MATVEC_3x3(vec, r_shake, fc1) 838 r1(:) = pos(:, index_a) + imass1*dtsqby2*vec 839 CALL MATVEC_3x3(vec, r_shake, fc2) 840 r2(:) = pos(:, index_b) + imass2*dtsqby2*vec 841 CALL MATVEC_3x3(vec, r_shake, fc3) 842 r3(:) = pos(:, index_c) + imass3*dtsqby2*vec 843 CALL MATVEC_3x3(vec, r_shake, fc4) 844 r4(:) = pos(:, index_d) + imass4*dtsqby2*vec 845 CALL MATVEC_3x3(vec, r_shake, fc1) 846 v1(:) = vel(:, index_a) + imass1*dtby2*vec 847 CALL MATVEC_3x3(vec, r_shake, fc2) 848 v2(:) = vel(:, index_b) + imass2*dtby2*vec 849 CALL MATVEC_3x3(vec, r_shake, fc3) 850 v3(:) = vel(:, index_c) + imass3*dtby2*vec 851 CALL MATVEC_3x3(vec, r_shake, fc4) 852 v4(:) = vel(:, index_d) + imass4*dtby2*vec 853 854 r12 = r1 - r2 855 r13 = r1 - r3 856 r23 = r2 - r3 857 r14 = r1 - r4 858 r24 = r2 - r4 859 r34 = r3 - r4 860 ! compute the tolerance: 861 sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* & 862 g4x6_list(iconst)%dab 863 max_sigma = MAX(max_sigma, ABS(sigma)) 864 sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* & 865 g4x6_list(iconst)%dac 866 max_sigma = MAX(max_sigma, ABS(sigma)) 867 sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* & 868 g4x6_list(iconst)%dad 869 max_sigma = MAX(max_sigma, ABS(sigma)) 870 sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* & 871 g4x6_list(iconst)%dbc 872 max_sigma = MAX(max_sigma, ABS(sigma)) 873 sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* & 874 g4x6_list(iconst)%dbd 875 max_sigma = MAX(max_sigma, ABS(sigma)) 876 sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* & 877 g4x6_list(iconst)%dcd 878 max_sigma = MAX(max_sigma, ABS(sigma)) 879 880 ! update positions with full multiplier 881 pos(:, index_a) = r1(:) 882 pos(:, index_b) = r2(:) 883 pos(:, index_c) = r3(:) 884 pos(:, index_d) = r4(:) 885 886 ! update velocites with full multiplier 887 vel(:, index_a) = v1(:) 888 vel(:, index_b) = v2(:) 889 vel(:, index_c) = v3(:) 890 vel(:, index_d) = v4(:) 891 END DO 892 END SUBROUTINE shake_roll_4x6_low 893 894! ************************************************************************************************** 895!> \brief ... 896!> \param fixd_list ... 897!> \param g4x6_list ... 898!> \param lg4x6 ... 899!> \param first_atom ... 900!> \param particle_set ... 901!> \param vel ... 902!> \param dt ... 903!> \par History 904!> none 905! ************************************************************************************************** 906 SUBROUTINE rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 907 particle_set, vel, dt) 908 909 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 910 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 911 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 912 INTEGER, INTENT(IN) :: first_atom 913 TYPE(particle_type), POINTER :: particle_set(:) 914 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 915 REAL(kind=dp), INTENT(in) :: dt 916 917 INTEGER :: iconst, index_a, index_b, index_c, & 918 index_d 919 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, & 920 imass4, mass 921 REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r12, r13, r14, r23, & 922 r24, r34, v12, v13, v14, v23, v24, v34 923 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 924 REAL(KIND=dp), DIMENSION(6, 6) :: amat 925 TYPE(atomic_kind_type), POINTER :: atomic_kind 926 927 idt = 1.0_dp/dt 928 dtby2 = dt*.5_dp 929 DO iconst = 1, SIZE(g4x6_list) 930 IF (g4x6_list(iconst)%restraint%active) CYCLE 931 index_a = g4x6_list(iconst)%a + first_atom - 1 932 index_b = g4x6_list(iconst)%b + first_atom - 1 933 index_c = g4x6_list(iconst)%c + first_atom - 1 934 index_d = g4x6_list(iconst)%d + first_atom - 1 935 v12(:) = vel(:, index_a) - vel(:, index_b) 936 v13(:) = vel(:, index_a) - vel(:, index_c) 937 v14(:) = vel(:, index_a) - vel(:, index_d) 938 v23(:) = vel(:, index_b) - vel(:, index_c) 939 v24(:) = vel(:, index_b) - vel(:, index_d) 940 v34(:) = vel(:, index_c) - vel(:, index_d) 941 942 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:) 943 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:) 944 r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:) 945 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:) 946 r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:) 947 r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:) 948 atomic_kind => particle_set(index_a)%atomic_kind 949 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 950 imass1 = 1.0_dp/mass 951 atomic_kind => particle_set(index_b)%atomic_kind 952 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 953 imass2 = 1.0_dp/mass 954 atomic_kind => particle_set(index_c)%atomic_kind 955 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 956 imass3 = 1.0_dp/mass 957 atomic_kind => particle_set(index_d)%atomic_kind 958 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 959 imass4 = 1.0_dp/mass 960 lg4x6(iconst)%fa = -2.0_dp*r12 961 lg4x6(iconst)%fb = -2.0_dp*r13 962 lg4x6(iconst)%fc = -2.0_dp*r14 963 lg4x6(iconst)%fd = -2.0_dp*r23 964 lg4x6(iconst)%fe = -2.0_dp*r24 965 lg4x6(iconst)%ff = -2.0_dp*r34 966 ! Check for fixed atom constraints 967 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, & 968 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 969 ! construct matrix 970 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, lg4x6(iconst)%fa) 971 amat(1, 2) = imass1*DOTPROD_3D(r12, lg4x6(iconst)%fb) 972 amat(1, 3) = imass1*DOTPROD_3D(r12, lg4x6(iconst)%fc) 973 amat(1, 4) = -imass2*DOTPROD_3D(r12, lg4x6(iconst)%fd) 974 amat(1, 5) = -imass2*DOTPROD_3D(r12, lg4x6(iconst)%fe) 975 amat(1, 6) = 0.0_dp 976 977 amat(2, 1) = imass1*DOTPROD_3D(r13, lg4x6(iconst)%fa) 978 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, lg4x6(iconst)%fb) 979 amat(2, 3) = imass1*DOTPROD_3D(r13, lg4x6(iconst)%fc) 980 amat(2, 4) = imass3*DOTPROD_3D(r13, lg4x6(iconst)%fd) 981 amat(2, 5) = 0.0_dp 982 amat(2, 6) = -imass3*DOTPROD_3D(r13, lg4x6(iconst)%ff) 983 984 amat(3, 1) = imass1*DOTPROD_3D(r14, lg4x6(iconst)%fa) 985 amat(3, 2) = imass1*DOTPROD_3D(r14, lg4x6(iconst)%fb) 986 amat(3, 3) = (imass1 + imass4)*DOTPROD_3D(r14, lg4x6(iconst)%fc) 987 amat(3, 4) = 0.0_dp 988 amat(3, 5) = imass4*DOTPROD_3D(r14, lg4x6(iconst)%fe) 989 amat(3, 6) = imass4*DOTPROD_3D(r14, lg4x6(iconst)%ff) 990 991 amat(4, 1) = -imass2*DOTPROD_3D(r23, lg4x6(iconst)%fa) 992 amat(4, 2) = imass3*DOTPROD_3D(r23, lg4x6(iconst)%fb) 993 amat(4, 3) = 0.0_dp 994 amat(4, 4) = (imass3 + imass2)*DOTPROD_3D(r23, lg4x6(iconst)%fd) 995 amat(4, 5) = imass2*DOTPROD_3D(r23, lg4x6(iconst)%fe) 996 amat(4, 6) = -imass3*DOTPROD_3D(r23, lg4x6(iconst)%ff) 997 998 amat(5, 1) = -imass2*DOTPROD_3D(r24, lg4x6(iconst)%fa) 999 amat(5, 2) = 0.0_dp 1000 amat(5, 3) = imass4*DOTPROD_3D(r24, lg4x6(iconst)%fc) 1001 amat(5, 4) = imass2*DOTPROD_3D(r24, lg4x6(iconst)%fd) 1002 amat(5, 5) = (imass4 + imass2)*DOTPROD_3D(r24, lg4x6(iconst)%fe) 1003 amat(5, 6) = imass4*DOTPROD_3D(r24, lg4x6(iconst)%ff) 1004 1005 amat(6, 1) = 0.0_dp 1006 amat(6, 2) = -imass3*DOTPROD_3D(r34, lg4x6(iconst)%fb) 1007 amat(6, 3) = imass4*DOTPROD_3D(r34, lg4x6(iconst)%fc) 1008 amat(6, 4) = -imass3*DOTPROD_3D(r34, lg4x6(iconst)%fd) 1009 amat(6, 5) = imass4*DOTPROD_3D(r34, lg4x6(iconst)%fe) 1010 amat(6, 6) = (imass3 + imass4)*DOTPROD_3D(r34, lg4x6(iconst)%ff) 1011 1012 ! construct solution vector 1013 bvec(1, 1) = DOTPROD_3D(r12, v12) 1014 bvec(2, 1) = DOTPROD_3D(r13, v13) 1015 bvec(3, 1) = DOTPROD_3D(r14, v14) 1016 bvec(4, 1) = DOTPROD_3D(r23, v23) 1017 bvec(5, 1) = DOTPROD_3D(r24, v24) 1018 bvec(6, 1) = DOTPROD_3D(r34, v34) 1019 bvec = -bvec*2.0_dp*idt 1020 1021 ! get lambda 1022 CALL solve_system(amat, 6, bvec) 1023 lg4x6(iconst)%lambda(:) = bvec(:, 1) 1024 1025 fc1 = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + & 1026 lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb + & 1027 lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc 1028 fc2 = -lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + & 1029 lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + & 1030 lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe 1031 fc3 = -lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb - & 1032 lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + & 1033 lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff 1034 fc4 = -lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc - & 1035 lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe - & 1036 lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff 1037 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:) 1038 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:) 1039 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:) 1040 vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:) 1041 END DO 1042 END SUBROUTINE rattle_4x6_low 1043 1044! ************************************************************************************************** 1045!> \brief ... 1046!> \param fixd_list ... 1047!> \param g4x6_list ... 1048!> \param lg4x6 ... 1049!> \param first_atom ... 1050!> \param particle_set ... 1051!> \param vel ... 1052!> \param r_rattle ... 1053!> \param dt ... 1054!> \param veps ... 1055!> \par History 1056!> none 1057! ************************************************************************************************** 1058 SUBROUTINE rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, & 1059 particle_set, vel, r_rattle, dt, veps) 1060 1061 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 1062 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 1063 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:) 1064 INTEGER, INTENT(IN) :: first_atom 1065 TYPE(particle_type), POINTER :: particle_set(:) 1066 REAL(KIND=dp), INTENT(INOUT) :: vel(:, :) 1067 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle 1068 REAL(kind=dp), INTENT(in) :: dt 1069 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps 1070 1071 INTEGER :: iconst, index_a, index_b, index_c, & 1072 index_d 1073 REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, & 1074 imass4, mass 1075 REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, & 1076 f_roll5, f_roll6, fc1, fc2, fc3, fc4, & 1077 r12, r13, r14, r23, r24, r34, v12, & 1078 v13, v14, v23, v24, v34, vec 1079 REAL(KIND=dp), DIMENSION(6) :: lambda 1080 REAL(KIND=dp), DIMENSION(6, 1) :: bvec 1081 REAL(KIND=dp), DIMENSION(6, 6) :: amat 1082 TYPE(atomic_kind_type), POINTER :: atomic_kind 1083 1084 idt = 1.0_dp/dt 1085 dtby2 = dt*.5_dp 1086 DO iconst = 1, SIZE(g4x6_list) 1087 IF (g4x6_list(iconst)%restraint%active) CYCLE 1088 index_a = g4x6_list(iconst)%a + first_atom - 1 1089 index_b = g4x6_list(iconst)%b + first_atom - 1 1090 index_c = g4x6_list(iconst)%c + first_atom - 1 1091 index_d = g4x6_list(iconst)%d + first_atom - 1 1092 v12(:) = vel(:, index_a) - vel(:, index_b) 1093 v13(:) = vel(:, index_a) - vel(:, index_c) 1094 v14(:) = vel(:, index_a) - vel(:, index_d) 1095 v23(:) = vel(:, index_b) - vel(:, index_c) 1096 v24(:) = vel(:, index_b) - vel(:, index_d) 1097 v34(:) = vel(:, index_c) - vel(:, index_d) 1098 1099 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:) 1100 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:) 1101 r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:) 1102 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:) 1103 r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:) 1104 r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:) 1105 atomic_kind => particle_set(index_a)%atomic_kind 1106 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 1107 imass1 = 1.0_dp/mass 1108 atomic_kind => particle_set(index_b)%atomic_kind 1109 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 1110 imass2 = 1.0_dp/mass 1111 atomic_kind => particle_set(index_c)%atomic_kind 1112 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 1113 imass3 = 1.0_dp/mass 1114 atomic_kind => particle_set(index_d)%atomic_kind 1115 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass) 1116 imass4 = 1.0_dp/mass 1117 lg4x6(iconst)%fa = -2.0_dp*r12 1118 lg4x6(iconst)%fb = -2.0_dp*r13 1119 lg4x6(iconst)%fc = -2.0_dp*r14 1120 lg4x6(iconst)%fd = -2.0_dp*r23 1121 lg4x6(iconst)%fe = -2.0_dp*r24 1122 lg4x6(iconst)%ff = -2.0_dp*r34 1123 ! Check for fixed atom constraints 1124 CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, & 1125 index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst)) 1126 ! roll the fc 1127 CALL MATVEC_3x3(f_roll1, r_rattle, lg4x6(iconst)%fa) 1128 CALL MATVEC_3x3(f_roll2, r_rattle, lg4x6(iconst)%fb) 1129 CALL MATVEC_3x3(f_roll3, r_rattle, lg4x6(iconst)%fc) 1130 CALL MATVEC_3x3(f_roll4, r_rattle, lg4x6(iconst)%fd) 1131 CALL MATVEC_3x3(f_roll5, r_rattle, lg4x6(iconst)%fe) 1132 CALL MATVEC_3x3(f_roll6, r_rattle, lg4x6(iconst)%ff) 1133 ! construct matrix 1134 amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, f_roll1) 1135 amat(1, 2) = imass1*DOTPROD_3D(r12, f_roll2) 1136 amat(1, 3) = imass1*DOTPROD_3D(r12, f_roll3) 1137 amat(1, 4) = -imass2*DOTPROD_3D(r12, f_roll4) 1138 amat(1, 5) = -imass2*DOTPROD_3D(r12, f_roll5) 1139 amat(1, 6) = 0.0_dp 1140 1141 amat(2, 1) = imass1*DOTPROD_3D(r13, f_roll1) 1142 amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, f_roll2) 1143 amat(2, 3) = imass1*DOTPROD_3D(r13, f_roll3) 1144 amat(2, 4) = imass3*DOTPROD_3D(r13, f_roll4) 1145 amat(2, 5) = 0.0_dp 1146 amat(2, 6) = -imass3*DOTPROD_3D(r13, f_roll6) 1147 1148 amat(3, 1) = imass1*DOTPROD_3D(r14, f_roll1) 1149 amat(3, 2) = imass1*DOTPROD_3D(r14, f_roll2) 1150 amat(3, 3) = (imass1 + imass4)*DOTPROD_3D(r14, f_roll3) 1151 amat(3, 4) = 0.0_dp 1152 amat(3, 5) = imass4*DOTPROD_3D(r14, f_roll5) 1153 amat(3, 6) = imass4*DOTPROD_3D(r14, f_roll6) 1154 1155 amat(4, 1) = -imass2*DOTPROD_3D(r23, f_roll1) 1156 amat(4, 2) = imass3*DOTPROD_3D(r23, f_roll2) 1157 amat(4, 3) = 0.0_dp 1158 amat(4, 4) = (imass3 + imass2)*DOTPROD_3D(r23, f_roll4) 1159 amat(4, 5) = imass2*DOTPROD_3D(r23, f_roll5) 1160 amat(4, 6) = -imass3*DOTPROD_3D(r23, f_roll6) 1161 1162 amat(5, 1) = -imass2*DOTPROD_3D(r24, f_roll1) 1163 amat(5, 2) = 0.0_dp 1164 amat(5, 3) = imass4*DOTPROD_3D(r24, f_roll3) 1165 amat(5, 4) = imass2*DOTPROD_3D(r24, f_roll4) 1166 amat(5, 5) = (imass4 + imass2)*DOTPROD_3D(r24, f_roll5) 1167 amat(5, 6) = imass4*DOTPROD_3D(r24, f_roll6) 1168 1169 amat(6, 1) = 0.0_dp 1170 amat(6, 2) = -imass3*DOTPROD_3D(r34, f_roll2) 1171 amat(6, 3) = imass4*DOTPROD_3D(r34, f_roll3) 1172 amat(6, 4) = -imass3*DOTPROD_3D(r34, f_roll4) 1173 amat(6, 5) = imass4*DOTPROD_3D(r34, f_roll5) 1174 amat(6, 6) = (imass3 + imass4)*DOTPROD_3D(r34, f_roll6) 1175 1176 ! construct solution vector 1177 CALL matvec_3x3(vec, veps, r12) 1178 bvec(1, 1) = DOTPROD_3D(r12, v12 + vec) 1179 CALL matvec_3x3(vec, veps, r13) 1180 bvec(2, 1) = DOTPROD_3D(r13, v13 + vec) 1181 CALL matvec_3x3(vec, veps, r14) 1182 bvec(3, 1) = DOTPROD_3D(r14, v14 + vec) 1183 CALL matvec_3x3(vec, veps, r23) 1184 bvec(4, 1) = DOTPROD_3D(r23, v23 + vec) 1185 CALL matvec_3x3(vec, veps, r24) 1186 bvec(5, 1) = DOTPROD_3D(r24, v24 + vec) 1187 CALL matvec_3x3(vec, veps, r34) 1188 bvec(6, 1) = DOTPROD_3D(r34, v34 + vec) 1189 bvec = -bvec*2.0_dp*idt 1190 1191 ! get lambda 1192 CALL solve_system(amat, 6, bvec) 1193 lambda(:) = bvec(:, 1) 1194 lg4x6(iconst)%lambda(:) = lambda 1195 1196 fc1 = lambda(1)*f_roll1 + & 1197 lambda(2)*f_roll2 + & 1198 lambda(3)*f_roll3 1199 fc2 = -lambda(1)*f_roll1 + & 1200 lambda(4)*f_roll4 + & 1201 lambda(5)*f_roll5 1202 fc3 = -lambda(2)*f_roll2 - & 1203 lambda(4)*f_roll4 + & 1204 lambda(6)*f_roll6 1205 fc4 = -lambda(3)*f_roll3 - & 1206 lambda(5)*f_roll5 - & 1207 lambda(6)*f_roll6 1208 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:) 1209 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:) 1210 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:) 1211 vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:) 1212 END DO 1213 END SUBROUTINE rattle_roll_4x6_low 1214 1215END MODULE constraint_4x6 1216