(git:ed6f26b)
Loading...
Searching...
No Matches
colvar_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
14MODULE 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
59CONTAINS
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) THEN
312 cpwarn("Error in inverting the Gmatrix larger than 1.0E-8!")
313 END IF
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 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
387 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
388 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
389 TYPE(molecule_kind_type), POINTER :: molecule_kind
390
391 NULLIFY (fixd_list)
392
393 molecule_kind => molecule%molecule_kind
394 CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
395 CALL get_molecule(molecule, lcolv=lcolv)
396 CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
397 coords, cell, cvalues, bmatrix, offset, n_tot, map)
398
399 END SUBROUTINE eval_colv_int
400
401! **************************************************************************************************
402!> \brief Computes Intermolecular colvar
403!> \param gci ...
404!> \param particle_set ...
405!> \param coords ...
406!> \param cell ...
407!> \param cvalues ...
408!> \param Bmatrix ...
409!> \param offset ...
410!> \param n_tot ...
411!> \param map ...
412!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
413! **************************************************************************************************
414 SUBROUTINE eval_colv_ext(gci, particle_set, coords, cell, cvalues, &
415 Bmatrix, offset, n_tot, map)
416 TYPE(global_constraint_type), POINTER :: gci
417 TYPE(particle_type), POINTER :: particle_set(:)
418 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
419 TYPE(cell_type), POINTER :: cell
420 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cvalues
421 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
422 INTEGER, INTENT(IN) :: offset
423 INTEGER, INTENT(INOUT) :: n_tot
424 INTEGER, DIMENSION(:), POINTER :: map
425
426 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
427 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
428 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
429
430 colv_list => gci%colv_list
431 fixd_list => gci%fixd_list
432 lcolv => gci%lcolv
433 CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
434 coords, cell, cvalues, bmatrix, offset, n_tot, map)
435
436 END SUBROUTINE eval_colv_ext
437
438! **************************************************************************************************
439!> \brief Real evaluation of colvar and of the Wilson-Eliashevich Matrix
440!> B_ik : i: internal coordinates
441!> k: cartesian coordinates
442!> \param colv_list ...
443!> \param fixd_list ...
444!> \param lcolv ...
445!> \param particle_set ...
446!> \param coords ...
447!> \param cell ...
448!> \param cvalues ...
449!> \param Bmatrix ...
450!> \param offset ...
451!> \param n_tot ...
452!> \param map ...
453!> \author Teodoro Laino 05.2007 [tlaino] - Zurich University
454! **************************************************************************************************
455 SUBROUTINE eval_colv_low(colv_list, fixd_list, lcolv, particle_set, coords, &
456 cell, cvalues, Bmatrix, offset, n_tot, map)
457
458 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
459 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
460 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
461 TYPE(particle_type), POINTER :: particle_set(:)
462 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: coords
463 TYPE(cell_type), POINTER :: cell
464 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: cvalues
465 REAL(kind=dp), DIMENSION(:, :), OPTIONAL, POINTER :: bmatrix
466 INTEGER, INTENT(IN) :: offset
467 INTEGER, INTENT(INOUT) :: n_tot
468 INTEGER, DIMENSION(:), POINTER :: map
469
470 INTEGER :: iatm, iconst, ind, ival
471
472 ival = offset
473 DO iconst = 1, SIZE(colv_list)
474 n_tot = n_tot + 1
475 ival = ival + 1
476 ! Update colvar
477 IF (PRESENT(coords)) THEN
478 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
479 pos=reshape(coords, (/3, SIZE(particle_set)/)), fixd_list=fixd_list)
480 ELSE
481 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
482 fixd_list=fixd_list)
483 END IF
484 cvalues(ival) = lcolv(iconst)%colvar%ss
485 map(ival) = colv_list(iconst)%inp_seq_num
486 ! Build the Wilson-Eliashevich Matrix
487 IF (PRESENT(bmatrix)) THEN
488 DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom)
489 ind = (lcolv(iconst)%colvar%i_atom(iatm) - 1)*3
490 bmatrix(ind + 1, ival) = lcolv(iconst)%colvar%dsdr(1, iatm)
491 bmatrix(ind + 2, ival) = lcolv(iconst)%colvar%dsdr(2, iatm)
492 bmatrix(ind + 3, ival) = lcolv(iconst)%colvar%dsdr(3, iatm)
493 END DO
494 END IF
495 END DO
496
497 END SUBROUTINE eval_colv_low
498
499! **************************************************************************************************
500!> \brief Computes the forces in the frame of collective variables, and additional
501!> also the local metric tensor
502!> \param force_env ...
503!> \param forces ...
504!> \param coords ...
505!> \param nsize_xyz ...
506!> \param nsize_int ...
507!> \param cvalues ...
508!> \param Mmatrix ...
509!> \author Teodoro Laino 05.2007
510! **************************************************************************************************
511 SUBROUTINE get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, &
512 Mmatrix)
513 TYPE(force_env_type), POINTER :: force_env
514 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
515 OPTIONAL :: forces, coords
516 INTEGER, INTENT(IN) :: nsize_xyz, nsize_int
517 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: cvalues, mmatrix
518
519 INTEGER :: i, j, k
520 REAL(kind=dp) :: tmp
521 REAL(kind=dp), DIMENSION(:), POINTER :: massi, wrk
522 REAL(kind=dp), DIMENSION(:, :), POINTER :: amatrix, bmatrix
523
524 ALLOCATE (bmatrix(nsize_xyz, nsize_int))
525 ALLOCATE (massi(nsize_xyz))
526 ! Transform gradients if requested
527 IF (PRESENT(forces)) THEN
528 ALLOCATE (wrk(nsize_int))
529 ALLOCATE (amatrix(nsize_int, nsize_xyz))
530 ! Compute the transformation matrices and the inverse mass diagonal Matrix
531 CALL eval_colvar(force_env, coords, cvalues, bmatrix, massi, amatrix)
532 wrk = matmul(amatrix, forces)
533 forces = 0.0_dp
534 forces(1:nsize_int) = wrk
535 DEALLOCATE (amatrix)
536 DEALLOCATE (wrk)
537 ELSE
538 ! Compute the transformation matrices and the inverse mass diagonal Matrix
539 CALL eval_colvar(force_env, coords, cvalues, bmatrix, massi)
540 END IF
541 ! Compute the Metric Tensor
542 DO i = 1, nsize_int
543 DO j = 1, i
544 tmp = 0.0_dp
545 DO k = 1, nsize_xyz
546 tmp = tmp + bmatrix(k, j)*massi(k)*bmatrix(k, i)
547 END DO
548 mmatrix((i - 1)*nsize_int + j) = tmp
549 mmatrix((j - 1)*nsize_int + i) = tmp
550 END DO
551 END DO
552 DEALLOCATE (massi)
553 DEALLOCATE (bmatrix)
554 END SUBROUTINE get_clv_force
555
556! **************************************************************************************************
557!> \brief Complete the description of the COORDINATION colvar when
558!> defined using KINDS
559!> \param colvar ...
560!> \param particles ...
561!> \par History
562!> 1.2009 Fabio Sterpone : Added a part for population
563!> 10.2014 Moved out of colvar_types.F [Ole Schuett]
564!> \author Teodoro Laino - 07.2007
565! **************************************************************************************************
566 SUBROUTINE post_process_colvar(colvar, particles)
567 TYPE(colvar_type), POINTER :: colvar
568 TYPE(particle_type), DIMENSION(:), OPTIONAL, &
569 POINTER :: particles
570
571 CHARACTER(len=3) :: name_kind
572 INTEGER :: i, ii, j, natoms, nkinds, nr_frame, stat
573 REAL(kind=dp) :: mtot
574
575 natoms = SIZE(particles)
576 IF (colvar%type_id == coord_colvar_id) THEN
577 IF (colvar%coord_param%use_kinds_from .OR. colvar%coord_param%use_kinds_to) THEN
578 ! Atoms from
579 IF (colvar%coord_param%use_kinds_from) THEN
580 colvar%coord_param%use_kinds_from = .false.
581 nkinds = SIZE(colvar%coord_param%c_kinds_from)
582 DO i = 1, natoms
583 DO j = 1, nkinds
584 name_kind = trim(particles(i)%atomic_kind%name)
585 CALL uppercase(name_kind)
586 IF (trim(colvar%coord_param%c_kinds_from(j)) == name_kind) THEN
587 CALL reallocate(colvar%coord_param%i_at_from, 1, colvar%coord_param%n_atoms_from + 1)
588 colvar%coord_param%n_atoms_from = colvar%coord_param%n_atoms_from + 1
589 colvar%coord_param%i_at_from(colvar%coord_param%n_atoms_from) = i
590 END IF
591 END DO
592 END DO
593 stat = colvar%coord_param%n_atoms_from
594 cpassert(stat /= 0)
595 END IF
596 ! Atoms to
597 IF (colvar%coord_param%use_kinds_to) THEN
598 colvar%coord_param%use_kinds_to = .false.
599 nkinds = SIZE(colvar%coord_param%c_kinds_to)
600 DO i = 1, natoms
601 DO j = 1, nkinds
602 name_kind = trim(particles(i)%atomic_kind%name)
603 CALL uppercase(name_kind)
604 IF (trim(colvar%coord_param%c_kinds_to(j)) == name_kind) THEN
605 CALL reallocate(colvar%coord_param%i_at_to, 1, colvar%coord_param%n_atoms_to + 1)
606 colvar%coord_param%n_atoms_to = colvar%coord_param%n_atoms_to + 1
607 colvar%coord_param%i_at_to(colvar%coord_param%n_atoms_to) = i
608 END IF
609 END DO
610 END DO
611 stat = colvar%coord_param%n_atoms_to
612 cpassert(stat /= 0)
613 END IF
614 ! Atoms to b
615 IF (colvar%coord_param%use_kinds_to_b) THEN
616 colvar%coord_param%use_kinds_to_b = .false.
617 nkinds = SIZE(colvar%coord_param%c_kinds_to_b)
618 DO i = 1, natoms
619 DO j = 1, nkinds
620 name_kind = trim(particles(i)%atomic_kind%name)
621 CALL uppercase(name_kind)
622 IF (trim(colvar%coord_param%c_kinds_to_b(j)) == name_kind) THEN
623 CALL reallocate(colvar%coord_param%i_at_to_b, 1, colvar%coord_param%n_atoms_to_b + 1)
624 colvar%coord_param%n_atoms_to_b = colvar%coord_param%n_atoms_to_b + 1
625 colvar%coord_param%i_at_to_b(colvar%coord_param%n_atoms_to_b) = i
626 END IF
627 END DO
628 END DO
629 stat = colvar%coord_param%n_atoms_to_b
630 cpassert(stat /= 0)
631 END IF
632
633 ! Setup the colvar
634 CALL colvar_setup(colvar)
635 END IF
636 END IF
637
638 IF (colvar%type_id == mindist_colvar_id) THEN
639 IF (colvar%mindist_param%use_kinds_from .OR. colvar%mindist_param%use_kinds_to) THEN
640 ! Atoms from
641 IF (colvar%mindist_param%use_kinds_from) THEN
642 colvar%mindist_param%use_kinds_from = .false.
643 nkinds = SIZE(colvar%mindist_param%k_coord_from)
644 DO i = 1, natoms
645 DO j = 1, nkinds
646 name_kind = trim(particles(i)%atomic_kind%name)
647 CALL uppercase(name_kind)
648 IF (trim(colvar%mindist_param%k_coord_from(j)) == name_kind) THEN
649 CALL reallocate(colvar%mindist_param%i_coord_from, 1, colvar%mindist_param%n_coord_from + 1)
650 colvar%mindist_param%n_coord_from = colvar%mindist_param%n_coord_from + 1
651 colvar%mindist_param%i_coord_from(colvar%mindist_param%n_coord_from) = i
652 END IF
653 END DO
654 END DO
655 stat = colvar%mindist_param%n_coord_from
656 cpassert(stat /= 0)
657 END IF
658 ! Atoms to
659 IF (colvar%mindist_param%use_kinds_to) THEN
660 colvar%mindist_param%use_kinds_to = .false.
661 nkinds = SIZE(colvar%mindist_param%k_coord_to)
662 DO i = 1, natoms
663 DO j = 1, nkinds
664 name_kind = trim(particles(i)%atomic_kind%name)
665 CALL uppercase(name_kind)
666 IF (trim(colvar%mindist_param%k_coord_to(j)) == name_kind) THEN
667 CALL reallocate(colvar%mindist_param%i_coord_to, 1, colvar%mindist_param%n_coord_to + 1)
668 colvar%mindist_param%n_coord_to = colvar%mindist_param%n_coord_to + 1
669 colvar%mindist_param%i_coord_to(colvar%mindist_param%n_coord_to) = i
670 END IF
671 END DO
672 END DO
673 stat = colvar%mindist_param%n_coord_to
674 cpassert(stat /= 0)
675 END IF
676 ! Setup the colvar
677 CALL colvar_setup(colvar)
678 END IF
679 END IF
680
681 IF (colvar%type_id == population_colvar_id) THEN
682
683 IF (colvar%population_param%use_kinds_from .OR. colvar%population_param%use_kinds_to) THEN
684 ! Atoms from
685 IF (colvar%population_param%use_kinds_from) THEN
686 colvar%population_param%use_kinds_from = .false.
687 nkinds = SIZE(colvar%population_param%c_kinds_from)
688 DO i = 1, natoms
689 DO j = 1, nkinds
690 name_kind = trim(particles(i)%atomic_kind%name)
691 CALL uppercase(name_kind)
692 IF (trim(colvar%population_param%c_kinds_from(j)) == name_kind) THEN
693 CALL reallocate(colvar%population_param%i_at_from, 1, colvar%population_param%n_atoms_from + 1)
694 colvar%population_param%n_atoms_from = colvar%population_param%n_atoms_from + 1
695 colvar%population_param%i_at_from(colvar%population_param%n_atoms_from) = i
696 END IF
697 END DO
698 END DO
699 stat = colvar%population_param%n_atoms_from
700 cpassert(stat /= 0)
701 END IF
702 ! Atoms to
703 IF (colvar%population_param%use_kinds_to) THEN
704 colvar%population_param%use_kinds_to = .false.
705 nkinds = SIZE(colvar%population_param%c_kinds_to)
706 DO i = 1, natoms
707 DO j = 1, nkinds
708 name_kind = trim(particles(i)%atomic_kind%name)
709 CALL uppercase(name_kind)
710 IF (trim(colvar%population_param%c_kinds_to(j)) == name_kind) THEN
711 CALL reallocate(colvar%population_param%i_at_to, 1, colvar%population_param%n_atoms_to + 1)
712 colvar%population_param%n_atoms_to = colvar%population_param%n_atoms_to + 1
713 colvar%population_param%i_at_to(colvar%population_param%n_atoms_to) = i
714 END IF
715 END DO
716 END DO
717 stat = colvar%population_param%n_atoms_to
718 cpassert(stat /= 0)
719 END IF
720 ! Setup the colvar
721 CALL colvar_setup(colvar)
722 END IF
723
724 END IF
725
726 IF (colvar%type_id == gyration_colvar_id) THEN
727
728 IF (colvar%gyration_param%use_kinds) THEN
729 ! Atoms from
730 IF (colvar%gyration_param%use_kinds) THEN
731 colvar%gyration_param%use_kinds = .false.
732 nkinds = SIZE(colvar%gyration_param%c_kinds)
733 DO i = 1, natoms
734 DO j = 1, nkinds
735 name_kind = trim(particles(i)%atomic_kind%name)
736 CALL uppercase(name_kind)
737 IF (trim(colvar%gyration_param%c_kinds(j)) == name_kind) THEN
738 CALL reallocate(colvar%gyration_param%i_at, 1, colvar%gyration_param%n_atoms + 1)
739 colvar%gyration_param%n_atoms = colvar%gyration_param%n_atoms + 1
740 colvar%gyration_param%i_at(colvar%gyration_param%n_atoms) = i
741 END IF
742 END DO
743 END DO
744 stat = colvar%gyration_param%n_atoms
745 cpassert(stat /= 0)
746 END IF
747 ! Setup the colvar
748 CALL colvar_setup(colvar)
749 END IF
750 END IF
751
752 IF (colvar%type_id == rmsd_colvar_id) THEN
753 IF (colvar%rmsd_param%subset == rmsd_all) THEN
754 ! weights are masses
755 DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
756 ii = colvar%rmsd_param%i_rmsd(i)
757 colvar%rmsd_param%weights(ii) = particles(ii)%atomic_kind%mass
758 END DO
759 mtot = minval(colvar%rmsd_param%weights)
760 colvar%rmsd_param%weights = colvar%rmsd_param%weights/mtot
761 ELSE IF (colvar%rmsd_param%subset == rmsd_list) THEN
762 ! all weights are = 1
763 DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
764 ii = colvar%rmsd_param%i_rmsd(i)
765 colvar%rmsd_param%weights(ii) = 1.0_dp ! particles(ii)%atomic_kind%mass
766 END DO
767
768 END IF
769
770 IF (colvar%rmsd_param%align_frames) THEN
771 nr_frame = SIZE(colvar%rmsd_param%r_ref, 2)
772 DO i = 2, nr_frame
773 CALL rmsd3(particles, colvar%rmsd_param%r_ref(:, i), colvar%rmsd_param%r_ref(:, 1), -1, &
774 rotate=.true.)
775 END DO
776 END IF
777
778 END IF
779
780 IF (colvar%type_id == distance_from_path_colvar_id .OR. colvar%type_id == reaction_path_colvar_id) THEN
781 IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd) THEN
782 IF (colvar%reaction_path_param%align_frames) THEN
783 nr_frame = colvar%reaction_path_param%nr_frames
784 DO i = 2, nr_frame
785 CALL rmsd3(particles, colvar%reaction_path_param%r_ref(:, i), colvar%reaction_path_param%r_ref(:, 1), -1, &
786 rotate=.true.)
787 END DO
788 END IF
789 END IF
790 END IF
791
792 END SUBROUTINE post_process_colvar
793
794END MODULE colvar_utils
mathlib::invert_matrix
Definition mathlib.F:70
memory_utilities::reallocate
Definition memory_utilities.F:27
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:1630
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::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:567
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::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:513
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, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
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:88
input_constants::rmsd_all
integer, parameter, public rmsd_all
Definition input_constants.F:88
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:3362
util
All kind of helpful little routines.
Definition util.F:14
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
colvar_types::colvar_counters
Definition colvar_types.F:388
colvar_types::colvar_type
parameters for a collective variable
Definition colvar_types.F:346
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::colvar_constraint_type
Definition molecule_kind_types.F:117
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:149
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::local_colvar_constraint_type
Definition molecule_types.F:41
molecule_types::molecule_type
Definition molecule_types.F:110
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35