1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief evaluations of colvar for internal coordinates schemes 8!> \par History 9!> 05-2007 created [tlaino] 10!> \author Teodoro Laino - Zurich University (2007) [tlaino] 11! ************************************************************************************************** 12MODULE colvar_utils 13 USE cell_types, ONLY: cell_type 14 USE colvar_methods, ONLY: colvar_eval_mol_f 15 USE colvar_types, ONLY: & 16 colvar_counters, colvar_setup, colvar_type, coord_colvar_id, distance_from_path_colvar_id, & 17 gyration_colvar_id, mindist_colvar_id, population_colvar_id, reaction_path_colvar_id, & 18 rmsd_colvar_id 19 USE cp_subsys_types, ONLY: cp_subsys_get,& 20 cp_subsys_type 21 USE distribution_1d_types, ONLY: distribution_1d_type 22 USE force_env_types, ONLY: force_env_get,& 23 force_env_type 24 USE input_constants, ONLY: rmsd_all,& 25 rmsd_list 26 USE kinds, ONLY: dp 27 USE mathlib, ONLY: invert_matrix 28 USE memory_utilities, ONLY: reallocate 29 USE message_passing, ONLY: mp_sum 30 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 31 USE molecule_kind_types, ONLY: colvar_constraint_type,& 32 fixd_constraint_type,& 33 get_molecule_kind,& 34 molecule_kind_type 35 USE molecule_list_types, ONLY: molecule_list_type 36 USE molecule_types, ONLY: get_molecule,& 37 global_constraint_type,& 38 local_colvar_constraint_type,& 39 molecule_type 40 USE particle_list_types, ONLY: particle_list_type 41 USE particle_types, ONLY: particle_type 42 USE rmsd, ONLY: rmsd3 43 USE string_utilities, ONLY: uppercase 44 USE util, ONLY: sort 45#include "./base/base_uses.f90" 46 47 IMPLICIT NONE 48 49 PRIVATE 50 PUBLIC :: number_of_colvar, & 51 eval_colvar, & 52 set_colvars_target, & 53 get_clv_force, & 54 post_process_colvar 55 56 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'colvar_utils' 57 58CONTAINS 59 60! ************************************************************************************************** 61!> \brief Gives back the number of colvar defined for a force_eval 62!> \param force_env ... 63!> \param only_intra_colvar ... 64!> \param unique ... 65!> \return ... 66!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 67! ************************************************************************************************** 68 FUNCTION number_of_colvar(force_env, only_intra_colvar, unique) RESULT(ntot) 69 TYPE(force_env_type), POINTER :: force_env 70 LOGICAL, INTENT(IN), OPTIONAL :: only_intra_colvar, unique 71 INTEGER :: ntot 72 73 CHARACTER(LEN=*), PARAMETER :: routineN = 'number_of_colvar', & 74 routineP = moduleN//':'//routineN 75 76 INTEGER :: handle, ikind, imol 77 LOGICAL :: my_unique, skip_inter_colvar 78 TYPE(colvar_counters) :: ncolv 79 TYPE(cp_subsys_type), POINTER :: subsys 80 TYPE(global_constraint_type), POINTER :: gci 81 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 82 TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:) 83 TYPE(molecule_list_type), POINTER :: molecules 84 TYPE(molecule_type), POINTER :: molecule, molecule_set(:) 85 86 NULLIFY (subsys, molecules, molecule_kind, molecule, molecule_set, gci) 87 CALL timeset(routineN, handle) 88 skip_inter_colvar = .FALSE. 89 my_unique = .FALSE. 90 IF (PRESENT(only_intra_colvar)) skip_inter_colvar = only_intra_colvar 91 IF (PRESENT(unique)) my_unique = unique 92 ntot = 0 93 CALL force_env_get(force_env=force_env, subsys=subsys) 94 CALL cp_subsys_get(subsys=subsys, molecules=molecules, gci=gci, & 95 molecule_kinds=molecule_kinds) 96 97 molecule_set => molecules%els 98 ! Intramolecular Colvar 99 IF (my_unique) THEN 100 molecule_kind_set => molecule_kinds%els 101 DO ikind = 1, molecule_kinds%n_els 102 molecule_kind => molecule_kind_set(ikind) 103 CALL get_molecule_kind(molecule_kind, ncolv=ncolv) 104 ntot = ntot + ncolv%ntot 105 END DO 106 ELSE 107 MOL: DO imol = 1, SIZE(molecule_set) 108 molecule => molecule_set(imol) 109 molecule_kind => molecule%molecule_kind 110 111 CALL get_molecule_kind(molecule_kind, & 112 ncolv=ncolv) 113 ntot = ntot + ncolv%ntot 114 END DO MOL 115 END IF 116 ! Intermolecular Colvar 117 IF (.NOT. skip_inter_colvar) THEN 118 IF (ASSOCIATED(gci)) THEN 119 ntot = ntot + gci%ncolv%ntot 120 END IF 121 END IF 122 CALL timestop(handle) 123 124 END FUNCTION number_of_colvar 125 126! ************************************************************************************************** 127!> \brief Set the value of target for constraints/restraints 128!> \param targets ... 129!> \param force_env ... 130!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 131! ************************************************************************************************** 132 SUBROUTINE set_colvars_target(targets, force_env) 133 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: targets 134 TYPE(force_env_type), POINTER :: force_env 135 136 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_colvars_target', & 137 routineP = moduleN//':'//routineN 138 139 INTEGER :: handle, i, ikind, ind, nkind 140 TYPE(cell_type), POINTER :: cell 141 TYPE(colvar_constraint_type), DIMENSION(:), & 142 POINTER :: colv_list 143 TYPE(colvar_counters) :: ncolv 144 TYPE(cp_subsys_type), POINTER :: subsys 145 TYPE(global_constraint_type), POINTER :: gci 146 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 147 TYPE(molecule_kind_type), POINTER :: molecule_kind 148 149 NULLIFY (cell, subsys, molecule_kinds, molecule_kind, gci, colv_list) 150 CALL timeset(routineN, handle) 151 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 152 CALL cp_subsys_get(subsys=subsys, gci=gci, molecule_kinds=molecule_kinds) 153 154 nkind = molecule_kinds%n_els 155 ! Set Target for Intramolecular Colvars 156 MOL: DO ikind = 1, nkind 157 molecule_kind => molecule_kinds%els(ikind) 158 CALL get_molecule_kind(molecule_kind, & 159 colv_list=colv_list, & 160 ncolv=ncolv) 161 IF (ncolv%ntot /= 0) THEN 162 DO i = 1, SIZE(colv_list) 163 ind = colv_list(i)%inp_seq_num 164 colv_list(i)%expected_value = targets(ind) 165 END DO 166 END IF 167 END DO MOL 168 ! Set Target for Intermolecular Colvars 169 IF (ASSOCIATED(gci)) THEN 170 IF (gci%ncolv%ntot /= 0) THEN 171 colv_list => gci%colv_list 172 DO i = 1, SIZE(colv_list) 173 ind = colv_list(i)%inp_seq_num 174 colv_list(i)%expected_value = targets(ind) 175 END DO 176 ENDIF 177 END IF 178 CALL timestop(handle) 179 180 END SUBROUTINE set_colvars_target 181 182! ************************************************************************************************** 183!> \brief Computes the values of colvars and the Wilson matrix B and its invers A 184!> \param force_env ... 185!> \param coords ... 186!> \param cvalues ... 187!> \param Bmatrix ... 188!> \param MassI ... 189!> \param Amatrix ... 190!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 191! ************************************************************************************************** 192 SUBROUTINE eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix) 193 194 TYPE(force_env_type), POINTER :: force_env 195 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords 196 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cvalues 197 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: Bmatrix 198 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: MassI 199 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: Amatrix 200 201 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_colvar', routineP = moduleN//':'//routineN 202 203 INTEGER :: handle, i, ikind, imol, n_tot, natom, & 204 nkind, nmol_per_kind, offset 205 INTEGER, DIMENSION(:), POINTER :: map, wrk 206 LOGICAL :: check 207 REAL(KIND=dp) :: inv_error 208 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bwrk, Gmatrix, Gmatrix_i 209 REAL(KIND=dp), DIMENSION(:), POINTER :: rwrk 210 TYPE(cell_type), POINTER :: cell 211 TYPE(colvar_counters) :: ncolv 212 TYPE(cp_subsys_type), POINTER :: subsys 213 TYPE(distribution_1d_type), POINTER :: local_molecules 214 TYPE(global_constraint_type), POINTER :: gci 215 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 216 TYPE(molecule_kind_type), POINTER :: molecule_kind 217 TYPE(molecule_list_type), POINTER :: molecules 218 TYPE(molecule_type), POINTER :: molecule, molecule_set(:) 219 TYPE(particle_list_type), POINTER :: particles 220 TYPE(particle_type), POINTER :: particle_set(:) 221 222 NULLIFY (cell, subsys, local_molecules, molecule_kinds, & 223 molecules, molecule_kind, molecule, & 224 molecule_set, particles, particle_set, gci) 225 IF (PRESENT(Bmatrix)) THEN 226 check = ASSOCIATED(Bmatrix) 227 CPASSERT(check) 228 Bmatrix = 0.0_dp 229 END IF 230 CALL timeset(routineN, handle) 231 ALLOCATE (map(SIZE(cvalues))) 232 map = HUGE(0) ! init all, since used in a sort, but not all set in parallel. 233 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell) 234 n_tot = 0 235 cvalues = 0.0_dp 236 CALL cp_subsys_get(subsys=subsys, & 237 particles=particles, & 238 molecules=molecules, & 239 local_molecules=local_molecules, & 240 gci=gci, & 241 molecule_kinds=molecule_kinds) 242 243 nkind = molecule_kinds%n_els 244 particle_set => particles%els 245 molecule_set => molecules%els 246 ! Intramolecular Colvars 247 IF (number_of_colvar(force_env, only_intra_colvar=.TRUE.) /= 0) THEN 248 MOL: DO ikind = 1, nkind 249 nmol_per_kind = local_molecules%n_el(ikind) 250 DO imol = 1, nmol_per_kind 251 i = local_molecules%list(ikind)%array(imol) 252 molecule => molecule_set(i) 253 molecule_kind => molecule%molecule_kind 254 255 CALL get_molecule_kind(molecule_kind, & 256 ncolv=ncolv) 257 offset = get_colvar_offset(i, molecule_set) 258 ! Collective variables 259 IF (ncolv%ntot /= 0) THEN 260 CALL eval_colv_int(molecule, particle_set, coords, cell, cvalues, & 261 Bmatrix, offset, n_tot, map) 262 ENDIF 263 END DO 264 END DO MOL 265 CALL mp_sum(n_tot, force_env%para_env%group) 266 CALL mp_sum(cvalues, force_env%para_env%group) 267 IF (PRESENT(Bmatrix)) CALL mp_sum(Bmatrix, force_env%para_env%group) 268 END IF 269 offset = n_tot 270 ! Intermolecular Colvars 271 IF (ASSOCIATED(gci)) THEN 272 IF (gci%ncolv%ntot /= 0) THEN 273 CALL eval_colv_ext(gci, particle_set, coords, cell, cvalues, & 274 Bmatrix, offset, n_tot, map) 275 ENDIF 276 END IF 277 CPASSERT(n_tot == SIZE(cvalues)) 278 ! Sort values of Collective Variables according the order of the input 279 ! sections 280 ALLOCATE (wrk(SIZE(cvalues))) 281 ALLOCATE (rwrk(SIZE(cvalues))) 282 CALL sort(map, SIZE(map), wrk) 283 rwrk = cvalues 284 DO i = 1, SIZE(wrk) 285 cvalues(i) = rwrk(wrk(i)) 286 END DO 287 ! check and sort on Bmatrix 288 IF (PRESENT(Bmatrix)) THEN 289 check = n_tot == SIZE(Bmatrix, 2) 290 CPASSERT(check) 291 ALLOCATE (bwrk(SIZE(Bmatrix, 1), SIZE(Bmatrix, 2))) 292 bwrk(:, :) = Bmatrix 293 DO i = 1, SIZE(wrk) 294 Bmatrix(:, i) = bwrk(:, wrk(i)) 295 END DO 296 DEALLOCATE (bwrk) 297 END IF 298 DEALLOCATE (rwrk) 299 DEALLOCATE (wrk) 300 DEALLOCATE (map) 301 ! Construction of the Amatrix 302 IF (PRESENT(Bmatrix) .AND. PRESENT(Amatrix)) THEN 303 CPASSERT(ASSOCIATED(Amatrix)) 304 check = SIZE(Bmatrix, 1) == SIZE(Amatrix, 2) 305 CPASSERT(check) 306 check = SIZE(Bmatrix, 2) == SIZE(Amatrix, 1) 307 CPASSERT(check) 308 ALLOCATE (Gmatrix(n_tot, n_tot)) 309 ALLOCATE (Gmatrix_i(n_tot, n_tot)) 310 Gmatrix(:, :) = MATMUL(TRANSPOSE(Bmatrix), Bmatrix) 311 CALL invert_matrix(Gmatrix, Gmatrix_i, inv_error) 312 IF (ABS(inv_error) > 1.0E-8_dp) & 313 CPWARN("Error in inverting the Gmatrix larger than 1.0E-8!") 314 Amatrix = MATMUL(Gmatrix_i, TRANSPOSE(Bmatrix)) 315 DEALLOCATE (Gmatrix_i) 316 DEALLOCATE (Gmatrix) 317 END IF 318 IF (PRESENT(MassI)) THEN 319 natom = SIZE(particle_set) 320 CPASSERT(ASSOCIATED(MassI)) 321 CPASSERT(SIZE(MassI) == natom*3) 322 DO i = 1, natom 323 MassI((i - 1)*3 + 1) = 1.0_dp/particle_set(i)%atomic_kind%mass 324 MassI((i - 1)*3 + 2) = 1.0_dp/particle_set(i)%atomic_kind%mass 325 MassI((i - 1)*3 + 3) = 1.0_dp/particle_set(i)%atomic_kind%mass 326 END DO 327 END IF 328 CALL timestop(handle) 329 330 END SUBROUTINE eval_colvar 331 332! ************************************************************************************************** 333!> \brief Computes the offset of the colvar for the specific molecule 334!> \param i ... 335!> \param molecule_set ... 336!> \return ... 337!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 338! ************************************************************************************************** 339 FUNCTION get_colvar_offset(i, molecule_set) RESULT(offset) 340 INTEGER, INTENT(IN) :: i 341 TYPE(molecule_type), POINTER :: molecule_set(:) 342 INTEGER :: offset 343 344 INTEGER :: j 345 TYPE(colvar_counters) :: ncolv 346 TYPE(molecule_kind_type), POINTER :: molecule_kind 347 TYPE(molecule_type), POINTER :: molecule 348 349 offset = 0 350 DO j = 1, i - 1 351 molecule => molecule_set(j) 352 molecule_kind => molecule%molecule_kind 353 CALL get_molecule_kind(molecule_kind, & 354 ncolv=ncolv) 355 offset = offset + ncolv%ntot 356 END DO 357 358 END FUNCTION get_colvar_offset 359 360! ************************************************************************************************** 361!> \brief Computes Intramolecular colvar 362!> \param molecule ... 363!> \param particle_set ... 364!> \param coords ... 365!> \param cell ... 366!> \param cvalues ... 367!> \param Bmatrix ... 368!> \param offset ... 369!> \param n_tot ... 370!> \param map ... 371!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 372! ************************************************************************************************** 373 SUBROUTINE eval_colv_int(molecule, particle_set, coords, cell, cvalues, & 374 Bmatrix, offset, n_tot, map) 375 376 TYPE(molecule_type), POINTER :: molecule 377 TYPE(particle_type), POINTER :: particle_set(:) 378 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords 379 TYPE(cell_type), POINTER :: cell 380 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cvalues 381 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: Bmatrix 382 INTEGER, INTENT(IN) :: offset 383 INTEGER, INTENT(INOUT) :: n_tot 384 INTEGER, DIMENSION(:), POINTER :: map 385 386 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_colv_int', routineP = moduleN//':'//routineN 387 388 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 389 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 390 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 391 TYPE(molecule_kind_type), POINTER :: molecule_kind 392 393 NULLIFY (fixd_list) 394 395 molecule_kind => molecule%molecule_kind 396 CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list) 397 CALL get_molecule(molecule, lcolv=lcolv) 398 CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, & 399 coords, cell, cvalues, Bmatrix, offset, n_tot, map) 400 401 END SUBROUTINE eval_colv_int 402 403! ************************************************************************************************** 404!> \brief Computes Intermolecular colvar 405!> \param gci ... 406!> \param particle_set ... 407!> \param coords ... 408!> \param cell ... 409!> \param cvalues ... 410!> \param Bmatrix ... 411!> \param offset ... 412!> \param n_tot ... 413!> \param map ... 414!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 415! ************************************************************************************************** 416 SUBROUTINE eval_colv_ext(gci, particle_set, coords, cell, cvalues, & 417 Bmatrix, offset, n_tot, map) 418 TYPE(global_constraint_type), POINTER :: gci 419 TYPE(particle_type), POINTER :: particle_set(:) 420 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords 421 TYPE(cell_type), POINTER :: cell 422 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cvalues 423 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: Bmatrix 424 INTEGER, INTENT(IN) :: offset 425 INTEGER, INTENT(INOUT) :: n_tot 426 INTEGER, DIMENSION(:), POINTER :: map 427 428 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_colv_ext', routineP = moduleN//':'//routineN 429 430 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 431 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 432 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 433 434 colv_list => gci%colv_list 435 fixd_list => gci%fixd_list 436 lcolv => gci%lcolv 437 CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, & 438 coords, cell, cvalues, Bmatrix, offset, n_tot, map) 439 440 END SUBROUTINE eval_colv_ext 441 442! ************************************************************************************************** 443!> \brief Real evaluation of colvar and of the Wilson-Eliashevich Matrix 444!> B_ik : i: internal coordinates 445!> k: cartesian coordinates 446!> \param colv_list ... 447!> \param fixd_list ... 448!> \param lcolv ... 449!> \param particle_set ... 450!> \param coords ... 451!> \param cell ... 452!> \param cvalues ... 453!> \param Bmatrix ... 454!> \param offset ... 455!> \param n_tot ... 456!> \param map ... 457!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University 458! ************************************************************************************************** 459 SUBROUTINE eval_colv_low(colv_list, fixd_list, lcolv, particle_set, coords, & 460 cell, cvalues, Bmatrix, offset, n_tot, map) 461 462 TYPE(colvar_constraint_type), POINTER :: colv_list(:) 463 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list 464 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:) 465 TYPE(particle_type), POINTER :: particle_set(:) 466 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords 467 TYPE(cell_type), POINTER :: cell 468 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cvalues 469 REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: Bmatrix 470 INTEGER, INTENT(IN) :: offset 471 INTEGER, INTENT(INOUT) :: n_tot 472 INTEGER, DIMENSION(:), POINTER :: map 473 474 CHARACTER(LEN=*), PARAMETER :: routineN = 'eval_colv_low', routineP = moduleN//':'//routineN 475 476 INTEGER :: iatm, iconst, ind, ival 477 478 ival = offset 479 DO iconst = 1, SIZE(colv_list) 480 n_tot = n_tot + 1 481 ival = ival + 1 482 ! Update colvar 483 IF (PRESENT(coords)) THEN 484 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, & 485 pos=RESHAPE(coords, (/3, SIZE(particle_set)/)), fixd_list=fixd_list) 486 ELSE 487 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, & 488 fixd_list=fixd_list) 489 END IF 490 cvalues(ival) = lcolv(iconst)%colvar%ss 491 map(ival) = colv_list(iconst)%inp_seq_num 492 ! Build the Wilson-Eliashevich Matrix 493 IF (PRESENT(Bmatrix)) THEN 494 DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom) 495 ind = (lcolv(iconst)%colvar%i_atom(iatm) - 1)*3 496 Bmatrix(ind + 1, ival) = lcolv(iconst)%colvar%dsdr(1, iatm) 497 Bmatrix(ind + 2, ival) = lcolv(iconst)%colvar%dsdr(2, iatm) 498 Bmatrix(ind + 3, ival) = lcolv(iconst)%colvar%dsdr(3, iatm) 499 END DO 500 END IF 501 END DO 502 503 END SUBROUTINE eval_colv_low 504 505! ************************************************************************************************** 506!> \brief Computes the forces in the frame of collective variables, and additional 507!> also the local metric tensor 508!> \param force_env ... 509!> \param forces ... 510!> \param coords ... 511!> \param nsize_xyz ... 512!> \param nsize_int ... 513!> \param cvalues ... 514!> \param Mmatrix ... 515!> \author Teodoro Laino 05.2007 516! ************************************************************************************************** 517 SUBROUTINE get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, & 518 Mmatrix) 519 TYPE(force_env_type), POINTER :: force_env 520 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & 521 OPTIONAL :: forces, coords 522 INTEGER, INTENT(IN) :: nsize_xyz, nsize_int 523 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cvalues, Mmatrix 524 525 CHARACTER(len=*), PARAMETER :: routineN = 'get_clv_force', routineP = moduleN//':'//routineN 526 527 INTEGER :: i, j, k 528 REAL(KIND=dp) :: tmp 529 REAL(KIND=dp), DIMENSION(:), POINTER :: MassI, wrk 530 REAL(KIND=dp), DIMENSION(:, :), POINTER :: Amatrix, Bmatrix 531 532 ALLOCATE (Bmatrix(nsize_xyz, nsize_int)) 533 ALLOCATE (MassI(nsize_xyz)) 534 ! Transform gradients if requested 535 IF (PRESENT(forces)) THEN 536 ALLOCATE (wrk(nsize_int)) 537 ALLOCATE (Amatrix(nsize_int, nsize_xyz)) 538 ! Compute the transformation matrices and the inverse mass diagonal Matrix 539 CALL eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix) 540 wrk = MATMUL(Amatrix, forces) 541 forces = 0.0_dp 542 forces(1:nsize_int) = wrk 543 DEALLOCATE (Amatrix) 544 DEALLOCATE (wrk) 545 ELSE 546 ! Compute the transformation matrices and the inverse mass diagonal Matrix 547 CALL eval_colvar(force_env, coords, cvalues, Bmatrix, MassI) 548 END IF 549 ! Compute the Metric Tensor 550 DO i = 1, nsize_int 551 DO j = 1, i 552 tmp = 0.0_dp 553 DO k = 1, nsize_xyz 554 tmp = tmp + Bmatrix(k, j)*MassI(k)*Bmatrix(k, i) 555 END DO 556 Mmatrix((i - 1)*nsize_int + j) = tmp 557 Mmatrix((j - 1)*nsize_int + i) = tmp 558 END DO 559 END DO 560 DEALLOCATE (MassI) 561 DEALLOCATE (Bmatrix) 562 END SUBROUTINE get_clv_force 563 564! ************************************************************************************************** 565!> \brief Complete the description of the COORDINATION colvar when 566!> defined using KINDS 567!> \param colvar ... 568!> \param particles ... 569!> \par History 570!> 1.2009 Fabio Sterpone : Added a part for population 571!> 10.2014 Moved out of colvar_types.F [Ole Schuett] 572!> \author Teodoro Laino - 07.2007 573! ************************************************************************************************** 574 SUBROUTINE post_process_colvar(colvar, particles) 575 TYPE(colvar_type), POINTER :: colvar 576 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 577 POINTER :: particles 578 579 CHARACTER(len=*), PARAMETER :: routineN = 'post_process_colvar', & 580 routineP = moduleN//':'//routineN 581 582 CHARACTER(len=3) :: name_kind 583 INTEGER :: i, ii, j, natoms, nkinds, nr_frame, stat 584 585 natoms = SIZE(particles) 586 IF (colvar%type_id == coord_colvar_id) THEN 587 IF (colvar%coord_param%use_kinds_from .OR. colvar%coord_param%use_kinds_to) THEN 588 ! Atoms from 589 IF (colvar%coord_param%use_kinds_from) THEN 590 colvar%coord_param%use_kinds_from = .FALSE. 591 nkinds = SIZE(colvar%coord_param%c_kinds_from) 592 DO i = 1, natoms 593 DO j = 1, nkinds 594 name_kind = TRIM(particles(i)%atomic_kind%name) 595 CALL uppercase(name_kind) 596 IF (TRIM(colvar%coord_param%c_kinds_from(j)) == name_kind) THEN 597 CALL reallocate(colvar%coord_param%i_at_from, 1, colvar%coord_param%n_atoms_from + 1) 598 colvar%coord_param%n_atoms_from = colvar%coord_param%n_atoms_from + 1 599 colvar%coord_param%i_at_from(colvar%coord_param%n_atoms_from) = i 600 END IF 601 END DO 602 END DO 603 stat = colvar%coord_param%n_atoms_from 604 CPASSERT(stat /= 0) 605 END IF 606 ! Atoms to 607 IF (colvar%coord_param%use_kinds_to) THEN 608 colvar%coord_param%use_kinds_to = .FALSE. 609 nkinds = SIZE(colvar%coord_param%c_kinds_to) 610 DO i = 1, natoms 611 DO j = 1, nkinds 612 name_kind = TRIM(particles(i)%atomic_kind%name) 613 CALL uppercase(name_kind) 614 IF (TRIM(colvar%coord_param%c_kinds_to(j)) == name_kind) THEN 615 CALL reallocate(colvar%coord_param%i_at_to, 1, colvar%coord_param%n_atoms_to + 1) 616 colvar%coord_param%n_atoms_to = colvar%coord_param%n_atoms_to + 1 617 colvar%coord_param%i_at_to(colvar%coord_param%n_atoms_to) = i 618 END IF 619 END DO 620 END DO 621 stat = colvar%coord_param%n_atoms_to 622 CPASSERT(stat /= 0) 623 END IF 624 ! Atoms to b 625 IF (colvar%coord_param%use_kinds_to_b) THEN 626 colvar%coord_param%use_kinds_to_b = .FALSE. 627 nkinds = SIZE(colvar%coord_param%c_kinds_to_b) 628 DO i = 1, natoms 629 DO j = 1, nkinds 630 name_kind = TRIM(particles(i)%atomic_kind%name) 631 CALL uppercase(name_kind) 632 IF (TRIM(colvar%coord_param%c_kinds_to_b(j)) == name_kind) THEN 633 CALL reallocate(colvar%coord_param%i_at_to_b, 1, colvar%coord_param%n_atoms_to_b + 1) 634 colvar%coord_param%n_atoms_to_b = colvar%coord_param%n_atoms_to_b + 1 635 colvar%coord_param%i_at_to_b(colvar%coord_param%n_atoms_to_b) = i 636 END IF 637 END DO 638 END DO 639 stat = colvar%coord_param%n_atoms_to_b 640 CPASSERT(stat /= 0) 641 END IF 642 643 ! Setup the colvar 644 CALL colvar_setup(colvar) 645 END IF 646 END IF 647 648 IF (colvar%type_id == mindist_colvar_id) THEN 649 IF (colvar%mindist_param%use_kinds_from .OR. colvar%mindist_param%use_kinds_to) THEN 650 ! Atoms from 651 IF (colvar%mindist_param%use_kinds_from) THEN 652 colvar%mindist_param%use_kinds_from = .FALSE. 653 nkinds = SIZE(colvar%mindist_param%k_coord_from) 654 DO i = 1, natoms 655 DO j = 1, nkinds 656 name_kind = TRIM(particles(i)%atomic_kind%name) 657 CALL uppercase(name_kind) 658 IF (TRIM(colvar%mindist_param%k_coord_from(j)) == name_kind) THEN 659 CALL reallocate(colvar%mindist_param%i_coord_from, 1, colvar%mindist_param%n_coord_from + 1) 660 colvar%mindist_param%n_coord_from = colvar%mindist_param%n_coord_from + 1 661 colvar%mindist_param%i_coord_from(colvar%mindist_param%n_coord_from) = i 662 END IF 663 END DO 664 END DO 665 stat = colvar%mindist_param%n_coord_from 666 CPASSERT(stat /= 0) 667 END IF 668 ! Atoms to 669 IF (colvar%mindist_param%use_kinds_to) THEN 670 colvar%mindist_param%use_kinds_to = .FALSE. 671 nkinds = SIZE(colvar%mindist_param%k_coord_to) 672 DO i = 1, natoms 673 DO j = 1, nkinds 674 name_kind = TRIM(particles(i)%atomic_kind%name) 675 CALL uppercase(name_kind) 676 IF (TRIM(colvar%mindist_param%k_coord_to(j)) == name_kind) THEN 677 CALL reallocate(colvar%mindist_param%i_coord_to, 1, colvar%mindist_param%n_coord_to + 1) 678 colvar%mindist_param%n_coord_to = colvar%mindist_param%n_coord_to + 1 679 colvar%mindist_param%i_coord_to(colvar%mindist_param%n_coord_to) = i 680 END IF 681 END DO 682 END DO 683 stat = colvar%mindist_param%n_coord_to 684 CPASSERT(stat /= 0) 685 END IF 686 ! Setup the colvar 687 CALL colvar_setup(colvar) 688 END IF 689 END IF 690 691 IF (colvar%type_id == population_colvar_id) THEN 692 693 IF (colvar%population_param%use_kinds_from .OR. colvar%population_param%use_kinds_to) THEN 694 ! Atoms from 695 IF (colvar%population_param%use_kinds_from) THEN 696 colvar%population_param%use_kinds_from = .FALSE. 697 nkinds = SIZE(colvar%population_param%c_kinds_from) 698 DO i = 1, natoms 699 DO j = 1, nkinds 700 name_kind = TRIM(particles(i)%atomic_kind%name) 701 CALL uppercase(name_kind) 702 IF (TRIM(colvar%population_param%c_kinds_from(j)) == name_kind) THEN 703 CALL reallocate(colvar%population_param%i_at_from, 1, colvar%population_param%n_atoms_from + 1) 704 colvar%population_param%n_atoms_from = colvar%population_param%n_atoms_from + 1 705 colvar%population_param%i_at_from(colvar%population_param%n_atoms_from) = i 706 END IF 707 END DO 708 END DO 709 stat = colvar%population_param%n_atoms_from 710 CPASSERT(stat /= 0) 711 END IF 712 ! Atoms to 713 IF (colvar%population_param%use_kinds_to) THEN 714 colvar%population_param%use_kinds_to = .FALSE. 715 nkinds = SIZE(colvar%population_param%c_kinds_to) 716 DO i = 1, natoms 717 DO j = 1, nkinds 718 name_kind = TRIM(particles(i)%atomic_kind%name) 719 CALL uppercase(name_kind) 720 IF (TRIM(colvar%population_param%c_kinds_to(j)) == name_kind) THEN 721 CALL reallocate(colvar%population_param%i_at_to, 1, colvar%population_param%n_atoms_to + 1) 722 colvar%population_param%n_atoms_to = colvar%population_param%n_atoms_to + 1 723 colvar%population_param%i_at_to(colvar%population_param%n_atoms_to) = i 724 END IF 725 END DO 726 END DO 727 stat = colvar%population_param%n_atoms_to 728 CPASSERT(stat /= 0) 729 END IF 730 ! Setup the colvar 731 CALL colvar_setup(colvar) 732 END IF 733 734 END IF 735 736 IF (colvar%type_id == gyration_colvar_id) THEN 737 738 IF (colvar%gyration_param%use_kinds) THEN 739 ! Atoms from 740 IF (colvar%gyration_param%use_kinds) THEN 741 colvar%gyration_param%use_kinds = .FALSE. 742 nkinds = SIZE(colvar%gyration_param%c_kinds) 743 DO i = 1, natoms 744 DO j = 1, nkinds 745 name_kind = TRIM(particles(i)%atomic_kind%name) 746 CALL uppercase(name_kind) 747 IF (TRIM(colvar%gyration_param%c_kinds(j)) == name_kind) THEN 748 CALL reallocate(colvar%gyration_param%i_at, 1, colvar%gyration_param%n_atoms + 1) 749 colvar%gyration_param%n_atoms = colvar%gyration_param%n_atoms + 1 750 colvar%gyration_param%i_at(colvar%gyration_param%n_atoms) = i 751 END IF 752 END DO 753 END DO 754 stat = colvar%gyration_param%n_atoms 755 CPASSERT(stat /= 0) 756 END IF 757 ! Setup the colvar 758 CALL colvar_setup(colvar) 759 END IF 760 END IF 761 762 IF (colvar%type_id == rmsd_colvar_id) THEN 763 IF (colvar%rmsd_param%subset == rmsd_all .OR. colvar%rmsd_param%subset == rmsd_list) THEN 764 ! weights are masses 765 DO i = 1, SIZE(colvar%rmsd_param%i_rmsd) 766 ii = colvar%rmsd_param%i_rmsd(i) 767 colvar%rmsd_param%weights(ii) = particles(ii)%atomic_kind%mass 768 END DO 769 END IF 770 771 IF (colvar%rmsd_param%align_frames) THEN 772 nr_frame = SIZE(colvar%rmsd_param%r_ref, 2) 773 DO i = 2, nr_frame 774 CALL rmsd3(particles, colvar%rmsd_param%r_ref(:, i), colvar%rmsd_param%r_ref(:, 1), -1, & 775 rotate=.TRUE.) 776 END DO 777 END IF 778 779 END IF 780 781 IF (colvar%type_id == distance_from_path_colvar_id .OR. colvar%type_id == reaction_path_colvar_id) THEN 782 IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd) THEN 783 IF (colvar%reaction_path_param%align_frames) THEN 784 nr_frame = colvar%reaction_path_param%nr_frames 785 DO i = 2, nr_frame 786 CALL rmsd3(particles, colvar%reaction_path_param%r_ref(:, i), colvar%reaction_path_param%r_ref(:, 1), -1, & 787 rotate=.TRUE.) 788 END DO 789 END IF 790 END IF 791 END IF 792 793 END SUBROUTINE post_process_colvar 794 795END MODULE colvar_utils 796