Line data Source code
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 198 : 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 198 : TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:)
83 : TYPE(molecule_list_type), POINTER :: molecules
84 198 : TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
85 :
86 198 : NULLIFY (subsys, molecules, molecule_kind, molecule, molecule_set, gci)
87 198 : CALL timeset(routineN, handle)
88 198 : skip_inter_colvar = .FALSE.
89 198 : my_unique = .FALSE.
90 198 : IF (PRESENT(only_intra_colvar)) skip_inter_colvar = only_intra_colvar
91 198 : IF (PRESENT(unique)) my_unique = unique
92 198 : ntot = 0
93 198 : CALL force_env_get(force_env=force_env, subsys=subsys)
94 : CALL cp_subsys_get(subsys=subsys, molecules=molecules, gci=gci, &
95 198 : molecule_kinds=molecule_kinds)
96 :
97 198 : molecule_set => molecules%els
98 : ! Intramolecular Colvar
99 198 : IF (my_unique) THEN
100 34 : molecule_kind_set => molecule_kinds%els
101 472 : DO ikind = 1, molecule_kinds%n_els
102 438 : molecule_kind => molecule_kind_set(ikind)
103 438 : CALL get_molecule_kind(molecule_kind, ncolv=ncolv)
104 472 : ntot = ntot + ncolv%ntot
105 : END DO
106 : ELSE
107 2456 : MOL: DO imol = 1, SIZE(molecule_set)
108 2292 : molecule => molecule_set(imol)
109 2292 : molecule_kind => molecule%molecule_kind
110 :
111 : CALL get_molecule_kind(molecule_kind, &
112 2292 : ncolv=ncolv)
113 2456 : ntot = ntot + ncolv%ntot
114 : END DO MOL
115 : END IF
116 : ! Intermolecular Colvar
117 198 : IF (.NOT. skip_inter_colvar) THEN
118 104 : IF (ASSOCIATED(gci)) THEN
119 104 : ntot = ntot + gci%ncolv%ntot
120 : END IF
121 : END IF
122 198 : CALL timestop(handle)
123 :
124 198 : 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 60 : 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 60 : 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 60 : NULLIFY (cell, subsys, molecule_kinds, molecule_kind, gci, colv_list)
149 60 : CALL timeset(routineN, handle)
150 60 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
151 60 : CALL cp_subsys_get(subsys=subsys, gci=gci, molecule_kinds=molecule_kinds)
152 :
153 60 : nkind = molecule_kinds%n_els
154 : ! Set Target for Intramolecular Colvars
155 120 : MOL: DO ikind = 1, nkind
156 60 : molecule_kind => molecule_kinds%els(ikind)
157 : CALL get_molecule_kind(molecule_kind, &
158 : colv_list=colv_list, &
159 60 : ncolv=ncolv)
160 120 : IF (ncolv%ntot /= 0) THEN
161 120 : DO i = 1, SIZE(colv_list)
162 60 : ind = colv_list(i)%inp_seq_num
163 120 : colv_list(i)%expected_value = targets(ind)
164 : END DO
165 : END IF
166 : END DO MOL
167 : ! Set Target for Intermolecular Colvars
168 60 : IF (ASSOCIATED(gci)) THEN
169 60 : IF (gci%ncolv%ntot /= 0) THEN
170 0 : colv_list => gci%colv_list
171 0 : DO i = 1, SIZE(colv_list)
172 0 : ind = colv_list(i)%inp_seq_num
173 0 : colv_list(i)%expected_value = targets(ind)
174 : END DO
175 : END IF
176 : END IF
177 60 : CALL timestop(handle)
178 :
179 60 : 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 94 : 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 94 : INTEGER, DIMENSION(:), POINTER :: map, wrk
205 : LOGICAL :: check
206 : REAL(KIND=dp) :: inv_error
207 94 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bwrk, Gmatrix, Gmatrix_i
208 94 : 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 94 : TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
218 : TYPE(particle_list_type), POINTER :: particles
219 94 : TYPE(particle_type), POINTER :: particle_set(:)
220 :
221 94 : NULLIFY (cell, subsys, local_molecules, molecule_kinds, &
222 94 : molecules, molecule_kind, molecule, &
223 94 : molecule_set, particles, particle_set, gci)
224 94 : IF (PRESENT(Bmatrix)) THEN
225 86 : check = ASSOCIATED(Bmatrix)
226 86 : CPASSERT(check)
227 423782 : Bmatrix = 0.0_dp
228 : END IF
229 94 : CALL timeset(routineN, handle)
230 282 : ALLOCATE (map(SIZE(cvalues)))
231 1410 : map = HUGE(0) ! init all, since used in a sort, but not all set in parallel.
232 94 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
233 94 : n_tot = 0
234 1410 : 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 94 : molecule_kinds=molecule_kinds)
241 :
242 94 : nkind = molecule_kinds%n_els
243 94 : particle_set => particles%els
244 94 : molecule_set => molecules%els
245 : ! Intramolecular Colvars
246 94 : IF (number_of_colvar(force_env, only_intra_colvar=.TRUE.) /= 0) THEN
247 474 : MOL: DO ikind = 1, nkind
248 380 : nmol_per_kind = local_molecules%n_el(ikind)
249 698 : DO imol = 1, nmol_per_kind
250 224 : i = local_molecules%list(ikind)%array(imol)
251 224 : molecule => molecule_set(i)
252 224 : molecule_kind => molecule%molecule_kind
253 :
254 : CALL get_molecule_kind(molecule_kind, &
255 224 : ncolv=ncolv)
256 224 : offset = get_colvar_offset(i, molecule_set)
257 : ! Collective variables
258 828 : IF (ncolv%ntot /= 0) THEN
259 : CALL eval_colv_int(molecule, particle_set, coords, cell, cvalues, &
260 328 : Bmatrix, offset, n_tot, map)
261 : END IF
262 : END DO
263 : END DO MOL
264 94 : CALL force_env%para_env%sum(n_tot)
265 2726 : CALL force_env%para_env%sum(cvalues)
266 847486 : IF (PRESENT(Bmatrix)) CALL force_env%para_env%sum(Bmatrix)
267 : END IF
268 94 : offset = n_tot
269 : ! Intermolecular Colvars
270 94 : IF (ASSOCIATED(gci)) THEN
271 94 : IF (gci%ncolv%ntot /= 0) THEN
272 : CALL eval_colv_ext(gci, particle_set, coords, cell, cvalues, &
273 0 : Bmatrix, offset, n_tot, map)
274 : END IF
275 : END IF
276 94 : CPASSERT(n_tot == SIZE(cvalues))
277 : ! Sort values of Collective Variables according the order of the input
278 : ! sections
279 188 : ALLOCATE (wrk(SIZE(cvalues)))
280 282 : ALLOCATE (rwrk(SIZE(cvalues)))
281 94 : CALL sort(map, SIZE(map), wrk)
282 1410 : rwrk = cvalues
283 1410 : DO i = 1, SIZE(wrk)
284 1410 : cvalues(i) = rwrk(wrk(i))
285 : END DO
286 : ! check and sort on Bmatrix
287 94 : IF (PRESENT(Bmatrix)) THEN
288 86 : check = n_tot == SIZE(Bmatrix, 2)
289 86 : CPASSERT(check)
290 344 : ALLOCATE (bwrk(SIZE(Bmatrix, 1), SIZE(Bmatrix, 2)))
291 423782 : bwrk(:, :) = Bmatrix
292 1394 : DO i = 1, SIZE(wrk)
293 423782 : Bmatrix(:, i) = bwrk(:, wrk(i))
294 : END DO
295 86 : DEALLOCATE (bwrk)
296 : END IF
297 94 : DEALLOCATE (rwrk)
298 94 : DEALLOCATE (wrk)
299 94 : DEALLOCATE (map)
300 : ! Construction of the Amatrix
301 94 : IF (PRESENT(Bmatrix) .AND. PRESENT(Amatrix)) THEN
302 60 : CPASSERT(ASSOCIATED(Amatrix))
303 60 : check = SIZE(Bmatrix, 1) == SIZE(Amatrix, 2)
304 60 : CPASSERT(check)
305 60 : check = SIZE(Bmatrix, 2) == SIZE(Amatrix, 1)
306 60 : CPASSERT(check)
307 240 : ALLOCATE (Gmatrix(n_tot, n_tot))
308 180 : ALLOCATE (Gmatrix_i(n_tot, n_tot))
309 6540 : Gmatrix(:, :) = MATMUL(TRANSPOSE(Bmatrix), Bmatrix)
310 60 : CALL invert_matrix(Gmatrix, Gmatrix_i, inv_error)
311 60 : IF (ABS(inv_error) > 1.0E-8_dp) THEN
312 0 : CPWARN("Error in inverting the Gmatrix larger than 1.0E-8!")
313 : END IF
314 30720 : Amatrix = MATMUL(Gmatrix_i, TRANSPOSE(Bmatrix))
315 60 : DEALLOCATE (Gmatrix_i)
316 120 : DEALLOCATE (Gmatrix)
317 : END IF
318 94 : IF (PRESENT(MassI)) THEN
319 86 : natom = SIZE(particle_set)
320 86 : CPASSERT(ASSOCIATED(MassI))
321 86 : CPASSERT(SIZE(MassI) == natom*3)
322 4018 : DO i = 1, natom
323 3932 : MassI((i - 1)*3 + 1) = 1.0_dp/particle_set(i)%atomic_kind%mass
324 3932 : MassI((i - 1)*3 + 2) = 1.0_dp/particle_set(i)%atomic_kind%mass
325 4018 : MassI((i - 1)*3 + 3) = 1.0_dp/particle_set(i)%atomic_kind%mass
326 : END DO
327 : END IF
328 94 : CALL timestop(handle)
329 :
330 274 : 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 224 : 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 224 : offset = 0
350 1082 : DO j = 1, i - 1
351 858 : molecule => molecule_set(j)
352 858 : molecule_kind => molecule%molecule_kind
353 : CALL get_molecule_kind(molecule_kind, &
354 858 : ncolv=ncolv)
355 1082 : offset = offset + ncolv%ntot
356 : END DO
357 :
358 224 : 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 198 : 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 198 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
387 198 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
388 198 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
389 : TYPE(molecule_kind_type), POINTER :: molecule_kind
390 :
391 198 : NULLIFY (fixd_list)
392 :
393 198 : molecule_kind => molecule%molecule_kind
394 198 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
395 198 : CALL get_molecule(molecule, lcolv=lcolv)
396 : CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
397 328 : coords, cell, cvalues, Bmatrix, offset, n_tot, map)
398 :
399 198 : 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 0 : 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 0 : colv_list => gci%colv_list
431 0 : fixd_list => gci%fixd_list
432 0 : lcolv => gci%lcolv
433 : CALL eval_colv_low(colv_list, fixd_list, lcolv, particle_set, &
434 0 : coords, cell, cvalues, Bmatrix, offset, n_tot, map)
435 :
436 0 : 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 198 : SUBROUTINE eval_colv_low(colv_list, fixd_list, lcolv, particle_set, coords, &
456 198 : 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 198 : ival = offset
473 890 : DO iconst = 1, SIZE(colv_list)
474 692 : n_tot = n_tot + 1
475 692 : ival = ival + 1
476 : ! Update colvar
477 692 : IF (PRESENT(coords)) THEN
478 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
479 204 : 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 624 : fixd_list=fixd_list)
483 : END IF
484 692 : cvalues(ival) = lcolv(iconst)%colvar%ss
485 692 : map(ival) = colv_list(iconst)%inp_seq_num
486 : ! Build the Wilson-Eliashevich Matrix
487 890 : IF (PRESENT(Bmatrix)) THEN
488 2172 : DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom)
489 1488 : ind = (lcolv(iconst)%colvar%i_atom(iatm) - 1)*3
490 1488 : Bmatrix(ind + 1, ival) = lcolv(iconst)%colvar%dsdr(1, iatm)
491 1488 : Bmatrix(ind + 2, ival) = lcolv(iconst)%colvar%dsdr(2, iatm)
492 2172 : Bmatrix(ind + 3, ival) = lcolv(iconst)%colvar%dsdr(3, iatm)
493 : END DO
494 : END IF
495 : END DO
496 :
497 198 : 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 86 : SUBROUTINE get_clv_force(force_env, forces, coords, nsize_xyz, nsize_int, cvalues, &
512 86 : 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 86 : REAL(KIND=dp), DIMENSION(:), POINTER :: MassI, wrk
522 86 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: Amatrix, Bmatrix
523 :
524 344 : ALLOCATE (Bmatrix(nsize_xyz, nsize_int))
525 258 : ALLOCATE (MassI(nsize_xyz))
526 : ! Transform gradients if requested
527 86 : IF (PRESENT(forces)) THEN
528 180 : ALLOCATE (wrk(nsize_int))
529 180 : ALLOCATE (Amatrix(nsize_int, nsize_xyz))
530 : ! Compute the transformation matrices and the inverse mass diagonal Matrix
531 60 : CALL eval_colvar(force_env, coords, cvalues, Bmatrix, MassI, Amatrix)
532 15540 : wrk = MATMUL(Amatrix, forces)
533 3120 : forces = 0.0_dp
534 120 : forces(1:nsize_int) = wrk
535 60 : DEALLOCATE (Amatrix)
536 60 : DEALLOCATE (wrk)
537 : ELSE
538 : ! Compute the transformation matrices and the inverse mass diagonal Matrix
539 52 : CALL eval_colvar(force_env, coords, cvalues, Bmatrix, MassI)
540 : END IF
541 : ! Compute the Metric Tensor
542 1394 : DO i = 1, nsize_int
543 32030 : DO j = 1, i
544 : tmp = 0.0_dp
545 10307232 : DO k = 1, nsize_xyz
546 10307232 : tmp = tmp + Bmatrix(k, j)*MassI(k)*Bmatrix(k, i)
547 : END DO
548 30636 : Mmatrix((i - 1)*nsize_int + j) = tmp
549 31944 : Mmatrix((j - 1)*nsize_int + i) = tmp
550 : END DO
551 : END DO
552 86 : DEALLOCATE (MassI)
553 86 : DEALLOCATE (Bmatrix)
554 86 : 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 912 : 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 912 : natoms = SIZE(particles)
576 912 : IF (colvar%type_id == coord_colvar_id) THEN
577 54 : IF (colvar%coord_param%use_kinds_from .OR. colvar%coord_param%use_kinds_to) THEN
578 : ! Atoms from
579 6 : IF (colvar%coord_param%use_kinds_from) THEN
580 6 : colvar%coord_param%use_kinds_from = .FALSE.
581 6 : nkinds = SIZE(colvar%coord_param%c_kinds_from)
582 30 : DO i = 1, natoms
583 54 : DO j = 1, nkinds
584 24 : name_kind = TRIM(particles(i)%atomic_kind%name)
585 24 : CALL uppercase(name_kind)
586 48 : IF (TRIM(colvar%coord_param%c_kinds_from(j)) == name_kind) THEN
587 8 : CALL reallocate(colvar%coord_param%i_at_from, 1, colvar%coord_param%n_atoms_from + 1)
588 8 : colvar%coord_param%n_atoms_from = colvar%coord_param%n_atoms_from + 1
589 8 : colvar%coord_param%i_at_from(colvar%coord_param%n_atoms_from) = i
590 : END IF
591 : END DO
592 : END DO
593 6 : stat = colvar%coord_param%n_atoms_from
594 6 : CPASSERT(stat /= 0)
595 : END IF
596 : ! Atoms to
597 6 : IF (colvar%coord_param%use_kinds_to) THEN
598 6 : colvar%coord_param%use_kinds_to = .FALSE.
599 6 : nkinds = SIZE(colvar%coord_param%c_kinds_to)
600 30 : DO i = 1, natoms
601 54 : DO j = 1, nkinds
602 24 : name_kind = TRIM(particles(i)%atomic_kind%name)
603 24 : CALL uppercase(name_kind)
604 48 : IF (TRIM(colvar%coord_param%c_kinds_to(j)) == name_kind) THEN
605 10 : CALL reallocate(colvar%coord_param%i_at_to, 1, colvar%coord_param%n_atoms_to + 1)
606 10 : colvar%coord_param%n_atoms_to = colvar%coord_param%n_atoms_to + 1
607 10 : colvar%coord_param%i_at_to(colvar%coord_param%n_atoms_to) = i
608 : END IF
609 : END DO
610 : END DO
611 6 : stat = colvar%coord_param%n_atoms_to
612 6 : CPASSERT(stat /= 0)
613 : END IF
614 : ! Atoms to b
615 6 : IF (colvar%coord_param%use_kinds_to_b) THEN
616 2 : colvar%coord_param%use_kinds_to_b = .FALSE.
617 2 : nkinds = SIZE(colvar%coord_param%c_kinds_to_b)
618 8 : DO i = 1, natoms
619 14 : DO j = 1, nkinds
620 6 : name_kind = TRIM(particles(i)%atomic_kind%name)
621 6 : CALL uppercase(name_kind)
622 12 : IF (TRIM(colvar%coord_param%c_kinds_to_b(j)) == name_kind) THEN
623 2 : CALL reallocate(colvar%coord_param%i_at_to_b, 1, colvar%coord_param%n_atoms_to_b + 1)
624 2 : colvar%coord_param%n_atoms_to_b = colvar%coord_param%n_atoms_to_b + 1
625 2 : 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 2 : stat = colvar%coord_param%n_atoms_to_b
630 2 : CPASSERT(stat /= 0)
631 : END IF
632 :
633 : ! Setup the colvar
634 6 : CALL colvar_setup(colvar)
635 : END IF
636 : END IF
637 :
638 912 : IF (colvar%type_id == mindist_colvar_id) THEN
639 0 : IF (colvar%mindist_param%use_kinds_from .OR. colvar%mindist_param%use_kinds_to) THEN
640 : ! Atoms from
641 0 : IF (colvar%mindist_param%use_kinds_from) THEN
642 0 : colvar%mindist_param%use_kinds_from = .FALSE.
643 0 : nkinds = SIZE(colvar%mindist_param%k_coord_from)
644 0 : DO i = 1, natoms
645 0 : DO j = 1, nkinds
646 0 : name_kind = TRIM(particles(i)%atomic_kind%name)
647 0 : CALL uppercase(name_kind)
648 0 : IF (TRIM(colvar%mindist_param%k_coord_from(j)) == name_kind) THEN
649 0 : CALL reallocate(colvar%mindist_param%i_coord_from, 1, colvar%mindist_param%n_coord_from + 1)
650 0 : colvar%mindist_param%n_coord_from = colvar%mindist_param%n_coord_from + 1
651 0 : colvar%mindist_param%i_coord_from(colvar%mindist_param%n_coord_from) = i
652 : END IF
653 : END DO
654 : END DO
655 0 : stat = colvar%mindist_param%n_coord_from
656 0 : CPASSERT(stat /= 0)
657 : END IF
658 : ! Atoms to
659 0 : IF (colvar%mindist_param%use_kinds_to) THEN
660 0 : colvar%mindist_param%use_kinds_to = .FALSE.
661 0 : nkinds = SIZE(colvar%mindist_param%k_coord_to)
662 0 : DO i = 1, natoms
663 0 : DO j = 1, nkinds
664 0 : name_kind = TRIM(particles(i)%atomic_kind%name)
665 0 : CALL uppercase(name_kind)
666 0 : IF (TRIM(colvar%mindist_param%k_coord_to(j)) == name_kind) THEN
667 0 : CALL reallocate(colvar%mindist_param%i_coord_to, 1, colvar%mindist_param%n_coord_to + 1)
668 0 : colvar%mindist_param%n_coord_to = colvar%mindist_param%n_coord_to + 1
669 0 : colvar%mindist_param%i_coord_to(colvar%mindist_param%n_coord_to) = i
670 : END IF
671 : END DO
672 : END DO
673 0 : stat = colvar%mindist_param%n_coord_to
674 0 : CPASSERT(stat /= 0)
675 : END IF
676 : ! Setup the colvar
677 0 : CALL colvar_setup(colvar)
678 : END IF
679 : END IF
680 :
681 912 : IF (colvar%type_id == population_colvar_id) THEN
682 :
683 8 : IF (colvar%population_param%use_kinds_from .OR. colvar%population_param%use_kinds_to) THEN
684 : ! Atoms from
685 8 : IF (colvar%population_param%use_kinds_from) THEN
686 0 : colvar%population_param%use_kinds_from = .FALSE.
687 0 : nkinds = SIZE(colvar%population_param%c_kinds_from)
688 0 : DO i = 1, natoms
689 0 : DO j = 1, nkinds
690 0 : name_kind = TRIM(particles(i)%atomic_kind%name)
691 0 : CALL uppercase(name_kind)
692 0 : IF (TRIM(colvar%population_param%c_kinds_from(j)) == name_kind) THEN
693 0 : CALL reallocate(colvar%population_param%i_at_from, 1, colvar%population_param%n_atoms_from + 1)
694 0 : colvar%population_param%n_atoms_from = colvar%population_param%n_atoms_from + 1
695 0 : colvar%population_param%i_at_from(colvar%population_param%n_atoms_from) = i
696 : END IF
697 : END DO
698 : END DO
699 0 : stat = colvar%population_param%n_atoms_from
700 0 : CPASSERT(stat /= 0)
701 : END IF
702 : ! Atoms to
703 8 : IF (colvar%population_param%use_kinds_to) THEN
704 8 : colvar%population_param%use_kinds_to = .FALSE.
705 8 : nkinds = SIZE(colvar%population_param%c_kinds_to)
706 32 : DO i = 1, natoms
707 56 : DO j = 1, nkinds
708 24 : name_kind = TRIM(particles(i)%atomic_kind%name)
709 24 : CALL uppercase(name_kind)
710 48 : IF (TRIM(colvar%population_param%c_kinds_to(j)) == name_kind) THEN
711 16 : CALL reallocate(colvar%population_param%i_at_to, 1, colvar%population_param%n_atoms_to + 1)
712 16 : colvar%population_param%n_atoms_to = colvar%population_param%n_atoms_to + 1
713 16 : colvar%population_param%i_at_to(colvar%population_param%n_atoms_to) = i
714 : END IF
715 : END DO
716 : END DO
717 8 : stat = colvar%population_param%n_atoms_to
718 8 : CPASSERT(stat /= 0)
719 : END IF
720 : ! Setup the colvar
721 8 : CALL colvar_setup(colvar)
722 : END IF
723 :
724 : END IF
725 :
726 912 : IF (colvar%type_id == gyration_colvar_id) THEN
727 :
728 2 : IF (colvar%gyration_param%use_kinds) THEN
729 : ! Atoms from
730 : IF (colvar%gyration_param%use_kinds) THEN
731 2 : colvar%gyration_param%use_kinds = .FALSE.
732 2 : nkinds = SIZE(colvar%gyration_param%c_kinds)
733 28 : DO i = 1, natoms
734 54 : DO j = 1, nkinds
735 26 : name_kind = TRIM(particles(i)%atomic_kind%name)
736 26 : CALL uppercase(name_kind)
737 52 : IF (TRIM(colvar%gyration_param%c_kinds(j)) == name_kind) THEN
738 26 : CALL reallocate(colvar%gyration_param%i_at, 1, colvar%gyration_param%n_atoms + 1)
739 26 : colvar%gyration_param%n_atoms = colvar%gyration_param%n_atoms + 1
740 26 : colvar%gyration_param%i_at(colvar%gyration_param%n_atoms) = i
741 : END IF
742 : END DO
743 : END DO
744 2 : stat = colvar%gyration_param%n_atoms
745 2 : CPASSERT(stat /= 0)
746 : END IF
747 : ! Setup the colvar
748 2 : CALL colvar_setup(colvar)
749 : END IF
750 : END IF
751 :
752 912 : IF (colvar%type_id == rmsd_colvar_id) THEN
753 4 : IF (colvar%rmsd_param%subset == rmsd_all) THEN
754 : ! weights are masses
755 0 : DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
756 0 : ii = colvar%rmsd_param%i_rmsd(i)
757 0 : colvar%rmsd_param%weights(ii) = particles(ii)%atomic_kind%mass
758 : END DO
759 0 : mtot = MINVAL(colvar%rmsd_param%weights)
760 0 : colvar%rmsd_param%weights = colvar%rmsd_param%weights/mtot
761 4 : ELSE IF (colvar%rmsd_param%subset == rmsd_list) THEN
762 : ! all weights are = 1
763 28 : DO i = 1, SIZE(colvar%rmsd_param%i_rmsd)
764 24 : ii = colvar%rmsd_param%i_rmsd(i)
765 28 : colvar%rmsd_param%weights(ii) = 1.0_dp ! particles(ii)%atomic_kind%mass
766 : END DO
767 :
768 : END IF
769 :
770 4 : IF (colvar%rmsd_param%align_frames) THEN
771 0 : nr_frame = SIZE(colvar%rmsd_param%r_ref, 2)
772 0 : DO i = 2, nr_frame
773 : CALL rmsd3(particles, colvar%rmsd_param%r_ref(:, i), colvar%rmsd_param%r_ref(:, 1), -1, &
774 0 : rotate=.TRUE.)
775 : END DO
776 : END IF
777 :
778 : END IF
779 :
780 912 : IF (colvar%type_id == distance_from_path_colvar_id .OR. colvar%type_id == reaction_path_colvar_id) THEN
781 20 : IF (colvar%reaction_path_param%dist_rmsd .OR. colvar%reaction_path_param%rmsd) THEN
782 8 : IF (colvar%reaction_path_param%align_frames) THEN
783 8 : nr_frame = colvar%reaction_path_param%nr_frames
784 40 : 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 40 : rotate=.TRUE.)
787 : END DO
788 : END IF
789 : END IF
790 : END IF
791 :
792 912 : END SUBROUTINE post_process_colvar
793 :
794 60 : END MODULE colvar_utils
|
