(git:ccc2433)
colvar_utils.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief evaluations of colvar for internal coordinates schemes
10 !> \par History
11 !> 05-2007 created [tlaino]
12 !> \author Teodoro Laino - Zurich University (2007) [tlaino]
13 ! **************************************************************************************************
14 MODULE colvar_utils
15  USE cell_types, ONLY: cell_type
16  USE colvar_methods, ONLY: colvar_eval_mol_f
17  USE colvar_types, ONLY: &
18  colvar_counters, colvar_setup, colvar_type, coord_colvar_id, distance_from_path_colvar_id, &
19  gyration_colvar_id, mindist_colvar_id, population_colvar_id, reaction_path_colvar_id, &
20  rmsd_colvar_id
21  USE cp_subsys_types, ONLY: cp_subsys_get,&
22  cp_subsys_type
23  USE distribution_1d_types, ONLY: distribution_1d_type
24  USE force_env_types, ONLY: force_env_get,&
25  force_env_type
26  USE input_constants, ONLY: rmsd_all,&
27  rmsd_list
28  USE kinds, ONLY: dp
29  USE mathlib, ONLY: invert_matrix
30  USE memory_utilities, ONLY: reallocate
31  USE molecule_kind_list_types, ONLY: molecule_kind_list_type
32  USE molecule_kind_types, ONLY: colvar_constraint_type,&
33  fixd_constraint_type,&
34  get_molecule_kind,&
35  molecule_kind_type
36  USE molecule_list_types, ONLY: molecule_list_type
37  USE molecule_types, ONLY: get_molecule,&
38  global_constraint_type,&
39  local_colvar_constraint_type,&
40  molecule_type
41  USE particle_list_types, ONLY: particle_list_type
42  USE particle_types, ONLY: particle_type
43  USE rmsd, ONLY: rmsd3
44  USE string_utilities, ONLY: uppercase
45  USE util, ONLY: sort
46 #include "./base/base_uses.f90"
47 
48  IMPLICIT NONE
49 
50  PRIVATE
51  PUBLIC :: number_of_colvar, &
52  eval_colvar, &
53  set_colvars_target, &
54  get_clv_force, &
55  post_process_colvar
56 
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'colvar_utils'
58 
59 CONTAINS
60 
61 ! **************************************************************************************************
62 !> \brief Gives back the number of colvar defined for a force_eval
63 !> \param force_env ...
64 !> \param only_intra_colvar ...
65 !> \param unique ...
66 !> \return ...
67 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
68 ! **************************************************************************************************
69  FUNCTION number_of_colvar(force_env, only_intra_colvar, unique) RESULT(ntot)
70  TYPE(force_env_type), POINTER :: force_env
71  LOGICAL, INTENT(IN), OPTIONAL :: only_intra_colvar, unique
72  INTEGER :: ntot
73 
74  CHARACTER(LEN=*), PARAMETER :: routinen = 'number_of_colvar'
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 
138  INTEGER :: handle, i, ikind, ind, nkind
139  TYPE(cell_type), POINTER :: cell
140  TYPE(colvar_constraint_type), DIMENSION(:), &
141  POINTER :: colv_list
142  TYPE(colvar_counters) :: ncolv
143  TYPE(cp_subsys_type), POINTER :: subsys
144  TYPE(global_constraint_type), POINTER :: gci
145  TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
146  TYPE(molecule_kind_type), POINTER :: molecule_kind
147 
148  NULLIFY (cell, subsys, molecule_kinds, molecule_kind, gci, colv_list)
149  CALL timeset(routinen, handle)
150  CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
151  CALL cp_subsys_get(subsys=subsys, gci=gci, molecule_kinds=molecule_kinds)
152 
153  nkind = molecule_kinds%n_els
154  ! Set Target for Intramolecular Colvars
155  mol: DO ikind = 1, nkind
156  molecule_kind => molecule_kinds%els(ikind)
157  CALL get_molecule_kind(molecule_kind, &
158  colv_list=colv_list, &
159  ncolv=ncolv)
160  IF (ncolv%ntot /= 0) THEN
161  DO i = 1, SIZE(colv_list)
162  ind = colv_list(i)%inp_seq_num
163  colv_list(i)%expected_value = targets(ind)
164  END DO
165  END IF
166  END DO mol
167  ! Set Target for Intermolecular Colvars
168  IF (ASSOCIATED(gci)) THEN
169  IF (gci%ncolv%ntot /= 0) THEN
170  colv_list => gci%colv_list
171  DO i = 1, SIZE(colv_list)
172  ind = colv_list(i)%inp_seq_num
173  colv_list(i)%expected_value = targets(ind)
174  END DO
175  END IF
176  END IF
177  CALL timestop(handle)
178 
179  END SUBROUTINE set_colvars_target
180 
181 ! **************************************************************************************************
182 !> \brief Computes the values of colvars and the Wilson matrix B and its invers A
183 !> \param force_env ...
184 !> \param coords ...
185 !> \param cvalues ...
186 !> \param Bmatrix ...
187 !> \param MassI ...
188 !> \param Amatrix ...
189 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
190 ! **************************************************************************************************
191  SUBROUTINE eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix)
192 
193  TYPE(force_env_type), POINTER :: force_env
194  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
195  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cvalues
196  REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
197  REAL(kind=dp), DIMENSION(:), OPTIONAL, POINTER :: massi
198  REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: amatrix
199 
200  CHARACTER(LEN=*), PARAMETER :: routinen = 'eval_colvar'
201 
202  INTEGER :: handle, i, ikind, imol, n_tot, natom, &
203  nkind, nmol_per_kind, offset
204  INTEGER, DIMENSION(:), POINTER :: map, wrk
205  LOGICAL :: check
206  REAL(kind=dp) :: inv_error
207  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: bwrk, gmatrix, gmatrix_i
208  REAL(kind=dp), DIMENSION(:), POINTER :: rwrk
209  TYPE(cell_type), POINTER :: cell
210  TYPE(colvar_counters) :: ncolv
211  TYPE(cp_subsys_type), POINTER :: subsys
212  TYPE(distribution_1d_type), POINTER :: local_molecules
213  TYPE(global_constraint_type), POINTER :: gci
214  TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
215  TYPE(molecule_kind_type), POINTER :: molecule_kind
216  TYPE(molecule_list_type), POINTER :: molecules
217  TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
218  TYPE(particle_list_type), POINTER :: particles
219  TYPE(particle_type), POINTER :: particle_set(:)
220 
221  NULLIFY (cell, subsys, local_molecules, molecule_kinds, &
222  molecules, molecule_kind, molecule, &
223  molecule_set, particles, particle_set, gci)
224  IF (PRESENT(bmatrix)) THEN
225  check = ASSOCIATED(bmatrix)
226  cpassert(check)
227  bmatrix = 0.0_dp
228  END IF
229  CALL timeset(routinen, handle)
230  ALLOCATE (map(SIZE(cvalues)))
231  map = huge(0) ! init all, since used in a sort, but not all set in parallel.
232  CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
233  n_tot = 0
234  cvalues = 0.0_dp
235  CALL cp_subsys_get(subsys=subsys, &
236  particles=particles, &
237  molecules=molecules, &
238  local_molecules=local_molecules, &
239  gci=gci, &
240  molecule_kinds=molecule_kinds)
241 
242  nkind = molecule_kinds%n_els
243  particle_set => particles%els
244  molecule_set => molecules%els
245  ! Intramolecular Colvars
246  IF (number_of_colvar(force_env, only_intra_colvar=.true.) /= 0) THEN
247  mol: DO ikind = 1, nkind
248  nmol_per_kind = local_molecules%n_el(ikind)
249  DO imol = 1, nmol_per_kind
250  i = local_molecules%list(ikind)%array(imol)
251  molecule => molecule_set(i)
252  molecule_kind => molecule%molecule_kind
253 
254  CALL get_molecule_kind(molecule_kind, &
255  ncolv=ncolv)
256  offset = get_colvar_offset(i, molecule_set)
257  ! Collective variables
258  IF (ncolv%ntot /= 0) THEN
259  CALL eval_colv_int(molecule, particle_set, coords, cell, cvalues, &
260  bmatrix, offset, n_tot, map)
261  END IF
262  END DO
263  END DO mol
264  CALL force_env%para_env%sum(n_tot)
265  CALL force_env%para_env%sum(cvalues)
266  IF (PRESENT(bmatrix)) CALL force_env%para_env%sum(bmatrix)
267  END IF
268  offset = n_tot
269  ! Intermolecular Colvars
270  IF (ASSOCIATED(gci)) THEN
271  IF (gci%ncolv%ntot /= 0) THEN
272  CALL eval_colv_ext(gci, particle_set, coords, cell, cvalues, &
273  bmatrix, offset, n_tot, map)
274  END IF
275  END IF
276  cpassert(n_tot == SIZE(cvalues))
277  ! Sort values of Collective Variables according the order of the input
278  ! sections
279  ALLOCATE (wrk(SIZE(cvalues)))
280  ALLOCATE (rwrk(SIZE(cvalues)))
281  CALL sort(map, SIZE(map), wrk)
282  rwrk = cvalues
283  DO i = 1, SIZE(wrk)
284  cvalues(i) = rwrk(wrk(i))
285  END DO
286  ! check and sort on Bmatrix
287  IF (PRESENT(bmatrix)) THEN
288  check = n_tot == SIZE(bmatrix, 2)
289  cpassert(check)
290  ALLOCATE (bwrk(SIZE(bmatrix, 1), SIZE(bmatrix, 2)))
291  bwrk(:, :) = bmatrix
292  DO i = 1, SIZE(wrk)
293  bmatrix(:, i) = bwrk(:, wrk(i))
294  END DO
295  DEALLOCATE (bwrk)
296  END IF
297  DEALLOCATE (rwrk)
298  DEALLOCATE (wrk)
299  DEALLOCATE (map)
300  ! Construction of the Amatrix
301  IF (PRESENT(bmatrix) .AND. PRESENT(amatrix)) THEN
302  cpassert(ASSOCIATED(amatrix))
303  check = SIZE(bmatrix, 1) == SIZE(amatrix, 2)
304  cpassert(check)
305  check = SIZE(bmatrix, 2) == SIZE(amatrix, 1)
306  cpassert(check)
307  ALLOCATE (gmatrix(n_tot, n_tot))
308  ALLOCATE (gmatrix_i(n_tot, n_tot))
309  gmatrix(:, :) = matmul(transpose(bmatrix), bmatrix)
310  CALL invert_matrix(gmatrix, gmatrix_i, inv_error)
311  IF (abs(inv_error) > 1.0e-8_dp) &
312  cpwarn("Error in inverting the Gmatrix larger than 1.0E-8!")
313  amatrix = matmul(gmatrix_i, transpose(bmatrix))
314  DEALLOCATE (gmatrix_i)
315  DEALLOCATE (gmatrix)
316  END IF
317  IF (PRESENT(massi)) THEN
318  natom = SIZE(particle_set)
319  cpassert(ASSOCIATED(massi))
320  cpassert(SIZE(massi) == natom*3)
321  DO i = 1, natom
322  massi((i - 1)*3 + 1) = 1.0_dp/particle_set(i)%atomic_kind%mass
323  massi((i - 1)*3 + 2) = 1.0_dp/particle_set(i)%atomic_kind%mass
324  massi((i - 1)*3 + 3) = 1.0_dp/particle_set(i)%atomic_kind%mass
325  END DO
326  END IF
327  CALL timestop(handle)
328 
329  END SUBROUTINE eval_colvar
330 
331 ! **************************************************************************************************
332 !> \brief Computes the offset of the colvar for the specific molecule
333 !> \param i ...
334 !> \param molecule_set ...
335 !> \return ...
336 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
337 ! **************************************************************************************************
338  FUNCTION get_colvar_offset(i, molecule_set) RESULT(offset)
339  INTEGER, INTENT(IN) :: i
340  TYPE(molecule_type), POINTER :: molecule_set(:)
341  INTEGER :: offset
342 
343  INTEGER :: j
344  TYPE(colvar_counters) :: ncolv
345  TYPE(molecule_kind_type), POINTER :: molecule_kind
346  TYPE(molecule_type), POINTER :: molecule
347 
348  offset = 0
349  DO j = 1, i - 1
350  molecule => molecule_set(j)
351  molecule_kind => molecule%molecule_kind
352  CALL get_molecule_kind(molecule_kind, &
353  ncolv=ncolv)
354  offset = offset + ncolv%ntot
355  END DO
356 
357  END FUNCTION get_colvar_offset
358 
359 ! **************************************************************************************************
360 !> \brief Computes Intramolecular colvar
361 !> \param molecule ...
362 !> \param particle_set ...
363 !> \param coords ...
364 !> \param cell ...
365 !> \param cvalues ...
366 !> \param Bmatrix ...
367 !> \param offset ...
368 !> \param n_tot ...
369 !> \param map ...
370 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
371 ! **************************************************************************************************
372  SUBROUTINE eval_colv_int(molecule, particle_set, coords, cell, cvalues, &
373  Bmatrix, offset, n_tot, map)
374 
375  TYPE(molecule_type), POINTER :: molecule
376  TYPE(particle_type), POINTER :: particle_set(:)
377  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
378  TYPE(cell_type), POINTER :: cell
379  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cvalues
380  REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
381  INTEGER, INTENT(IN) :: offset
382  INTEGER, INTENT(INOUT) :: n_tot
383  INTEGER, DIMENSION(:), POINTER :: map
384 
385  TYPE(colvar_constraint_type), POINTER :: colv_list(:)
386  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
387  TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
388  TYPE(molecule_kind_type), POINTER :: molecule_kind
389 
390  NULLIFY (fixd_list)
391 
392  molecule_kind => molecule%molecule_kind
393  CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
394  CALL get_molecule(molecule, lcolv=lcolv)
395  CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
396  coords, cell, cvalues, bmatrix, offset, n_tot, map)
397 
398  END SUBROUTINE eval_colv_int
399 
400 ! **************************************************************************************************
401 !> \brief Computes Intermolecular colvar
402 !> \param gci ...
403 !> \param particle_set ...
404 !> \param coords ...
405 !> \param cell ...
406 !> \param cvalues ...
407 !> \param Bmatrix ...
408 !> \param offset ...
409 !> \param n_tot ...
410 !> \param map ...
411 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
412 ! **************************************************************************************************
413  SUBROUTINE eval_colv_ext(gci, particle_set, coords, cell, cvalues, &
414  Bmatrix, offset, n_tot, map)
415  TYPE(global_constraint_type), POINTER :: gci
416  TYPE(particle_type), POINTER :: particle_set(:)
417  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
418  TYPE(cell_type), POINTER :: cell
419  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cvalues
420  REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
421  INTEGER, INTENT(IN) :: offset
422  INTEGER, INTENT(INOUT) :: n_tot
423  INTEGER, DIMENSION(:), POINTER :: map
424 
425  TYPE(colvar_constraint_type), POINTER :: colv_list(:)
426  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
427  TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
428 
429  colv_list => gci%colv_list
430  fixd_list => gci%fixd_list
431  lcolv => gci%lcolv
432  CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
433  coords, cell, cvalues, bmatrix, offset, n_tot, map)
434 
435  END SUBROUTINE eval_colv_ext
436 
437 ! **************************************************************************************************
438 !> \brief Real evaluation of colvar and of the Wilson-Eliashevich Matrix
439 !> B_ik : i: internal coordinates
440 !> k: cartesian coordinates
441 !> \param colv_list ...
442 !> \param fixd_list ...
443 !> \param lcolv ...
444 !> \param particle_set ...
445 !> \param coords ...
446 !> \param cell ...
447 !> \param cvalues ...
448 !> \param Bmatrix ...
449 !> \param offset ...
450 !> \param n_tot ...
451 !> \param map ...
452 !> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
453 ! **************************************************************************************************
454  SUBROUTINE eval_colv_low(colv_list, fixd_list, lcolv, particle_set, coords, &
455  cell, cvalues, Bmatrix, offset, n_tot, map)
456 
457  TYPE(colvar_constraint_type), POINTER :: colv_list(:)
458  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
459  TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
460  TYPE(particle_type), POINTER :: particle_set(:)
461  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
462  TYPE(cell_type), POINTER :: cell
463  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cvalues
464  REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
465  INTEGER, INTENT(IN) :: offset
466  INTEGER, INTENT(INOUT) :: n_tot
467  INTEGER, DIMENSION(:), POINTER :: map
468 
469  INTEGER :: iatm, iconst, ind, ival
470 
471  ival = offset
472  DO iconst = 1, SIZE(colv_list)
473  n_tot = n_tot + 1
474  ival = ival + 1
475  ! Update colvar
476  IF (PRESENT(coords)) THEN
477  CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
478  pos=reshape(coords, (/3, SIZE(particle_set)/)), fixd_list=fixd_list)
479  ELSE
480  CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
481  fixd_list=fixd_list)
482  END IF
483  cvalues(ival) = lcolv(iconst)%colvar%ss
484  map(ival) = colv_list(iconst)%inp_seq_num
485  ! Build the Wilson-Eliashevich Matrix
486  IF (PRESENT(bmatrix)) THEN
487  DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom)
488  ind = (lcolv(iconst)%colvar%i_atom(iatm) - 1)*3
489  bmatrix(ind + 1, ival) = lcolv(iconst)%colvar%dsdr(1, iatm)
490  bmatrix(ind + 2, ival) = lcolv(iconst)%colvar%dsdr(2, iatm)
491  bmatrix(ind + 3, ival) = lcolv(iconst)%colvar%dsdr(3, iatm)
492  END DO
493  END IF
494  END DO
495 
496  END SUBROUTINE eval_colv_low
497 
498 ! **************************************************************************************************
499 !> \brief Computes the forces in the frame of collective variables, and additional
500 !> also the local metric tensor
501 !> \param force_env ...
502 !> \param forces ...
503 !> \param coords ...
504 !> \param nsize_xyz ...
505 !> \param nsize_int ...
506 !> \param cvalues ...
507 !> \param Mmatrix ...
508 !> \author Teodoro Laino 05.2007
509 ! **************************************************************************************************
510  SUBROUTINE get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, &
511  Mmatrix)
512  TYPE(force_env_type), POINTER :: force_env
513  REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
514  OPTIONAL :: forces, coords
515  INTEGER, INTENT(IN) :: nsize_xyz, nsize_int
516  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cvalues, mmatrix
517 
518  INTEGER :: i, j, k
519  REAL(kind=dp) :: tmp
520  REAL(kind=dp), DIMENSION(:), POINTER :: massi, wrk
521  REAL(kind=dp), DIMENSION(:, :), POINTER :: amatrix, bmatrix
522 
523  ALLOCATE (bmatrix(nsize_xyz, nsize_int))
524  ALLOCATE (massi(nsize_xyz))
525  ! Transform gradients if requested
526  IF (PRESENT(forces)) THEN
527  ALLOCATE (wrk(nsize_int))
528  ALLOCATE (amatrix(nsize_int, nsize_xyz))
529  ! Compute the transformation matrices and the inverse mass diagonal Matrix
530  CALL eval_colvar(force_env, coords, cvalues, bmatrix, massi, amatrix)
531  wrk = matmul(amatrix, forces)
532  forces = 0.0_dp
533  forces(1:nsize_int) = wrk
534  DEALLOCATE (amatrix)
535  DEALLOCATE (wrk)
536  ELSE
537  ! Compute the transformation matrices and the inverse mass diagonal Matrix
538  CALL eval_colvar(force_env, coords, cvalues, bmatrix, massi)
539  END IF
540  ! Compute the Metric Tensor
541  DO i = 1, nsize_int
542  DO j = 1, i
543  tmp = 0.0_dp
544  DO k = 1, nsize_xyz
545  tmp = tmp + bmatrix(k, j)*massi(k)*bmatrix(k, i)
546  END DO
547  mmatrix((i - 1)*nsize_int + j) = tmp
548  mmatrix((j - 1)*nsize_int + i) = tmp
549  END DO
550  END DO
551  DEALLOCATE (massi)
552  DEALLOCATE (bmatrix)
553  END SUBROUTINE get_clv_force
554 
555 ! **************************************************************************************************
556 !> \brief Complete the description of the COORDINATION colvar when
557 !> defined using KINDS
558 !> \param colvar ...
559 !> \param particles ...
560 !> \par History
561 !> 1.2009 Fabio Sterpone : Added a part for population
562 !> 10.2014 Moved out of colvar_types.F [Ole Schuett]
563 !> \author Teodoro Laino - 07.2007
564 ! **************************************************************************************************
565  SUBROUTINE post_process_colvar(colvar, particles)
566  TYPE(colvar_type), POINTER :: colvar
567  TYPE(particle_type), DIMENSION(:), OPTIONAL, &
568  POINTER :: particles
569 
570  CHARACTER(len=3) :: name_kind
571  INTEGER :: i, ii, j, natoms, nkinds, nr_frame, stat
572 
573  natoms = SIZE(particles)
574  IF (colvar%type_id == coord_colvar_id) THEN
575  IF (colvar%coord_param%use_kinds_from .OR. colvar%coord_param%use_kinds_to) THEN
576  ! Atoms from
577  IF (colvar%coord_param%use_kinds_from) THEN
578  colvar%coord_param%use_kinds_from = .false.
579  nkinds = SIZE(colvar%coord_param%c_kinds_from)
580  DO i = 1, natoms
581  DO j = 1, nkinds
582  name_kind = trim(particles(i)%atomic_kind%name)
583  CALL uppercase(name_kind)
584  IF (trim(colvar%coord_param%c_kinds_from(j)) == name_kind) THEN
585  CALL reallocate(colvar%coord_param%i_at_from, 1, colvar%coord_param%n_atoms_from + 1)
586  colvar%coord_param%n_atoms_from = colvar%coord_param%n_atoms_from + 1
587  colvar%coord_param%i_at_from(colvar%coord_param%n_atoms_from) = i
588  END IF
589  END DO
590  END DO
591  stat = colvar%coord_param%n_atoms_from
592  cpassert(stat /= 0)
593  END IF
594  ! Atoms to
595  IF (colvar%coord_param%use_kinds_to) THEN
596  colvar%coord_param%use_kinds_to = .false.
597  nkinds = SIZE(colvar%coord_param%c_kinds_to)
598  DO i = 1, natoms
599  DO j = 1, nkinds
600  name_kind = trim(particles(i)%atomic_kind%name)
601  CALL uppercase(name_kind)
602  IF (trim(colvar%coord_param%c_kinds_to(j)) == name_kind) THEN
603  CALL reallocate(colvar%coord_param%i_at_to, 1, colvar%coord_param%n_atoms_to + 1)
604  colvar%coord_param%n_atoms_to = colvar%coord_param%n_atoms_to + 1
605  colvar%coord_param%i_at_to(colvar%coord_param%n_atoms_to) = i
606  END IF
607  END DO
608  END DO
609  stat = colvar%coord_param%n_atoms_to
610  cpassert(stat /= 0)
611  END IF
612  ! Atoms to b
613  IF (colvar%coord_param%use_kinds_to_b) THEN
614  colvar%coord_param%use_kinds_to_b = .false.
615  nkinds = SIZE(colvar%coord_param%c_kinds_to_b)
616  DO i = 1, natoms
617  DO j = 1, nkinds
618  name_kind = trim(particles(i)%atomic_kind%name)
619  CALL uppercase(name_kind)
620  IF (trim(colvar%coord_param%c_kinds_to_b(j)) == name_kind) THEN
621  CALL reallocate(colvar%coord_param%i_at_to_b, 1, colvar%coord_param%n_atoms_to_b + 1)
622  colvar%coord_param%n_atoms_to_b = colvar%coord_param%n_atoms_to_b + 1
623  colvar%coord_param%i_at_to_b(colvar%coord_param%n_atoms_to_b) = i
624  END IF
625  END DO
626  END DO
627  stat = colvar%coord_param%n_atoms_to_b
628  cpassert(stat /= 0)
629  END IF
630 
631  ! Setup the colvar
632  CALL colvar_setup(colvar)
633  END IF
634  END IF
635 
636  IF (colvar%type_id == mindist_colvar_id) THEN
637  IF (colvar%mindist_param%use_kinds_from .OR. colvar%mindist_param%use_kinds_to) THEN
638  ! Atoms from
639  IF (colvar%mindist_param%use_kinds_from) THEN
640  colvar%mindist_param%use_kinds_from = .false.
641  nkinds = SIZE(colvar%mindist_param%k_coord_from)
642  DO i = 1, natoms
643  DO j = 1, nkinds
644  name_kind = trim(particles(i)%atomic_kind%name)
645  CALL uppercase(name_kind)
646  IF (trim(colvar%mindist_param%k_coord_from(j)) == name_kind) THEN
647  CALL reallocate(colvar%mindist_param%i_coord_from, 1, colvar%mindist_param%n_coord_from + 1)
648  colvar%mindist_param%n_coord_from = colvar%mindist_param%n_coord_from + 1
649  colvar%mindist_param%i_coord_from(colvar%mindist_param%n_coord_from) = i
650  END IF
651  END DO
652  END DO
653  stat = colvar%mindist_param%n_coord_from
654  cpassert(stat /= 0)
655  END IF
656  ! Atoms to
657  IF (colvar%mindist_param%use_kinds_to) THEN
658  colvar%mindist_param%use_kinds_to = .false.
659  nkinds = SIZE(colvar%mindist_param%k_coord_to)
660  DO i = 1, natoms
661  DO j = 1, nkinds
662  name_kind = trim(particles(i)%atomic_kind%name)
663  CALL uppercase(name_kind)
664  IF (trim(colvar%mindist_param%k_coord_to(j)) == name_kind) THEN
665  CALL reallocate(colvar%mindist_param%i_coord_to, 1, colvar%mindist_param%n_coord_to + 1)
666  colvar%mindist_param%n_coord_to = colvar%mindist_param%n_coord_to + 1
667  colvar%mindist_param%i_coord_to(colvar%mindist_param%n_coord_to) = i
668  END IF
669  END DO
670  END DO
671  stat = colvar%mindist_param%n_coord_to
672  cpassert(stat /= 0)
673  END IF
674  ! Setup the colvar
675  CALL colvar_setup(colvar)
676  END IF
677  END IF
678 
679  IF (colvar%type_id == population_colvar_id) THEN
680 
681  IF (colvar%population_param%use_kinds_from .OR. colvar%population_param%use_kinds_to) THEN
682  ! Atoms from
683  IF (colvar%population_param%use_kinds_from) THEN
684  colvar%population_param%use_kinds_from = .false.
685  nkinds = SIZE(colvar%population_param%c_kinds_from)
686  DO i = 1, natoms
687  DO j = 1, nkinds
688  name_kind = trim(particles(i)%atomic_kind%name)
689  CALL uppercase(name_kind)
690  IF (trim(colvar%population_param%c_kinds_from(j)) == name_kind) THEN
691  CALL reallocate(colvar%population_param%i_at_from, 1, colvar%population_param%n_atoms_from + 1)
692  colvar%population_param%n_atoms_from = colvar%population_param%n_atoms_from + 1
693  colvar%population_param%i_at_from(colvar%population_param%n_atoms_from) = i
694  END IF
695  END DO
696  END DO
697  stat = colvar%population_param%n_atoms_from
698  cpassert(stat /= 0)
699  END IF
700  ! Atoms to
701  IF (colvar%population_param%use_kinds_to) THEN
702  colvar%population_param%use_kinds_to = .false.
703  nkinds = SIZE(colvar%population_param%c_kinds_to)
704  DO i = 1, natoms
705  DO j = 1, nkinds
706  name_kind = trim(particles(i)%atomic_kind%name)
707  CALL uppercase(name_kind)
708  IF (trim(colvar%population_param%c_kinds_to(j)) == name_kind) THEN
709  CALL reallocate(colvar%population_param%i_at_to, 1, colvar%population_param%n_atoms_to + 1)
710  colvar%population_param%n_atoms_to = colvar%population_param%n_atoms_to + 1
711  colvar%population_param%i_at_to(colvar%population_param%n_atoms_to) = i
712  END IF
713  END DO
714  END DO
715  stat = colvar%population_param%n_atoms_to
716  cpassert(stat /= 0)
717  END IF
718  ! Setup the colvar
719  CALL colvar_setup(colvar)
720  END IF
721 
722  END IF
723 
724  IF (colvar%type_id == gyration_colvar_id) THEN
725 
726  IF (colvar%gyration_param%use_kinds) THEN
727  ! Atoms from
728  IF (colvar%gyration_param%use_kinds) THEN
729  colvar%gyration_param%use_kinds = .false.
730  nkinds = SIZE(colvar%gyration_param%c_kinds)
731  DO i = 1, natoms
732  DO j = 1, nkinds
733  name_kind = trim(particles(i)%atomic_kind%name)
734  CALL uppercase(name_kind)
735  IF (trim(colvar%gyration_param%c_kinds(j)) == name_kind) THEN
736  CALL reallocate(colvar%gyration_param%i_at, 1, colvar%gyration_param%n_atoms + 1)
737  colvar%gyration_param%n_atoms = colvar%gyration_param%n_atoms + 1
738  colvar%gyration_param%i_at(colvar%gyration_param%n_atoms) = i
739  END IF
740  END DO
741  END DO
742  stat = colvar%gyration_param%n_atoms
743  cpassert(stat /= 0)
744  END IF
745  ! Setup the colvar
746  CALL colvar_setup(colvar)
747  END IF
748  END IF
749 
750  IF (colvar%type_id == rmsd_colvar_id) THEN
751  IF (colvar%rmsd_param%subset == rmsd_all .OR. colvar%rmsd_param%subset == rmsd_list) THEN
752  ! weights are masses
753  DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
754  ii = colvar%rmsd_param%i_rmsd(i)
755  colvar%rmsd_param%weights(ii) = particles(ii)%atomic_kind%mass
756  END DO
757  END IF
758 
759  IF (colvar%rmsd_param%align_frames) THEN
760  nr_frame = SIZE(colvar%rmsd_param%r_ref, 2)
761  DO i = 2, nr_frame
762  CALL rmsd3(particles, colvar%rmsd_param%r_ref(:, i), colvar%rmsd_param%r_ref(:, 1), -1, &
763  rotate=.true.)
764  END DO
765  END IF
766 
767  END IF
768 
769  IF (colvar%type_id == distance_from_path_colvar_id .OR. colvar%type_id == reaction_path_colvar_id) THEN
770  IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd) THEN
771  IF (colvar%reaction_path_param%align_frames) THEN
772  nr_frame = colvar%reaction_path_param%nr_frames
773  DO i = 2, nr_frame
774  CALL rmsd3(particles, colvar%reaction_path_param%r_ref(:, i), colvar%reaction_path_param%r_ref(:, 1), -1, &
775  rotate=.true.)
776  END DO
777  END IF
778  END IF
779  END IF
780 
781  END SUBROUTINE post_process_colvar
782 
783 END MODULE colvar_utils
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition: colvar_methods.F:16
colvar_methods::colvar_eval_mol_f
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Definition: colvar_methods.F:1631
colvar_types
Initialize the collective variables types.
Definition: colvar_types.F:15
colvar_types::population_colvar_id
integer, parameter, public population_colvar_id
Definition: colvar_types.F:56
colvar_types::distance_from_path_colvar_id
integer, parameter, public distance_from_path_colvar_id
Definition: colvar_types.F:56
colvar_types::rmsd_colvar_id
integer, parameter, public rmsd_colvar_id
Definition: colvar_types.F:56
colvar_types::mindist_colvar_id
integer, parameter, public mindist_colvar_id
Definition: colvar_types.F:56
colvar_types::gyration_colvar_id
integer, parameter, public gyration_colvar_id
Definition: colvar_types.F:56
colvar_types::coord_colvar_id
integer, parameter, public coord_colvar_id
Definition: colvar_types.F:56
colvar_types::colvar_setup
subroutine, public colvar_setup(colvar)
Finalize the setup of the collective variable.
Definition: colvar_types.F:502
colvar_types::reaction_path_colvar_id
integer, parameter, public reaction_path_colvar_id
Definition: colvar_types.F:56
colvar_utils
evaluations of colvar for internal coordinates schemes
Definition: colvar_utils.F:14
colvar_utils::eval_colvar
subroutine, public eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix)
Computes the values of colvars and the Wilson matrix B and its invers A.
Definition: colvar_utils.F:192
colvar_utils::get_clv_force
subroutine, public get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, Mmatrix)
Computes the forces in the frame of collective variables, and additional also the local metric tensor...
Definition: colvar_utils.F:512
colvar_utils::post_process_colvar
subroutine, public post_process_colvar(colvar, particles)
Complete the description of the COORDINATION colvar when defined using KINDS.
Definition: colvar_utils.F:566
colvar_utils::number_of_colvar
integer function, public number_of_colvar(force_env, only_intra_colvar, unique)
Gives back the number of colvar defined for a force_eval.
Definition: colvar_utils.F:70
colvar_utils::set_colvars_target
subroutine, public set_colvars_target(targets, force_env)
Set the value of target for constraints/restraints.
Definition: colvar_utils.F:133
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition: cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition: cp_subsys_types.F:345
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition: force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
Definition: force_env_types.F:329
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::rmsd_list
integer, parameter, public rmsd_list
Definition: input_constants.F:87
input_constants::rmsd_all
integer, parameter, public rmsd_all
Definition: input_constants.F:87
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
molecule_kind_list_types
represent a simple array based list of the given type
Definition: molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition: molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition: molecule_kind_types.F:510
molecule_list_types
represent a simple array based list of the given type
Definition: molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition: molecule_types.F:277
particle_list_types
represent a simple array based list of the given type
Definition: particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
rmsd
Defines functions to perform rmsd in 3D.
Definition: rmsd.F:12
rmsd::rmsd3
subroutine, public rmsd3(particle_set, r, r0, output_unit, weights, my_val, rotate, transl, rot, drmsd3)
Computes the RMSD in 3D. Provides also derivatives.
Definition: rmsd.F:53
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361
util
All kind of helpful little routines.
Definition: util.F:14