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_fxd 11 12 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 13 USE atomic_kind_types, ONLY: get_atomic_kind_set 14 USE cell_types, ONLY: 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_types, ONLY: colvar_type 22 USE cp_subsys_types, ONLY: cp_subsys_get,& 23 cp_subsys_type 24 USE distribution_1d_types, ONLY: distribution_1d_type 25 USE force_env_types, ONLY: force_env_get,& 26 force_env_type 27 USE kinds, ONLY: dp 28 USE message_passing, ONLY: mp_sum 29 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 30 USE molecule_kind_types, ONLY: fixd_constraint_type,& 31 get_molecule_kind,& 32 local_fixd_constraint_type,& 33 molecule_kind_type 34 USE molecule_types, ONLY: local_g3x3_constraint_type,& 35 local_g4x6_constraint_type 36 USE particle_list_types, ONLY: particle_list_type 37 USE particle_types, ONLY: particle_type 38 USE util, ONLY: sort 39#include "./base/base_uses.f90" 40 41 IMPLICIT NONE 42 43 PRIVATE 44 PUBLIC :: fix_atom_control, & 45 check_fixed_atom_cns_g3x3, & 46 check_fixed_atom_cns_g4x6, & 47 check_fixed_atom_cns_colv, & 48 create_local_fixd_list, & 49 release_local_fixd_list 50 51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_fxd' 52 53CONTAINS 54 55! ************************************************************************************************** 56!> \brief allows for fix atom constraints 57!> \param force_env ... 58!> \param w ... 59!> \par History 60!> - optionally apply fix atom constraint to random forces (Langevin) 61!> (04.10.206,MK) 62! ************************************************************************************************** 63 SUBROUTINE fix_atom_control(force_env, w) 64 TYPE(force_env_type), POINTER :: force_env 65 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL :: w 66 67 CHARACTER(len=*), PARAMETER :: routineN = 'fix_atom_control', & 68 routineP = moduleN//':'//routineN 69 70 INTEGER :: handle, i, ifixd, ii, ikind, iparticle, iparticle_local, my_atm_fixed, natom, & 71 ncore, nfixed_atoms, nkind, nparticle, nparticle_local, nshell, shell_index 72 LOGICAL :: shell_present 73 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: force 74 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 75 TYPE(cp_subsys_type), POINTER :: subsys 76 TYPE(distribution_1d_type), POINTER :: local_particles 77 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 78 TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:) 79 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 80 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 81 TYPE(molecule_kind_type), POINTER :: molecule_kind 82 TYPE(particle_list_type), POINTER :: core_particles, particles, & 83 shell_particles 84 TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, particle_set, & 85 shell_particle_set 86 87 CALL timeset(routineN, handle) 88 89 NULLIFY (atomic_kinds) 90 NULLIFY (core_particles) 91 NULLIFY (particles) 92 NULLIFY (shell_particles) 93 shell_present = .FALSE. 94 95 NULLIFY (lfixd_list) 96 CALL force_env_get(force_env=force_env, & 97 subsys=subsys) 98 CALL cp_subsys_get(subsys=subsys, & 99 atomic_kinds=atomic_kinds, & 100 core_particles=core_particles, & 101 local_particles=local_particles, & 102 molecule_kinds=molecule_kinds, & 103 natom=natom, & 104 ncore=ncore, & 105 nshell=nshell, & 106 particles=particles, & 107 shell_particles=shell_particles) 108 CALL get_atomic_kind_set(atomic_kind_set=atomic_kinds%els, & 109 shell_present=shell_present) 110 111 particle_set => particles%els 112 CPASSERT((SIZE(particle_set) == natom)) 113 IF (shell_present) THEN 114 core_particle_set => core_particles%els 115 CPASSERT((SIZE(core_particle_set) == ncore)) 116 shell_particle_set => shell_particles%els 117 CPASSERT((SIZE(shell_particle_set) == nshell)) 118 END IF 119 nparticle = natom + nshell 120 molecule_kind_set => molecule_kinds%els 121 122 nkind = molecule_kinds%n_els 123 my_atm_fixed = 0 124 DO ikind = 1, nkind 125 molecule_kind => molecule_kind_set(ikind) 126 CALL get_molecule_kind(molecule_kind, nfixd=nfixed_atoms) 127 my_atm_fixed = my_atm_fixed + nfixed_atoms 128 END DO 129 130 IF (my_atm_fixed /= 0) THEN 131 IF (.NOT. PRESENT(w)) THEN 132 ! Allocate scratch array 133 ALLOCATE (force(3, nparticle)) 134 force(:, :) = 0.0_dp 135 DO i = 1, SIZE(local_particles%n_el) 136 nparticle_local = local_particles%n_el(i) 137 DO iparticle_local = 1, nparticle_local 138 iparticle = local_particles%list(i)%array(iparticle_local) 139 shell_index = particle_set(iparticle)%shell_index 140 IF (shell_index == 0) THEN 141 force(:, iparticle) = particle_set(iparticle)%f(:) 142 ELSE 143 force(:, iparticle) = core_particle_set(shell_index)%f(:) 144 force(:, natom + shell_index) = shell_particle_set(shell_index)%f(:) 145 END IF 146 END DO 147 END DO 148 END IF 149 150 ! Create the list of locally fixed atoms 151 CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles) 152 153 ! Apply fixed atom constraint 154 DO ifixd = 1, SIZE(lfixd_list) 155 ikind = lfixd_list(ifixd)%ikind 156 ii = lfixd_list(ifixd)%ifixd_index 157 molecule_kind => molecule_kind_set(ikind) 158 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list) 159 IF (.NOT. fixd_list(ii)%restraint%active) THEN 160 iparticle = fixd_list(ii)%fixd 161 shell_index = particle_set(iparticle)%shell_index 162 ! Select constraint type 163 IF (PRESENT(w)) THEN 164 SELECT CASE (fixd_list(ii)%itype) 165 CASE (use_perd_x) 166 w(1, iparticle) = 0.0_dp 167 CASE (use_perd_y) 168 w(2, iparticle) = 0.0_dp 169 CASE (use_perd_z) 170 w(3, iparticle) = 0.0_dp 171 CASE (use_perd_xy) 172 w(1, iparticle) = 0.0_dp 173 w(2, iparticle) = 0.0_dp 174 CASE (use_perd_xz) 175 w(1, iparticle) = 0.0_dp 176 w(3, iparticle) = 0.0_dp 177 CASE (use_perd_yz) 178 w(2, iparticle) = 0.0_dp 179 w(3, iparticle) = 0.0_dp 180 CASE (use_perd_xyz) 181 w(:, iparticle) = 0.0_dp 182 END SELECT 183 ELSE 184 SELECT CASE (fixd_list(ii)%itype) 185 CASE (use_perd_x) 186 force(1, iparticle) = 0.0_dp 187 IF (shell_index /= 0) THEN 188 force(1, natom + shell_index) = 0.0_dp 189 END IF 190 CASE (use_perd_y) 191 force(2, iparticle) = 0.0_dp 192 IF (shell_index /= 0) THEN 193 force(2, natom + shell_index) = 0.0_dp 194 END IF 195 CASE (use_perd_z) 196 force(3, iparticle) = 0.0_dp 197 IF (shell_index /= 0) THEN 198 force(3, natom + shell_index) = 0.0_dp 199 END IF 200 CASE (use_perd_xy) 201 force(1, iparticle) = 0.0_dp 202 force(2, iparticle) = 0.0_dp 203 IF (shell_index /= 0) THEN 204 force(1, natom + shell_index) = 0.0_dp 205 force(2, natom + shell_index) = 0.0_dp 206 END IF 207 CASE (use_perd_xz) 208 force(1, iparticle) = 0.0_dp 209 force(3, iparticle) = 0.0_dp 210 IF (shell_index /= 0) THEN 211 force(1, natom + shell_index) = 0.0_dp 212 force(3, natom + shell_index) = 0.0_dp 213 END IF 214 CASE (use_perd_yz) 215 force(2, iparticle) = 0.0_dp 216 force(3, iparticle) = 0.0_dp 217 IF (shell_index /= 0) THEN 218 force(2, natom + shell_index) = 0.0_dp 219 force(3, natom + shell_index) = 0.0_dp 220 END IF 221 CASE (use_perd_xyz) 222 force(:, iparticle) = 0.0_dp 223 IF (shell_index /= 0) THEN 224 force(:, natom + shell_index) = 0.0_dp 225 END IF 226 END SELECT 227 END IF 228 END IF 229 END DO 230 CALL release_local_fixd_list(lfixd_list) 231 232 IF (.NOT. PRESENT(w)) THEN 233 CALL mp_sum(force, force_env%para_env%group) 234 DO iparticle = 1, natom 235 shell_index = particle_set(iparticle)%shell_index 236 IF (shell_index == 0) THEN 237 particle_set(iparticle)%f(:) = force(:, iparticle) 238 ELSE 239 core_particle_set(shell_index)%f(:) = force(:, iparticle) 240 shell_particle_set(shell_index)%f(:) = force(:, natom + shell_index) 241 END IF 242 END DO 243 DEALLOCATE (force) 244 END IF 245 END IF 246 247 CALL timestop(handle) 248 249 END SUBROUTINE fix_atom_control 250 251! ************************************************************************************************** 252!> \brief ... 253!> \param imass1 ... 254!> \param imass2 ... 255!> \param imass3 ... 256!> \param index_a ... 257!> \param index_b ... 258!> \param index_c ... 259!> \param fixd_list ... 260!> \param lg3x3 ... 261!> \par History 262!> none 263! ************************************************************************************************** 264 SUBROUTINE check_fixed_atom_cns_g3x3(imass1, imass2, imass3, & 265 index_a, index_b, index_c, fixd_list, lg3x3) 266 REAL(KIND=dp), INTENT(INOUT) :: imass1, imass2, imass3 267 INTEGER, INTENT(IN) :: index_a, index_b, index_c 268 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 269 TYPE(local_g3x3_constraint_type) :: lg3x3 270 271 INTEGER :: i 272 273 IF (lg3x3%init) THEN 274 imass1 = lg3x3%imass1 275 imass2 = lg3x3%imass2 276 imass3 = lg3x3%imass3 277 ELSE 278 IF (ASSOCIATED(fixd_list)) THEN 279 IF (SIZE(fixd_list) > 0) THEN 280 DO i = 1, SIZE(fixd_list) 281 IF (fixd_list(i)%fixd == index_a) THEN 282 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 283 IF (.NOT. fixd_list(i)%restraint%active) imass1 = 0.0_dp 284 EXIT 285 END IF 286 END DO 287 DO i = 1, SIZE(fixd_list) 288 IF (fixd_list(i)%fixd == index_b) THEN 289 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 290 IF (.NOT. fixd_list(i)%restraint%active) imass2 = 0.0_dp 291 EXIT 292 END IF 293 END DO 294 DO i = 1, SIZE(fixd_list) 295 IF (fixd_list(i)%fixd == index_c) THEN 296 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 297 IF (.NOT. fixd_list(i)%restraint%active) imass3 = 0.0_dp 298 EXIT 299 END IF 300 END DO 301 END IF 302 END IF 303 lg3x3%imass1 = imass1 304 lg3x3%imass2 = imass2 305 lg3x3%imass3 = imass3 306 lg3x3%init = .TRUE. 307 END IF 308 END SUBROUTINE check_fixed_atom_cns_g3x3 309 310! ************************************************************************************************** 311!> \brief ... 312!> \param imass1 ... 313!> \param imass2 ... 314!> \param imass3 ... 315!> \param imass4 ... 316!> \param index_a ... 317!> \param index_b ... 318!> \param index_c ... 319!> \param index_d ... 320!> \param fixd_list ... 321!> \param lg4x6 ... 322!> \par History 323!> none 324! ************************************************************************************************** 325 SUBROUTINE check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, & 326 index_a, index_b, index_c, index_d, fixd_list, lg4x6) 327 REAL(KIND=dp), INTENT(INOUT) :: imass1, imass2, imass3, imass4 328 INTEGER, INTENT(IN) :: index_a, index_b, index_c, index_d 329 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 330 TYPE(local_g4x6_constraint_type) :: lg4x6 331 332 INTEGER :: i 333 334 IF (lg4x6%init) THEN 335 imass1 = lg4x6%imass1 336 imass2 = lg4x6%imass2 337 imass3 = lg4x6%imass3 338 imass4 = lg4x6%imass4 339 ELSE 340 IF (ASSOCIATED(fixd_list)) THEN 341 IF (SIZE(fixd_list) > 0) THEN 342 DO i = 1, SIZE(fixd_list) 343 IF (fixd_list(i)%fixd == index_a) THEN 344 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 345 IF (.NOT. fixd_list(i)%restraint%active) imass1 = 0.0_dp 346 EXIT 347 END IF 348 END DO 349 DO i = 1, SIZE(fixd_list) 350 IF (fixd_list(i)%fixd == index_b) THEN 351 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 352 IF (.NOT. fixd_list(i)%restraint%active) imass2 = 0.0_dp 353 EXIT 354 END IF 355 END DO 356 DO i = 1, SIZE(fixd_list) 357 IF (fixd_list(i)%fixd == index_c) THEN 358 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 359 IF (.NOT. fixd_list(i)%restraint%active) imass3 = 0.0_dp 360 EXIT 361 END IF 362 END DO 363 DO i = 1, SIZE(fixd_list) 364 IF (fixd_list(i)%fixd == index_d) THEN 365 IF (fixd_list(i)%itype /= use_perd_xyz) CYCLE 366 IF (.NOT. fixd_list(i)%restraint%active) imass4 = 0.0_dp 367 EXIT 368 END IF 369 END DO 370 END IF 371 END IF 372 lg4x6%imass1 = imass1 373 lg4x6%imass2 = imass2 374 lg4x6%imass3 = imass3 375 lg4x6%imass4 = imass4 376 lg4x6%init = .TRUE. 377 END IF 378 END SUBROUTINE check_fixed_atom_cns_g4x6 379 380! ************************************************************************************************** 381!> \brief ... 382!> \param fixd_list ... 383!> \param colvar ... 384!> \par History 385!> none 386! ************************************************************************************************** 387 SUBROUTINE check_fixed_atom_cns_colv(fixd_list, colvar) 388 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 389 TYPE(colvar_type), POINTER :: colvar 390 391 INTEGER :: i, j, k 392 393 IF (ASSOCIATED(fixd_list)) THEN 394 IF (ASSOCIATED(fixd_list)) THEN 395 IF (SIZE(fixd_list) > 0) THEN 396 DO i = 1, SIZE(colvar%i_atom) 397 j = colvar%i_atom(i) 398 DO k = 1, SIZE(fixd_list) 399 IF (fixd_list(k)%fixd == j) THEN 400 IF (fixd_list(k)%itype /= use_perd_xyz) CYCLE 401 IF (.NOT. fixd_list(k)%restraint%active) & 402 colvar%dsdr(:, i) = 0.0_dp 403 EXIT 404 END IF 405 END DO 406 END DO 407 END IF 408 END IF 409 END IF 410 411 END SUBROUTINE check_fixed_atom_cns_colv 412 413! ************************************************************************************************** 414!> \brief setup a list of local atoms on which to apply constraints/restraints 415!> \param lfixd_list ... 416!> \param nkind ... 417!> \param molecule_kind_set ... 418!> \param local_particles ... 419!> \author Teodoro Laino [tlaino] - 11.2008 420! ************************************************************************************************** 421 SUBROUTINE create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, & 422 local_particles) 423 TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:) 424 INTEGER, INTENT(IN) :: nkind 425 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 426 TYPE(distribution_1d_type), POINTER :: local_particles 427 428 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_local_fixd_list', & 429 routineP = moduleN//':'//routineN 430 431 INTEGER :: handle, i, ikind, iparticle, & 432 iparticle_local, isize, jsize, ncnst, & 433 nparticle_local, nparticle_local_all, & 434 nsize 435 INTEGER, ALLOCATABLE, DIMENSION(:) :: fixed_atom_all, kind_index_all, & 436 local_particle_all, work0, work1, work2 437 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 438 TYPE(molecule_kind_type), POINTER :: molecule_kind 439 440 CALL timeset(routineN, handle) 441 CPASSERT(.NOT. ASSOCIATED(lfixd_list)) 442 nsize = 0 443 DO ikind = 1, nkind 444 molecule_kind => molecule_kind_set(ikind) 445 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list) 446 IF (ASSOCIATED(fixd_list)) THEN 447 nsize = nsize + SIZE(fixd_list) 448 END IF 449 END DO 450 IF (nsize /= 0) THEN 451 ALLOCATE (fixed_atom_all(nsize)) 452 ALLOCATE (work0(nsize)) 453 ALLOCATE (work1(nsize)) 454 ALLOCATE (kind_index_all(nsize)) 455 nsize = 0 456 DO ikind = 1, nkind 457 molecule_kind => molecule_kind_set(ikind) 458 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list) 459 IF (ASSOCIATED(fixd_list)) THEN 460 DO i = 1, SIZE(fixd_list) 461 nsize = nsize + 1 462 work0(nsize) = i 463 kind_index_all(nsize) = ikind 464 fixed_atom_all(nsize) = fixd_list(i)%fixd 465 END DO 466 END IF 467 END DO 468 ! Sort the number of all atoms to be constrained/restrained 469 CALL sort(fixed_atom_all, nsize, work1) 470 471 ! Sort the local particles 472 nparticle_local_all = 0 473 DO i = 1, SIZE(local_particles%n_el) 474 nparticle_local_all = nparticle_local_all + local_particles%n_el(i) 475 END DO 476 ALLOCATE (local_particle_all(nparticle_local_all)) 477 ALLOCATE (work2(nparticle_local_all)) 478 nparticle_local_all = 0 479 DO i = 1, SIZE(local_particles%n_el) 480 nparticle_local = local_particles%n_el(i) 481 DO iparticle_local = 1, nparticle_local 482 nparticle_local_all = nparticle_local_all + 1 483 iparticle = local_particles%list(i)%array(iparticle_local) 484 local_particle_all(nparticle_local_all) = iparticle 485 END DO 486 END DO 487 CALL sort(local_particle_all, nparticle_local_all, work2) 488 489 ! Count the amount of local constraints/restraints 490 ncnst = 0 491 jsize = 1 492 Loop_count: DO isize = 1, nparticle_local_all 493 DO WHILE (local_particle_all(isize) > fixed_atom_all(jsize)) 494 jsize = jsize + 1 495 IF (jsize > nsize) THEN 496 jsize = nsize 497 EXIT Loop_count 498 END IF 499 END DO 500 IF (local_particle_all(isize) == fixed_atom_all(jsize)) ncnst = ncnst + 1 501 END DO Loop_count 502 503 ! Allocate local fixed atom array 504 ALLOCATE (lfixd_list(ncnst)) 505 506 ! Fill array with constraints infos 507 ncnst = 0 508 jsize = 1 509 Loop_fill: DO isize = 1, nparticle_local_all 510 DO WHILE (local_particle_all(isize) > fixed_atom_all(jsize)) 511 jsize = jsize + 1 512 IF (jsize > nsize) THEN 513 jsize = nsize 514 EXIT Loop_fill 515 END IF 516 END DO 517 IF (local_particle_all(isize) == fixed_atom_all(jsize)) THEN 518 ncnst = ncnst + 1 519 lfixd_list(ncnst)%ifixd_index = work0(work1(jsize)) 520 lfixd_list(ncnst)%ikind = kind_index_all(work1(jsize)) 521 END IF 522 END DO Loop_fill 523 524 ! Deallocate working arrays 525 DEALLOCATE (local_particle_all) 526 DEALLOCATE (work2) 527 DEALLOCATE (fixed_atom_all) 528 DEALLOCATE (work1) 529 DEALLOCATE (kind_index_all) 530 ELSE 531 ! Allocate local fixed atom array with dimension 0 532 ALLOCATE (lfixd_list(0)) 533 END IF 534 CALL timestop(handle) 535 END SUBROUTINE create_local_fixd_list 536 537! ************************************************************************************************** 538!> \brief destroy the list of local atoms on which to apply constraints/restraints 539!> Teodoro Laino [tlaino] - 11.2008 540!> \param lfixd_list ... 541! ************************************************************************************************** 542 SUBROUTINE release_local_fixd_list(lfixd_list) 543 TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:) 544 545 CHARACTER(LEN=*), PARAMETER :: routineN = 'release_local_fixd_list', & 546 routineP = moduleN//':'//routineN 547 548 CPASSERT(ASSOCIATED(lfixd_list)) 549 DEALLOCATE (lfixd_list) 550 END SUBROUTINE release_local_fixd_list 551 552END MODULE constraint_fxd 553