1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Handles all possible kinds of restraints in CP2K 8!> \author Teodoro Laino 08.2006 [tlaino] - University of Zurich 9!> \par History 10!> Teodoro Laino [tlaino] - 11.2008 : Improved the fixd_list restraints 11! ************************************************************************************************** 12MODULE restraint 13 USE cell_types, ONLY: cell_type,& 14 use_perd_x,& 15 use_perd_xy,& 16 use_perd_xyz,& 17 use_perd_xz,& 18 use_perd_y,& 19 use_perd_yz,& 20 use_perd_z 21 USE colvar_methods, ONLY: colvar_eval_mol_f 22 USE colvar_types, ONLY: colvar_counters,& 23 diff_colvar 24 USE constraint_fxd, ONLY: create_local_fixd_list,& 25 release_local_fixd_list 26 USE cp_subsys_types, ONLY: cp_subsys_get,& 27 cp_subsys_type 28 USE distribution_1d_types, ONLY: distribution_1d_type 29 USE force_env_types, ONLY: force_env_get,& 30 force_env_set,& 31 force_env_type 32 USE kinds, ONLY: dp 33 USE message_passing, ONLY: mp_sum 34 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 35 USE molecule_kind_types, ONLY: colvar_constraint_type,& 36 fixd_constraint_type,& 37 g3x3_constraint_type,& 38 g4x6_constraint_type,& 39 get_molecule_kind,& 40 local_fixd_constraint_type,& 41 molecule_kind_type 42 USE molecule_list_types, ONLY: molecule_list_type 43 USE molecule_types, ONLY: get_molecule,& 44 global_constraint_type,& 45 local_colvar_constraint_type,& 46 molecule_type 47 USE particle_list_types, ONLY: particle_list_type 48 USE particle_types, ONLY: particle_type,& 49 update_particle_set 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 PRIVATE 55 PUBLIC :: restraint_control 56 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'restraint' 57 58CONTAINS 59 60! ************************************************************************************************** 61!> \brief Computes restraints 62!> \param force_env ... 63!> \author Teodoro Laino 08.2006 [tlaino] 64! ************************************************************************************************** 65 SUBROUTINE restraint_control(force_env) 66 67 TYPE(force_env_type), POINTER :: force_env 68 69 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_control', & 70 routineP = moduleN//':'//routineN 71 72 INTEGER :: handle, i, ifixd, ii, ikind, imol, & 73 iparticle, n3x3con_restraint, & 74 n4x6con_restraint, n_restraint, nkind, & 75 nmol_per_kind 76 REAL(KIND=dp) :: energy_3x3, energy_4x6, energy_colv, & 77 energy_fixd, extended_energies, k0, & 78 rab(3), rab2, targ(3) 79 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force 80 TYPE(cell_type), POINTER :: cell 81 TYPE(colvar_counters) :: ncolv 82 TYPE(cp_subsys_type), POINTER :: subsys 83 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles 84 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 85 TYPE(global_constraint_type), POINTER :: gci 86 TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:) 87 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 88 TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:) 89 TYPE(molecule_list_type), POINTER :: molecules 90 TYPE(molecule_type), POINTER :: molecule, molecule_set(:) 91 TYPE(particle_list_type), POINTER :: particles 92 TYPE(particle_type), POINTER :: particle_set(:) 93 94 NULLIFY (cell, subsys, local_molecules, local_particles, fixd_list, molecule_kinds, & 95 molecules, molecule_kind, molecule_kind_set, molecule, molecule_set, particles, & 96 particle_set, gci, lfixd_list) 97 CALL timeset(routineN, handle) 98 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 99 energy_4x6 = 0.0_dp 100 energy_3x3 = 0.0_dp 101 energy_colv = 0.0_dp 102 energy_fixd = 0.0_dp 103 n_restraint = 0 104 CALL cp_subsys_get(subsys=subsys, particles=particles, molecules=molecules, & 105 local_particles=local_particles, local_molecules=local_molecules, & 106 gci=gci, molecule_kinds=molecule_kinds) 107 108 nkind = molecule_kinds%n_els 109 particle_set => particles%els 110 molecule_set => molecules%els 111 molecule_kind_set => molecule_kinds%els 112 113 ! Intramolecular Restraints 114 ALLOCATE (force(3, SIZE(particle_set))) 115 force = 0.0_dp 116 117 ! Create the list of locally fixed atoms 118 CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles) 119 120 DO ifixd = 1, SIZE(lfixd_list) 121 ikind = lfixd_list(ifixd)%ikind 122 ii = lfixd_list(ifixd)%ifixd_index 123 molecule_kind => molecule_kind_set(ikind) 124 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list) 125 IF (fixd_list(ii)%restraint%active) THEN 126 n_restraint = n_restraint + 1 127 iparticle = fixd_list(ii)%fixd 128 k0 = fixd_list(ii)%restraint%k0 129 targ = fixd_list(ii)%coord 130 rab = 0.0_dp 131 SELECT CASE (fixd_list(ii)%itype) 132 CASE (use_perd_x) 133 rab(1) = particle_set(iparticle)%r(1) - targ(1) 134 CASE (use_perd_y) 135 rab(2) = particle_set(iparticle)%r(2) - targ(2) 136 CASE (use_perd_z) 137 rab(3) = particle_set(iparticle)%r(3) - targ(3) 138 CASE (use_perd_xy) 139 rab(1) = particle_set(iparticle)%r(1) - targ(1) 140 rab(2) = particle_set(iparticle)%r(2) - targ(2) 141 CASE (use_perd_xz) 142 rab(1) = particle_set(iparticle)%r(1) - targ(1) 143 rab(3) = particle_set(iparticle)%r(3) - targ(3) 144 CASE (use_perd_yz) 145 rab(2) = particle_set(iparticle)%r(2) - targ(2) 146 rab(3) = particle_set(iparticle)%r(3) - targ(3) 147 CASE (use_perd_xyz) 148 rab = particle_set(iparticle)%r - targ 149 END SELECT 150 rab2 = DOT_PRODUCT(rab, rab) 151 ! Energy 152 energy_fixd = energy_fixd + k0*rab2 153 ! Forces 154 force(:, iparticle) = force(:, iparticle) - 2.0_dp*k0*rab 155 END IF 156 END DO 157 CALL release_local_fixd_list(lfixd_list) 158 159 ! Loop over other kind of Restraints 160 MOL: DO ikind = 1, nkind 161 molecule_kind => molecule_kind_set(ikind) 162 nmol_per_kind = local_molecules%n_el(ikind) 163 DO imol = 1, nmol_per_kind 164 i = local_molecules%list(ikind)%array(imol) 165 molecule => molecule_set(i) 166 molecule_kind => molecule%molecule_kind 167 168 CALL get_molecule_kind(molecule_kind, & 169 ncolv=ncolv, & 170 ng3x3_restraint=n3x3con_restraint, & 171 ng4x6_restraint=n4x6con_restraint) 172 ! 3x3 173 IF (n3x3con_restraint /= 0) THEN 174 n_restraint = n_restraint + n3x3con_restraint 175 CALL restraint_3x3_int(molecule, particle_set, energy_3x3, force) 176 ENDIF 177 ! 4x6 178 IF (n4x6con_restraint /= 0) THEN 179 n_restraint = n_restraint + n4x6con_restraint 180 CALL restraint_4x6_int(molecule, particle_set, energy_4x6, force) 181 ENDIF 182 ! collective variables 183 IF (ncolv%nrestraint /= 0) THEN 184 n_restraint = n_restraint + ncolv%nrestraint 185 CALL restraint_colv_int(molecule, particle_set, cell, energy_colv, force) 186 ENDIF 187 END DO 188 END DO MOL 189 CALL mp_sum(n_restraint, force_env%para_env%group) 190 IF (n_restraint > 0) THEN 191 CALL mp_sum(energy_fixd, force_env%para_env%group) 192 CALL mp_sum(energy_3x3, force_env%para_env%group) 193 CALL mp_sum(energy_4x6, force_env%para_env%group) 194 CALL mp_sum(energy_colv, force_env%para_env%group) 195 CALL update_particle_set(particle_set, force_env%para_env%group, for=force, add=.TRUE.) 196 force = 0.0_dp 197 n_restraint = 0 198 ENDIF 199 ! Intermolecular Restraints 200 IF (ASSOCIATED(gci)) THEN 201 IF (gci%nrestraint > 0) THEN 202 ! 3x3 203 IF (gci%ng3x3_restraint /= 0) THEN 204 n_restraint = n_restraint + gci%ng3x3_restraint 205 CALL restraint_3x3_ext(gci, particle_set, energy_3x3, force) 206 ENDIF 207 ! 4x6 208 IF (gci%ng4x6_restraint /= 0) THEN 209 n_restraint = n_restraint + gci%ng4x6_restraint 210 CALL restraint_4x6_ext(gci, particle_set, energy_4x6, force) 211 ENDIF 212 ! collective variables 213 IF (gci%ncolv%nrestraint /= 0) THEN 214 n_restraint = n_restraint + gci%ncolv%nrestraint 215 CALL restraint_colv_ext(gci, particle_set, cell, energy_colv, force) 216 ENDIF 217 DO iparticle = 1, SIZE(particle_set) 218 particle_set(iparticle)%f = particle_set(iparticle)%f + force(:, iparticle) 219 END DO 220 END IF 221 END IF 222 DEALLOCATE (force) 223 224 ! Store restraint energies 225 CALL force_env_get(force_env=force_env, additional_potential=extended_energies) 226 extended_energies = extended_energies + energy_3x3 + & 227 energy_fixd + & 228 energy_4x6 + & 229 energy_colv 230 CALL force_env_set(force_env=force_env, additional_potential=extended_energies) 231 CALL timestop(handle) 232 233 END SUBROUTINE restraint_control 234 235! ************************************************************************************************** 236!> \brief Computes restraints 3x3 - Intramolecular 237!> \param molecule ... 238!> \param particle_set ... 239!> \param energy ... 240!> \param force ... 241!> \author Teodoro Laino 08.2006 [tlaino] 242! ************************************************************************************************** 243 SUBROUTINE restraint_3x3_int(molecule, particle_set, energy, force) 244 245 TYPE(molecule_type), POINTER :: molecule 246 TYPE(particle_type), POINTER :: particle_set(:) 247 REAL(KIND=dp), INTENT(INOUT) :: energy 248 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 249 250 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_3x3_int', & 251 routineP = moduleN//':'//routineN 252 253 INTEGER :: first_atom, ng3x3 254 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 255 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 256 TYPE(molecule_kind_type), POINTER :: molecule_kind 257 258 molecule_kind => molecule%molecule_kind 259 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, g3x3_list=g3x3_list, & 260 fixd_list=fixd_list) 261 CALL get_molecule(molecule, first_atom=first_atom) 262 CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, & 263 energy, force) 264 265 END SUBROUTINE restraint_3x3_int 266 267! ************************************************************************************************** 268!> \brief Computes restraints 4x6 - Intramolecular 269!> \param molecule ... 270!> \param particle_set ... 271!> \param energy ... 272!> \param force ... 273!> \author Teodoro Laino 08.2006 [tlaino] 274! ************************************************************************************************** 275 SUBROUTINE restraint_4x6_int(molecule, particle_set, energy, force) 276 277 TYPE(molecule_type), POINTER :: molecule 278 TYPE(particle_type), POINTER :: particle_set(:) 279 REAL(KIND=dp), INTENT(INOUT) :: energy 280 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 281 282 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_4x6_int', & 283 routineP = moduleN//':'//routineN 284 285 INTEGER :: first_atom, ng4x6 286 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 287 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 288 TYPE(molecule_kind_type), POINTER :: molecule_kind 289 290 molecule_kind => molecule%molecule_kind 291 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, g4x6_list=g4x6_list, & 292 fixd_list=fixd_list) 293 CALL get_molecule(molecule, first_atom=first_atom) 294 CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, & 295 energy, force) 296 297 END SUBROUTINE restraint_4x6_int 298 299! ************************************************************************************************** 300!> \brief Computes restraints colv - Intramolecular 301!> \param molecule ... 302!> \param particle_set ... 303!> \param cell ... 304!> \param energy ... 305!> \param force ... 306!> \author Teodoro Laino 08.2006 [tlaino] 307! ************************************************************************************************** 308 SUBROUTINE restraint_colv_int(molecule, particle_set, cell, energy, force) 309 310 TYPE(molecule_type), POINTER :: molecule 311 TYPE(particle_type), POINTER :: particle_set(:) 312 TYPE(cell_type), POINTER :: cell 313 REAL(KIND=dp), INTENT(INOUT) :: energy 314 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 315 316 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_colv_int', & 317 routineP = moduleN//':'//routineN 318 319 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 320 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 321 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 322 TYPE(molecule_kind_type), POINTER :: molecule_kind 323 324 NULLIFY (fixd_list) 325 326 molecule_kind => molecule%molecule_kind 327 CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list) 328 CALL get_molecule(molecule, lcolv=lcolv) 329 CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, & 330 cell, energy, force) 331 332 END SUBROUTINE restraint_colv_int 333 334! ************************************************************************************************** 335!> \brief Computes restraints 3x3 - Intermolecular 336!> \param gci ... 337!> \param particle_set ... 338!> \param energy ... 339!> \param force ... 340!> \author Teodoro Laino 08.2006 [tlaino] 341! ************************************************************************************************** 342 SUBROUTINE restraint_3x3_ext(gci, particle_set, energy, force) 343 344 TYPE(global_constraint_type), POINTER :: gci 345 TYPE(particle_type), POINTER :: particle_set(:) 346 REAL(KIND=dp), INTENT(INOUT) :: energy 347 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 348 349 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_3x3_ext', & 350 routineP = moduleN//':'//routineN 351 352 INTEGER :: first_atom, ng3x3 353 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 354 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 355 356 first_atom = 1 357 ng3x3 = gci%ng3x3 358 g3x3_list => gci%g3x3_list 359 fixd_list => gci%fixd_list 360 CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, & 361 energy, force) 362 363 END SUBROUTINE restraint_3x3_ext 364 365! ************************************************************************************************** 366!> \brief Computes restraints 4x6 - Intermolecular 367!> \param gci ... 368!> \param particle_set ... 369!> \param energy ... 370!> \param force ... 371!> \author Teodoro Laino 08.2006 [tlaino] 372! ************************************************************************************************** 373 SUBROUTINE restraint_4x6_ext(gci, particle_set, energy, force) 374 375 TYPE(global_constraint_type), POINTER :: gci 376 TYPE(particle_type), POINTER :: particle_set(:) 377 REAL(KIND=dp), INTENT(INOUT) :: energy 378 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 379 380 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_4x6_ext', & 381 routineP = moduleN//':'//routineN 382 383 INTEGER :: first_atom, ng4x6 384 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 385 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 386 387 first_atom = 1 388 ng4x6 = gci%ng4x6 389 g4x6_list => gci%g4x6_list 390 fixd_list => gci%fixd_list 391 CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, & 392 energy, force) 393 394 END SUBROUTINE restraint_4x6_ext 395 396! ************************************************************************************************** 397!> \brief Computes restraints colv - Intermolecular 398!> \param gci ... 399!> \param particle_set ... 400!> \param cell ... 401!> \param energy ... 402!> \param force ... 403!> \author Teodoro Laino 08.2006 [tlaino] 404! ************************************************************************************************** 405 SUBROUTINE restraint_colv_ext(gci, particle_set, cell, energy, force) 406 407 TYPE(global_constraint_type), POINTER :: gci 408 TYPE(particle_type), POINTER :: particle_set(:) 409 TYPE(cell_type), POINTER :: cell 410 REAL(KIND=dp), INTENT(INOUT) :: energy 411 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 412 413 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_colv_ext', & 414 routineP = moduleN//':'//routineN 415 416 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 417 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 418 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 419 420 colv_list => gci%colv_list 421 fixd_list => gci%fixd_list 422 lcolv => gci%lcolv 423 CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, & 424 cell, energy, force) 425 426 END SUBROUTINE restraint_colv_ext 427 428! ************************************************************************************************** 429!> \brief Computes restraints 3x3 - Real 3x3 restraints 430!> \param ng3x3 ... 431!> \param g3x3_list ... 432!> \param fixd_list ... 433!> \param first_atom ... 434!> \param particle_set ... 435!> \param energy ... 436!> \param force ... 437!> \author Teodoro Laino 08.2006 [tlaino] 438! ************************************************************************************************** 439 SUBROUTINE restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, & 440 particle_set, energy, force) 441 442 INTEGER :: ng3x3 443 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:) 444 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 445 INTEGER, INTENT(IN) :: first_atom 446 TYPE(particle_type), POINTER :: particle_set(:) 447 REAL(KIND=dp), INTENT(INOUT) :: energy 448 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 449 450 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_3x3_low', & 451 routineP = moduleN//':'//routineN 452 453 INTEGER :: iconst, index_a, index_b, index_c 454 REAL(KIND=dp) :: k, rab, rac, rbc, tab, tac, tbc 455 REAL(KIND=dp), DIMENSION(3) :: r0_12, r0_13, r0_23 456 457 DO iconst = 1, ng3x3 458 IF (.NOT. g3x3_list(iconst)%restraint%active) CYCLE 459 index_a = g3x3_list(iconst)%a + first_atom - 1 460 index_b = g3x3_list(iconst)%b + first_atom - 1 461 index_c = g3x3_list(iconst)%c + first_atom - 1 462 r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r 463 r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r 464 r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r 465 466 rab = SQRT(DOT_PRODUCT(r0_12, r0_12)) 467 rac = SQRT(DOT_PRODUCT(r0_13, r0_13)) 468 rbc = SQRT(DOT_PRODUCT(r0_23, r0_23)) 469 tab = rab - g3x3_list(ng3x3)%dab 470 tac = rac - g3x3_list(ng3x3)%dac 471 tbc = rbc - g3x3_list(ng3x3)%dbc 472 k = g3x3_list(iconst)%restraint%k0 473 ! Update Energy 474 energy = energy + k*(tab**2 + tac**2 + tbc**2) 475 ! Update Forces 476 force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac) 477 force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc) 478 force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc) 479 ! Fixed atoms 480 IF (ASSOCIATED(fixd_list)) THEN 481 IF (SIZE(fixd_list) > 0) THEN 482 IF (ANY(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp 483 IF (ANY(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp 484 IF (ANY(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp 485 END IF 486 END IF 487 END DO 488 489 END SUBROUTINE restraint_3x3_low 490 491! ************************************************************************************************** 492!> \brief Computes restraints 4x6 - Real 4x6 restraints 493!> \param ng4x6 ... 494!> \param g4x6_list ... 495!> \param fixd_list ... 496!> \param first_atom ... 497!> \param particle_set ... 498!> \param energy ... 499!> \param force ... 500!> \author Teodoro Laino 08.2006 [tlaino] 501! ************************************************************************************************** 502 SUBROUTINE restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, & 503 particle_set, energy, force) 504 505 INTEGER :: ng4x6 506 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:) 507 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 508 INTEGER, INTENT(IN) :: first_atom 509 TYPE(particle_type), POINTER :: particle_set(:) 510 REAL(KIND=dp), INTENT(INOUT) :: energy 511 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 512 513 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_4x6_low', & 514 routineP = moduleN//':'//routineN 515 516 INTEGER :: iconst, index_a, index_b, index_c, & 517 index_d 518 REAL(KIND=dp) :: k, rab, rac, rad, rbc, rbd, rcd, tab, & 519 tac, tad, tbc, tbd, tcd 520 REAL(KIND=dp), DIMENSION(3) :: r0_12, r0_13, r0_14, r0_23, r0_24, r0_34 521 522 DO iconst = 1, ng4x6 523 IF (.NOT. g4x6_list(iconst)%restraint%active) CYCLE 524 index_a = g4x6_list(iconst)%a + first_atom - 1 525 index_b = g4x6_list(iconst)%b + first_atom - 1 526 index_c = g4x6_list(iconst)%c + first_atom - 1 527 index_d = g4x6_list(iconst)%d + first_atom - 1 528 r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r 529 r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r 530 r0_14(:) = particle_set(index_a)%r - particle_set(index_d)%r 531 r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r 532 r0_24(:) = particle_set(index_b)%r - particle_set(index_d)%r 533 r0_34(:) = particle_set(index_c)%r - particle_set(index_d)%r 534 535 rab = SQRT(DOT_PRODUCT(r0_12, r0_12)) 536 rac = SQRT(DOT_PRODUCT(r0_13, r0_13)) 537 rad = SQRT(DOT_PRODUCT(r0_14, r0_14)) 538 rbc = SQRT(DOT_PRODUCT(r0_23, r0_23)) 539 rbd = SQRT(DOT_PRODUCT(r0_24, r0_24)) 540 rcd = SQRT(DOT_PRODUCT(r0_34, r0_34)) 541 542 tab = rab - g4x6_list(ng4x6)%dab 543 tac = rac - g4x6_list(ng4x6)%dac 544 tad = rad - g4x6_list(ng4x6)%dad 545 tbc = rbc - g4x6_list(ng4x6)%dbc 546 tbd = rbd - g4x6_list(ng4x6)%dbd 547 tcd = rcd - g4x6_list(ng4x6)%dcd 548 549 k = g4x6_list(iconst)%restraint%k0 550 ! Update Energy 551 energy = energy + k*(tab**2 + tac**2 + tad**2 + tbc**2 + tbd**2 + tcd**2) 552 ! Update Forces 553 force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac + r0_14/rad*tad) 554 force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc + r0_24/rbd*tbd) 555 force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc + r0_34/rcd*tcd) 556 force(:, index_d) = force(:, index_d) - 2.0_dp*k*(-r0_14/rad*tad - r0_24/rbd*tbd - r0_34/rcd*tcd) 557 ! Fixed atoms 558 IF (ASSOCIATED(fixd_list)) THEN 559 IF (SIZE(fixd_list) > 0) THEN 560 IF (ANY(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp 561 IF (ANY(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp 562 IF (ANY(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp 563 IF (ANY(fixd_list(:)%fixd == index_d)) force(:, index_d) = 0.0_dp 564 END IF 565 END IF 566 END DO 567 568 END SUBROUTINE restraint_4x6_low 569 570! ************************************************************************************************** 571!> \brief Computes restraints colv - Real COLVAR restraints 572!> \param colv_list ... 573!> \param fixd_list ... 574!> \param lcolv ... 575!> \param particle_set ... 576!> \param cell ... 577!> \param energy ... 578!> \param force ... 579!> \author Teodoro Laino 08.2006 [tlaino] 580! ************************************************************************************************** 581 SUBROUTINE restraint_colv_low(colv_list, fixd_list, lcolv, & 582 particle_set, cell, energy, force) 583 584 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 585 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 586 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 587 TYPE(particle_type), POINTER :: particle_set(:) 588 TYPE(cell_type), POINTER :: cell 589 REAL(KIND=dp), INTENT(INOUT) :: energy 590 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 591 592 CHARACTER(LEN=*), PARAMETER :: routineN = 'restraint_colv_low', & 593 routineP = moduleN//':'//routineN 594 595 INTEGER :: iatm, iconst, ind 596 REAL(KIND=dp) :: k, tab, targ 597 598 DO iconst = 1, SIZE(colv_list) 599 IF (.NOT. colv_list(iconst)%restraint%active) CYCLE 600 ! Update colvar 601 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, & 602 particles=particle_set, fixd_list=fixd_list) 603 604 k = colv_list(iconst)%restraint%k0 605 targ = colv_list(iconst)%expected_value 606 tab = diff_colvar(lcolv(iconst)%colvar, targ) 607 ! Update Energy 608 energy = energy + k*tab**2 609 ! Update Forces 610 DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom) 611 ind = lcolv(iconst)%colvar%i_atom(iatm) 612 force(:, ind) = force(:, ind) - 2.0_dp*k*tab*lcolv(iconst)%colvar%dsdr(:, iatm) 613 END DO 614 END DO 615 616 END SUBROUTINE restraint_colv_low 617 618END MODULE restraint 619