Line data Source code
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 Define the molecule kind structure types and the corresponding
10 : !> functionality
11 : !> \par History
12 : !> Teodoro Laino [tlaino] 12.2008 - Preparing for VIRTUAL SITE constraints
13 : !> (patch by Marcel Baer)
14 : !> \author Matthias Krack (22.08.2003)
15 : ! **************************************************************************************************
16 : MODULE molecule_kind_types
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE cell_types, ONLY: periodicity_string,&
20 : use_perd_x,&
21 : use_perd_xy,&
22 : use_perd_xyz,&
23 : use_perd_xz,&
24 : use_perd_y,&
25 : use_perd_yz,&
26 : use_perd_z
27 : USE colvar_types, ONLY: &
28 : Wc_colvar_id, acid_hyd_dist_colvar_id, acid_hyd_shell_colvar_id, angle_colvar_id, &
29 : colvar_counters, combine_colvar_id, coord_colvar_id, dfunct_colvar_id, dist_colvar_id, &
30 : distance_from_path_colvar_id, gyration_colvar_id, hbp_colvar_id, hydronium_dist_colvar_id, &
31 : hydronium_shell_colvar_id, mindist_colvar_id, no_colvar_id, plane_distance_colvar_id, &
32 : plane_plane_angle_colvar_id, population_colvar_id, qparm_colvar_id, &
33 : reaction_path_colvar_id, ring_puckering_colvar_id, rmsd_colvar_id, rotation_colvar_id, &
34 : torsion_colvar_id, u_colvar_id, xyz_diag_colvar_id, xyz_outerdiag_colvar_id
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_type
37 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
38 : cp_print_key_unit_nr
39 : USE cp_units, ONLY: cp_unit_from_cp2k
40 : USE force_field_kind_types, ONLY: &
41 : bend_kind_type, bond_kind_type, do_ff_undef, impr_kind_dealloc_ref, impr_kind_type, &
42 : opbend_kind_type, torsion_kind_dealloc_ref, torsion_kind_type, ub_kind_dealloc_ref, &
43 : ub_kind_type
44 : USE input_section_types, ONLY: section_vals_type
45 : USE kinds, ONLY: default_string_length,&
46 : dp
47 : USE shell_potential_types, ONLY: shell_kind_type
48 : #include "../base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 :
52 : PRIVATE
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecule_kind_types'
55 :
56 : ! Define the derived structure types
57 :
58 : TYPE atom_type
59 : TYPE(atomic_kind_type), POINTER :: atomic_kind => NULL()
60 : INTEGER :: id_name = 0
61 : END TYPE atom_type
62 :
63 : TYPE shell_type
64 : INTEGER :: a = 0
65 : CHARACTER(LEN=default_string_length) :: name = ""
66 : TYPE(shell_kind_type), POINTER :: shell_kind => NULL()
67 : END TYPE shell_type
68 :
69 : TYPE bond_type
70 : INTEGER :: a = 0, b = 0
71 : INTEGER :: id_type = do_ff_undef, itype = 0
72 : TYPE(bond_kind_type), POINTER :: bond_kind => NULL()
73 : END TYPE bond_type
74 :
75 : TYPE bend_type
76 : INTEGER :: a = 0, b = 0, c = 0
77 : INTEGER :: id_type = do_ff_undef, itype = 0
78 : TYPE(bend_kind_type), POINTER :: bend_kind => NULL()
79 : END TYPE bend_type
80 :
81 : TYPE ub_type
82 : INTEGER :: a = 0, b = 0, c = 0
83 : INTEGER :: id_type = do_ff_undef, itype = 0
84 : TYPE(ub_kind_type), POINTER :: ub_kind => NULL()
85 : END TYPE ub_type
86 :
87 : TYPE torsion_type
88 : INTEGER :: a = 0, b = 0, c = 0, d = 0
89 : INTEGER :: id_type = do_ff_undef, itype = 0
90 : TYPE(torsion_kind_type), POINTER :: torsion_kind => NULL()
91 : END TYPE torsion_type
92 :
93 : TYPE impr_type
94 : INTEGER :: a = 0, b = 0, c = 0, d = 0
95 : INTEGER :: id_type = do_ff_undef, itype = 0
96 : TYPE(impr_kind_type), POINTER :: impr_kind => NULL()
97 : END TYPE impr_type
98 :
99 : TYPE opbend_type
100 : INTEGER :: a = 0, b = 0, c = 0, d = 0
101 : INTEGER :: id_type = do_ff_undef, itype = 0
102 : TYPE(opbend_kind_type), POINTER :: opbend_kind => NULL()
103 : END TYPE opbend_type
104 :
105 : TYPE restraint_type
106 : LOGICAL :: active = .FALSE.
107 : REAL(KIND=dp) :: k0 = 0.0_dp
108 : END TYPE restraint_type
109 :
110 : ! Constraint types
111 : TYPE colvar_constraint_type
112 : INTEGER :: type_id = no_colvar_id
113 : INTEGER :: inp_seq_num = 0
114 : LOGICAL :: use_points = .FALSE.
115 : REAL(KIND=dp) :: expected_value = 0.0_dp
116 : REAL(KIND=dp) :: expected_value_growth_speed = 0.0_dp
117 : INTEGER, POINTER, DIMENSION(:) :: i_atoms => NULL()
118 : TYPE(restraint_type) :: restraint = restraint_type()
119 : END TYPE colvar_constraint_type
120 :
121 : TYPE g3x3_constraint_type
122 : INTEGER :: a = 0, b = 0, c = 0
123 : REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp
124 : TYPE(restraint_type) :: restraint = restraint_type()
125 : END TYPE g3x3_constraint_type
126 :
127 : TYPE g4x6_constraint_type
128 : INTEGER :: a = 0, b = 0, c = 0, d = 0
129 : REAL(KIND=dp) :: dab = 0.0_dp, dac = 0.0_dp, dbc = 0.0_dp, &
130 : dad = 0.0_dp, dbd = 0.0_dp, dcd = 0.0_dp
131 : TYPE(restraint_type) :: restraint = restraint_type()
132 : END TYPE g4x6_constraint_type
133 :
134 : TYPE vsite_constraint_type
135 : INTEGER :: a = 0, b = 0, c = 0, d = 0
136 : REAL(KIND=dp) :: wbc = 0.0_dp, wdc = 0.0_dp
137 : TYPE(restraint_type) :: restraint = restraint_type()
138 : END TYPE vsite_constraint_type
139 :
140 : TYPE fixd_constraint_type
141 : TYPE(restraint_type) :: restraint = restraint_type()
142 : INTEGER :: fixd = 0, itype = 0
143 : REAL(KIND=dp), DIMENSION(3) :: coord = 0.0_dp
144 : END TYPE fixd_constraint_type
145 :
146 : TYPE local_fixd_constraint_type
147 : INTEGER :: ifixd_index = 0, ikind = 0
148 : END TYPE local_fixd_constraint_type
149 :
150 : ! Molecule kind type
151 : TYPE molecule_kind_type
152 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list => NULL()
153 : TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set => NULL()
154 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list => NULL()
155 : TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set => NULL()
156 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list => NULL()
157 : TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set => NULL()
158 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list => NULL()
159 : TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set => NULL()
160 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list => NULL()
161 : TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set => NULL()
162 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list => NULL()
163 : TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set => NULL()
164 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list => NULL()
165 : TYPE(colvar_constraint_type), DIMENSION(:), &
166 : POINTER :: colv_list => NULL()
167 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list => NULL()
168 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list => NULL()
169 : TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list => NULL()
170 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list => NULL()
171 : TYPE(shell_type), DIMENSION(:), POINTER :: shell_list => NULL()
172 : CHARACTER(LEN=default_string_length) :: name = ""
173 : REAL(KIND=dp) :: charge = 0.0_dp, &
174 : mass = 0.0_dp
175 : INTEGER :: kind_number = 0, &
176 : natom = 0, &
177 : nbond = 0, &
178 : nbend = 0, &
179 : nimpr = 0, &
180 : nopbend = 0, &
181 : ntorsion = 0, &
182 : nub = 0, &
183 : ng3x3 = 0, &
184 : ng3x3_restraint = 0, &
185 : ng4x6 = 0, &
186 : ng4x6_restraint = 0, &
187 : nvsite = 0, &
188 : nvsite_restraint = 0, &
189 : nfixd = 0, &
190 : nfixd_restraint = 0, &
191 : nmolecule = 0, &
192 : nshell = 0
193 : TYPE(colvar_counters) :: ncolv = colvar_counters()
194 : INTEGER :: nsgf = 0, &
195 : nelectron = 0, &
196 : nelectron_alpha = 0, &
197 : nelectron_beta = 0
198 : INTEGER, DIMENSION(:), POINTER :: molecule_list => NULL()
199 : LOGICAL :: molname_generated = .FALSE.
200 : END TYPE molecule_kind_type
201 :
202 : ! Public subroutines
203 : PUBLIC :: allocate_molecule_kind_set, &
204 : deallocate_molecule_kind_set, &
205 : get_molecule_kind, &
206 : get_molecule_kind_set, &
207 : set_molecule_kind, &
208 : write_molecule_kind_set, &
209 : setup_colvar_counters, &
210 : write_colvar_constraint, &
211 : write_fixd_constraint, &
212 : write_g3x3_constraint, &
213 : write_g4x6_constraint, &
214 : write_vsite_constraint
215 :
216 : ! Public data types
217 : PUBLIC :: atom_type, &
218 : bend_type, &
219 : bond_type, &
220 : ub_type, &
221 : torsion_type, &
222 : impr_type, &
223 : opbend_type, &
224 : colvar_constraint_type, &
225 : g3x3_constraint_type, &
226 : g4x6_constraint_type, &
227 : vsite_constraint_type, &
228 : fixd_constraint_type, &
229 : local_fixd_constraint_type, &
230 : molecule_kind_type, &
231 : shell_type
232 :
233 : CONTAINS
234 :
235 : ! **************************************************************************************************
236 : !> \brief ...
237 : !> \param colv_list ...
238 : !> \param ncolv ...
239 : ! **************************************************************************************************
240 146795 : SUBROUTINE setup_colvar_counters(colv_list, ncolv)
241 : TYPE(colvar_constraint_type), DIMENSION(:), &
242 : POINTER :: colv_list
243 : TYPE(colvar_counters), INTENT(OUT) :: ncolv
244 :
245 : INTEGER :: k
246 :
247 146795 : IF (ASSOCIATED(colv_list)) THEN
248 1070 : DO k = 1, SIZE(colv_list)
249 448 : IF (colv_list(k)%restraint%active) ncolv%nrestraint = ncolv%nrestraint + 1
250 622 : SELECT CASE (colv_list(k)%type_id)
251 : CASE (angle_colvar_id)
252 50 : ncolv%nangle = ncolv%nangle + 1
253 : CASE (coord_colvar_id)
254 2 : ncolv%ncoord = ncolv%ncoord + 1
255 : CASE (population_colvar_id)
256 0 : ncolv%npopulation = ncolv%npopulation + 1
257 : CASE (gyration_colvar_id)
258 0 : ncolv%ngyration = ncolv%ngyration + 1
259 : CASE (rotation_colvar_id)
260 0 : ncolv%nrot = ncolv%nrot + 1
261 : CASE (dist_colvar_id)
262 334 : ncolv%ndist = ncolv%ndist + 1
263 : CASE (dfunct_colvar_id)
264 4 : ncolv%ndfunct = ncolv%ndfunct + 1
265 : CASE (plane_distance_colvar_id)
266 0 : ncolv%nplane_dist = ncolv%nplane_dist + 1
267 : CASE (plane_plane_angle_colvar_id)
268 4 : ncolv%nplane_angle = ncolv%nplane_angle + 1
269 : CASE (torsion_colvar_id)
270 38 : ncolv%ntorsion = ncolv%ntorsion + 1
271 : CASE (qparm_colvar_id)
272 0 : ncolv%nqparm = ncolv%nqparm + 1
273 : CASE (xyz_diag_colvar_id)
274 6 : ncolv%nxyz_diag = ncolv%nxyz_diag + 1
275 : CASE (xyz_outerdiag_colvar_id)
276 6 : ncolv%nxyz_outerdiag = ncolv%nxyz_outerdiag + 1
277 : CASE (hydronium_shell_colvar_id)
278 0 : ncolv%nhydronium_shell = ncolv%nhydronium_shell + 1
279 : CASE (hydronium_dist_colvar_id)
280 0 : ncolv%nhydronium_dist = ncolv%nhydronium_dist + 1
281 : CASE (acid_hyd_dist_colvar_id)
282 0 : ncolv%nacid_hyd_dist = ncolv%nacid_hyd_dist + 1
283 : CASE (acid_hyd_shell_colvar_id)
284 0 : ncolv%nacid_hyd_shell = ncolv%nacid_hyd_shell + 1
285 : CASE (reaction_path_colvar_id)
286 2 : ncolv%nreactionpath = ncolv%nreactionpath + 1
287 : CASE (combine_colvar_id)
288 2 : ncolv%ncombinecvs = ncolv%ncombinecvs + 1
289 : CASE DEFAULT
290 448 : CPABORT("")
291 : END SELECT
292 : END DO
293 : END IF
294 : ncolv%ntot = ncolv%ndist + &
295 : ncolv%nangle + &
296 : ncolv%ntorsion + &
297 : ncolv%ncoord + &
298 : ncolv%nplane_dist + &
299 : ncolv%nplane_angle + &
300 : ncolv%ndfunct + &
301 : ncolv%nrot + &
302 : ncolv%nqparm + &
303 : ncolv%nxyz_diag + &
304 : ncolv%nxyz_outerdiag + &
305 : ncolv%nhydronium_shell + &
306 : ncolv%nhydronium_dist + &
307 : ncolv%nacid_hyd_dist + &
308 : ncolv%nacid_hyd_shell + &
309 : ncolv%nreactionpath + &
310 : ncolv%ncombinecvs + &
311 : ncolv%npopulation + &
312 146795 : ncolv%ngyration
313 :
314 146795 : END SUBROUTINE setup_colvar_counters
315 :
316 : ! **************************************************************************************************
317 : !> \brief Allocate and initialize a molecule kind set.
318 : !> \param molecule_kind_set ...
319 : !> \param nmolecule_kind ...
320 : !> \date 22.08.2003
321 : !> \author Matthias Krack
322 : !> \version 1.0
323 : ! **************************************************************************************************
324 10278 : SUBROUTINE allocate_molecule_kind_set(molecule_kind_set, nmolecule_kind)
325 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
326 : INTEGER, INTENT(IN) :: nmolecule_kind
327 :
328 : INTEGER :: imolecule_kind
329 :
330 10278 : IF (ASSOCIATED(molecule_kind_set)) THEN
331 0 : CALL deallocate_molecule_kind_set(molecule_kind_set)
332 : END IF
333 :
334 167267 : ALLOCATE (molecule_kind_set(nmolecule_kind))
335 :
336 146711 : DO imolecule_kind = 1, nmolecule_kind
337 136433 : molecule_kind_set(imolecule_kind)%kind_number = imolecule_kind
338 : CALL setup_colvar_counters(molecule_kind_set(imolecule_kind)%colv_list, &
339 146711 : molecule_kind_set(imolecule_kind)%ncolv)
340 : END DO
341 :
342 10278 : END SUBROUTINE allocate_molecule_kind_set
343 :
344 : ! **************************************************************************************************
345 : !> \brief Deallocate a molecule kind set.
346 : !> \param molecule_kind_set ...
347 : !> \date 22.08.2003
348 : !> \author Matthias Krack
349 : !> \version 1.0
350 : ! **************************************************************************************************
351 10278 : SUBROUTINE deallocate_molecule_kind_set(molecule_kind_set)
352 :
353 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
354 :
355 : INTEGER :: i, imolecule_kind, j, nmolecule_kind
356 :
357 10278 : IF (ASSOCIATED(molecule_kind_set)) THEN
358 :
359 10278 : nmolecule_kind = SIZE(molecule_kind_set)
360 :
361 146711 : DO imolecule_kind = 1, nmolecule_kind
362 :
363 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%atom_list)) THEN
364 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%atom_list)
365 : END IF
366 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bend_kind_set)) THEN
367 122745 : DO i = 1, SIZE(molecule_kind_set(imolecule_kind)%bend_kind_set)
368 93644 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bend_kind_set(i)%legendre%coeffs)) &
369 31172 : DEALLOCATE (molecule_kind_set(imolecule_kind)%bend_kind_set(i)%legendre%coeffs)
370 : END DO
371 29101 : DEALLOCATE (molecule_kind_set(imolecule_kind)%bend_kind_set)
372 : END IF
373 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bend_list)) THEN
374 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%bend_list)
375 : END IF
376 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%ub_list)) THEN
377 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%ub_list)
378 : END IF
379 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%ub_kind_set)) THEN
380 29087 : CALL ub_kind_dealloc_ref(molecule_kind_set(imolecule_kind)%ub_kind_set)
381 : END IF
382 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%impr_list)) THEN
383 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%impr_list)
384 : END IF
385 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%impr_kind_set)) THEN
386 4882 : DO i = 1, SIZE(molecule_kind_set(imolecule_kind)%impr_kind_set)
387 4882 : CALL impr_kind_dealloc_ref() !This Subroutine doesn't deallocate anything, maybe needs to be implemented
388 : END DO
389 1672 : DEALLOCATE (molecule_kind_set(imolecule_kind)%impr_kind_set)
390 : END IF
391 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_list)) THEN
392 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_list)
393 : END IF
394 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%opbend_kind_set)) THEN
395 1672 : DEALLOCATE (molecule_kind_set(imolecule_kind)%opbend_kind_set)
396 : END IF
397 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bond_kind_set)) THEN
398 29431 : DEALLOCATE (molecule_kind_set(imolecule_kind)%bond_kind_set)
399 : END IF
400 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%bond_list)) THEN
401 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%bond_list)
402 : END IF
403 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%colv_list)) THEN
404 960 : DO j = 1, SIZE(molecule_kind_set(imolecule_kind)%colv_list)
405 960 : DEALLOCATE (molecule_kind_set(imolecule_kind)%colv_list(j)%i_atoms)
406 : END DO
407 578 : DEALLOCATE (molecule_kind_set(imolecule_kind)%colv_list)
408 : END IF
409 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%g3x3_list)) THEN
410 270 : DEALLOCATE (molecule_kind_set(imolecule_kind)%g3x3_list)
411 : END IF
412 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%g4x6_list)) THEN
413 20 : DEALLOCATE (molecule_kind_set(imolecule_kind)%g4x6_list)
414 : END IF
415 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%vsite_list)) THEN
416 10 : DEALLOCATE (molecule_kind_set(imolecule_kind)%vsite_list)
417 : END IF
418 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%fixd_list)) THEN
419 4926 : DEALLOCATE (molecule_kind_set(imolecule_kind)%fixd_list)
420 : END IF
421 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%torsion_kind_set)) THEN
422 83963 : DO i = 1, SIZE(molecule_kind_set(imolecule_kind)%torsion_kind_set)
423 83963 : CALL torsion_kind_dealloc_ref(molecule_kind_set(imolecule_kind)%torsion_kind_set(i))
424 : END DO
425 5534 : DEALLOCATE (molecule_kind_set(imolecule_kind)%torsion_kind_set)
426 : END IF
427 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%shell_list)) THEN
428 10874 : DEALLOCATE (molecule_kind_set(imolecule_kind)%shell_list)
429 : END IF
430 136433 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%torsion_list)) THEN
431 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%torsion_list)
432 : END IF
433 146711 : IF (ASSOCIATED(molecule_kind_set(imolecule_kind)%molecule_list)) THEN
434 136433 : DEALLOCATE (molecule_kind_set(imolecule_kind)%molecule_list)
435 : END IF
436 : END DO
437 :
438 10278 : DEALLOCATE (molecule_kind_set)
439 : END IF
440 10278 : NULLIFY (molecule_kind_set)
441 :
442 10278 : END SUBROUTINE deallocate_molecule_kind_set
443 :
444 : ! **************************************************************************************************
445 : !> \brief Get informations about a molecule kind.
446 : !> \param molecule_kind ...
447 : !> \param atom_list ...
448 : !> \param bond_list ...
449 : !> \param bend_list ...
450 : !> \param ub_list ...
451 : !> \param impr_list ...
452 : !> \param opbend_list ...
453 : !> \param colv_list ...
454 : !> \param fixd_list ...
455 : !> \param g3x3_list ...
456 : !> \param g4x6_list ...
457 : !> \param vsite_list ...
458 : !> \param torsion_list ...
459 : !> \param shell_list ...
460 : !> \param name ...
461 : !> \param mass ...
462 : !> \param charge ...
463 : !> \param kind_number ...
464 : !> \param natom ...
465 : !> \param nbend ...
466 : !> \param nbond ...
467 : !> \param nub ...
468 : !> \param nimpr ...
469 : !> \param nopbend ...
470 : !> \param nconstraint ...
471 : !> \param nconstraint_fixd ...
472 : !> \param nfixd ...
473 : !> \param ncolv ...
474 : !> \param ng3x3 ...
475 : !> \param ng4x6 ...
476 : !> \param nvsite ...
477 : !> \param nfixd_restraint ...
478 : !> \param ng3x3_restraint ...
479 : !> \param ng4x6_restraint ...
480 : !> \param nvsite_restraint ...
481 : !> \param nrestraints ...
482 : !> \param nmolecule ...
483 : !> \param nsgf ...
484 : !> \param nshell ...
485 : !> \param ntorsion ...
486 : !> \param molecule_list ...
487 : !> \param nelectron ...
488 : !> \param nelectron_alpha ...
489 : !> \param nelectron_beta ...
490 : !> \param bond_kind_set ...
491 : !> \param bend_kind_set ...
492 : !> \param ub_kind_set ...
493 : !> \param impr_kind_set ...
494 : !> \param opbend_kind_set ...
495 : !> \param torsion_kind_set ...
496 : !> \param molname_generated ...
497 : !> \date 27.08.2003
498 : !> \author Matthias Krack
499 : !> \version 1.0
500 : ! **************************************************************************************************
501 16466699 : SUBROUTINE get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, &
502 : ub_list, impr_list, opbend_list, colv_list, fixd_list, &
503 : g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, &
504 : name, mass, charge, kind_number, natom, nbend, nbond, nub, &
505 : nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, &
506 : nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, &
507 : nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, &
508 : molecule_list, nelectron, nelectron_alpha, nelectron_beta, &
509 : bond_kind_set, bend_kind_set, &
510 : ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, &
511 : molname_generated)
512 :
513 : TYPE(molecule_kind_type), INTENT(IN) :: molecule_kind
514 : TYPE(atom_type), DIMENSION(:), OPTIONAL, POINTER :: atom_list
515 : TYPE(bond_type), DIMENSION(:), OPTIONAL, POINTER :: bond_list
516 : TYPE(bend_type), DIMENSION(:), OPTIONAL, POINTER :: bend_list
517 : TYPE(ub_type), DIMENSION(:), OPTIONAL, POINTER :: ub_list
518 : TYPE(impr_type), DIMENSION(:), OPTIONAL, POINTER :: impr_list
519 : TYPE(opbend_type), DIMENSION(:), OPTIONAL, POINTER :: opbend_list
520 : TYPE(colvar_constraint_type), DIMENSION(:), &
521 : OPTIONAL, POINTER :: colv_list
522 : TYPE(fixd_constraint_type), DIMENSION(:), &
523 : OPTIONAL, POINTER :: fixd_list
524 : TYPE(g3x3_constraint_type), DIMENSION(:), &
525 : OPTIONAL, POINTER :: g3x3_list
526 : TYPE(g4x6_constraint_type), DIMENSION(:), &
527 : OPTIONAL, POINTER :: g4x6_list
528 : TYPE(vsite_constraint_type), DIMENSION(:), &
529 : OPTIONAL, POINTER :: vsite_list
530 : TYPE(torsion_type), DIMENSION(:), OPTIONAL, &
531 : POINTER :: torsion_list
532 : TYPE(shell_type), DIMENSION(:), OPTIONAL, POINTER :: shell_list
533 : CHARACTER(LEN=default_string_length), &
534 : INTENT(OUT), OPTIONAL :: name
535 : REAL(KIND=dp), OPTIONAL :: mass, charge
536 : INTEGER, INTENT(OUT), OPTIONAL :: kind_number, natom, nbend, nbond, nub, &
537 : nimpr, nopbend, nconstraint, &
538 : nconstraint_fixd, nfixd
539 : TYPE(colvar_counters), INTENT(out), OPTIONAL :: ncolv
540 : INTEGER, INTENT(OUT), OPTIONAL :: ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, &
541 : ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion
542 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list
543 : INTEGER, INTENT(OUT), OPTIONAL :: nelectron, nelectron_alpha, &
544 : nelectron_beta
545 : TYPE(bond_kind_type), DIMENSION(:), OPTIONAL, &
546 : POINTER :: bond_kind_set
547 : TYPE(bend_kind_type), DIMENSION(:), OPTIONAL, &
548 : POINTER :: bend_kind_set
549 : TYPE(ub_kind_type), DIMENSION(:), OPTIONAL, &
550 : POINTER :: ub_kind_set
551 : TYPE(impr_kind_type), DIMENSION(:), OPTIONAL, &
552 : POINTER :: impr_kind_set
553 : TYPE(opbend_kind_type), DIMENSION(:), OPTIONAL, &
554 : POINTER :: opbend_kind_set
555 : TYPE(torsion_kind_type), DIMENSION(:), OPTIONAL, &
556 : POINTER :: torsion_kind_set
557 : LOGICAL, INTENT(OUT), OPTIONAL :: molname_generated
558 :
559 : INTEGER :: i
560 :
561 16466699 : IF (PRESENT(atom_list)) atom_list => molecule_kind%atom_list
562 16466699 : IF (PRESENT(bend_list)) bend_list => molecule_kind%bend_list
563 16466699 : IF (PRESENT(bond_list)) bond_list => molecule_kind%bond_list
564 16466699 : IF (PRESENT(impr_list)) impr_list => molecule_kind%impr_list
565 16466699 : IF (PRESENT(opbend_list)) opbend_list => molecule_kind%opbend_list
566 16466699 : IF (PRESENT(ub_list)) ub_list => molecule_kind%ub_list
567 16466699 : IF (PRESENT(bond_kind_set)) bond_kind_set => molecule_kind%bond_kind_set
568 16466699 : IF (PRESENT(bend_kind_set)) bend_kind_set => molecule_kind%bend_kind_set
569 16466699 : IF (PRESENT(ub_kind_set)) ub_kind_set => molecule_kind%ub_kind_set
570 16466699 : IF (PRESENT(impr_kind_set)) impr_kind_set => molecule_kind%impr_kind_set
571 16466699 : IF (PRESENT(opbend_kind_set)) opbend_kind_set => molecule_kind%opbend_kind_set
572 16466699 : IF (PRESENT(torsion_kind_set)) torsion_kind_set => molecule_kind%torsion_kind_set
573 16466699 : IF (PRESENT(colv_list)) colv_list => molecule_kind%colv_list
574 16466699 : IF (PRESENT(g3x3_list)) g3x3_list => molecule_kind%g3x3_list
575 16466699 : IF (PRESENT(g4x6_list)) g4x6_list => molecule_kind%g4x6_list
576 16466699 : IF (PRESENT(vsite_list)) vsite_list => molecule_kind%vsite_list
577 16466699 : IF (PRESENT(fixd_list)) fixd_list => molecule_kind%fixd_list
578 16466699 : IF (PRESENT(torsion_list)) torsion_list => molecule_kind%torsion_list
579 16466699 : IF (PRESENT(shell_list)) shell_list => molecule_kind%shell_list
580 16466699 : IF (PRESENT(name)) name = molecule_kind%name
581 16466699 : IF (PRESENT(molname_generated)) molname_generated = molecule_kind%molname_generated
582 16466699 : IF (PRESENT(mass)) mass = molecule_kind%mass
583 16466699 : IF (PRESENT(charge)) charge = molecule_kind%charge
584 16466699 : IF (PRESENT(kind_number)) kind_number = molecule_kind%kind_number
585 16466699 : IF (PRESENT(natom)) natom = molecule_kind%natom
586 16466699 : IF (PRESENT(nbend)) nbend = molecule_kind%nbend
587 16466699 : IF (PRESENT(nbond)) nbond = molecule_kind%nbond
588 16466699 : IF (PRESENT(nub)) nub = molecule_kind%nub
589 16466699 : IF (PRESENT(nimpr)) nimpr = molecule_kind%nimpr
590 16466699 : IF (PRESENT(nopbend)) nopbend = molecule_kind%nopbend
591 16466699 : IF (PRESENT(nconstraint)) nconstraint = (molecule_kind%ncolv%ntot - molecule_kind%ncolv%nrestraint) + &
592 : 3*(molecule_kind%ng3x3 - molecule_kind%ng3x3_restraint) + &
593 : 6*(molecule_kind%ng4x6 - molecule_kind%ng4x6_restraint) + &
594 3180230 : 3*(molecule_kind%nvsite - molecule_kind%nvsite_restraint)
595 16466699 : IF (PRESENT(ncolv)) ncolv = molecule_kind%ncolv
596 16466699 : IF (PRESENT(ng3x3)) ng3x3 = molecule_kind%ng3x3
597 16466699 : IF (PRESENT(ng4x6)) ng4x6 = molecule_kind%ng4x6
598 16466699 : IF (PRESENT(nvsite)) nvsite = molecule_kind%nvsite
599 : ! Number of atoms that have one or more components fixed
600 16466699 : IF (PRESENT(nfixd)) nfixd = molecule_kind%nfixd
601 : ! Number of degrees of freedom fixed
602 16466699 : IF (PRESENT(nconstraint_fixd)) THEN
603 279728 : nconstraint_fixd = 0
604 279728 : IF (molecule_kind%nfixd /= 0) THEN
605 171910 : DO i = 1, SIZE(molecule_kind%fixd_list)
606 170016 : IF (molecule_kind%fixd_list(i)%restraint%active) CYCLE
607 1894 : SELECT CASE (molecule_kind%fixd_list(i)%itype)
608 : CASE (use_perd_x, use_perd_y, use_perd_z)
609 62976 : nconstraint_fixd = nconstraint_fixd + 1
610 : CASE (use_perd_xy, use_perd_xz, use_perd_yz)
611 20992 : nconstraint_fixd = nconstraint_fixd + 2
612 : CASE (use_perd_xyz)
613 169588 : nconstraint_fixd = nconstraint_fixd + 3
614 : END SELECT
615 : END DO
616 : END IF
617 : END IF
618 16466699 : IF (PRESENT(ng3x3_restraint)) ng3x3_restraint = molecule_kind%ng3x3_restraint
619 16466699 : IF (PRESENT(ng4x6_restraint)) ng4x6_restraint = molecule_kind%ng4x6_restraint
620 16466699 : IF (PRESENT(nvsite_restraint)) nvsite_restraint = molecule_kind%nvsite_restraint
621 16466699 : IF (PRESENT(nfixd_restraint)) nfixd_restraint = molecule_kind%nfixd_restraint
622 16466699 : IF (PRESENT(nrestraints)) nrestraints = molecule_kind%ncolv%nrestraint + &
623 : molecule_kind%ng3x3_restraint + &
624 : molecule_kind%ng4x6_restraint + &
625 264760 : molecule_kind%nvsite_restraint
626 16466699 : IF (PRESENT(nmolecule)) nmolecule = molecule_kind%nmolecule
627 16466699 : IF (PRESENT(nshell)) nshell = molecule_kind%nshell
628 16466699 : IF (PRESENT(ntorsion)) ntorsion = molecule_kind%ntorsion
629 16466699 : IF (PRESENT(nsgf)) nsgf = molecule_kind%nsgf
630 16466699 : IF (PRESENT(nelectron)) nelectron = molecule_kind%nelectron
631 16466699 : IF (PRESENT(nelectron_alpha)) nelectron_alpha = molecule_kind%nelectron_beta
632 16466699 : IF (PRESENT(nelectron_beta)) nelectron_beta = molecule_kind%nelectron_alpha
633 16466699 : IF (PRESENT(molecule_list)) molecule_list => molecule_kind%molecule_list
634 :
635 16466699 : END SUBROUTINE get_molecule_kind
636 :
637 : ! **************************************************************************************************
638 : !> \brief Get informations about a molecule kind set.
639 : !> \param molecule_kind_set ...
640 : !> \param maxatom ...
641 : !> \param natom ...
642 : !> \param nbond ...
643 : !> \param nbend ...
644 : !> \param nub ...
645 : !> \param ntorsion ...
646 : !> \param nimpr ...
647 : !> \param nopbend ...
648 : !> \param nconstraint ...
649 : !> \param nconstraint_fixd ...
650 : !> \param nmolecule ...
651 : !> \param nrestraints ...
652 : !> \date 27.08.2003
653 : !> \author Matthias Krack
654 : !> \version 1.0
655 : ! **************************************************************************************************
656 49036 : SUBROUTINE get_molecule_kind_set(molecule_kind_set, maxatom, natom, &
657 : nbond, nbend, nub, ntorsion, nimpr, nopbend, &
658 : nconstraint, nconstraint_fixd, nmolecule, &
659 : nrestraints)
660 :
661 : TYPE(molecule_kind_type), DIMENSION(:), INTENT(IN) :: molecule_kind_set
662 : INTEGER, INTENT(OUT), OPTIONAL :: maxatom, natom, nbond, nbend, nub, &
663 : ntorsion, nimpr, nopbend, nconstraint, &
664 : nconstraint_fixd, nmolecule, &
665 : nrestraints
666 :
667 : INTEGER :: ibend, ibond, iimpr, imolecule_kind, iopbend, itorsion, iub, na, nc, nc_fixd, &
668 : nfixd_restraint, nm, nmolecule_kind, nrestraints_tot
669 :
670 49036 : IF (PRESENT(maxatom)) maxatom = 0
671 49036 : IF (PRESENT(natom)) natom = 0
672 49036 : IF (PRESENT(nbond)) nbond = 0
673 49036 : IF (PRESENT(nbend)) nbend = 0
674 49036 : IF (PRESENT(nub)) nub = 0
675 49036 : IF (PRESENT(ntorsion)) ntorsion = 0
676 49036 : IF (PRESENT(nimpr)) nimpr = 0
677 49036 : IF (PRESENT(nopbend)) nopbend = 0
678 49036 : IF (PRESENT(nconstraint)) nconstraint = 0
679 49036 : IF (PRESENT(nconstraint_fixd)) nconstraint_fixd = 0
680 49036 : IF (PRESENT(nrestraints)) nrestraints = 0
681 49036 : IF (PRESENT(nmolecule)) nmolecule = 0
682 :
683 49036 : nmolecule_kind = SIZE(molecule_kind_set)
684 :
685 313796 : DO imolecule_kind = 1, nmolecule_kind
686 49036 : ASSOCIATE (molecule_kind => molecule_kind_set(imolecule_kind))
687 :
688 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
689 : natom=na, &
690 : nbond=ibond, &
691 : nbend=ibend, &
692 : nub=iub, &
693 : ntorsion=itorsion, &
694 : nimpr=iimpr, &
695 : nopbend=iopbend, &
696 : nconstraint=nc, &
697 : nconstraint_fixd=nc_fixd, &
698 : nfixd_restraint=nfixd_restraint, &
699 : nrestraints=nrestraints_tot, &
700 264760 : nmolecule=nm)
701 264760 : IF (PRESENT(maxatom)) maxatom = MAX(maxatom, na)
702 264760 : IF (PRESENT(natom)) natom = natom + na*nm
703 264760 : IF (PRESENT(nbond)) nbond = nbond + ibond*nm
704 264760 : IF (PRESENT(nbend)) nbend = nbend + ibend*nm
705 264760 : IF (PRESENT(nub)) nub = nub + iub*nm
706 264760 : IF (PRESENT(ntorsion)) ntorsion = ntorsion + itorsion*nm
707 264760 : IF (PRESENT(nimpr)) nimpr = nimpr + iimpr*nm
708 264760 : IF (PRESENT(nopbend)) nopbend = nopbend + iopbend*nm
709 264760 : IF (PRESENT(nconstraint)) nconstraint = nconstraint + nc*nm + nc_fixd
710 264760 : IF (PRESENT(nconstraint_fixd)) nconstraint_fixd = nconstraint_fixd + nc_fixd
711 264760 : IF (PRESENT(nmolecule)) nmolecule = nmolecule + nm
712 529520 : IF (PRESENT(nrestraints)) nrestraints = nrestraints + nm*nrestraints_tot + nfixd_restraint
713 :
714 : END ASSOCIATE
715 : END DO
716 :
717 49036 : END SUBROUTINE get_molecule_kind_set
718 :
719 : ! **************************************************************************************************
720 : !> \brief Set the components of a molecule kind.
721 : !> \param molecule_kind ...
722 : !> \param name ...
723 : !> \param mass ...
724 : !> \param charge ...
725 : !> \param kind_number ...
726 : !> \param molecule_list ...
727 : !> \param atom_list ...
728 : !> \param nbond ...
729 : !> \param bond_list ...
730 : !> \param nbend ...
731 : !> \param bend_list ...
732 : !> \param nub ...
733 : !> \param ub_list ...
734 : !> \param nimpr ...
735 : !> \param impr_list ...
736 : !> \param nopbend ...
737 : !> \param opbend_list ...
738 : !> \param ntorsion ...
739 : !> \param torsion_list ...
740 : !> \param fixd_list ...
741 : !> \param ncolv ...
742 : !> \param colv_list ...
743 : !> \param ng3x3 ...
744 : !> \param g3x3_list ...
745 : !> \param ng4x6 ...
746 : !> \param nfixd ...
747 : !> \param g4x6_list ...
748 : !> \param nvsite ...
749 : !> \param vsite_list ...
750 : !> \param ng3x3_restraint ...
751 : !> \param ng4x6_restraint ...
752 : !> \param nfixd_restraint ...
753 : !> \param nshell ...
754 : !> \param shell_list ...
755 : !> \param nvsite_restraint ...
756 : !> \param bond_kind_set ...
757 : !> \param bend_kind_set ...
758 : !> \param ub_kind_set ...
759 : !> \param torsion_kind_set ...
760 : !> \param impr_kind_set ...
761 : !> \param opbend_kind_set ...
762 : !> \param nelectron ...
763 : !> \param nsgf ...
764 : !> \param molname_generated ...
765 : !> \date 27.08.2003
766 : !> \author Matthias Krack
767 : !> \version 1.0
768 : ! **************************************************************************************************
769 2055811 : SUBROUTINE set_molecule_kind(molecule_kind, name, mass, charge, kind_number, &
770 : molecule_list, atom_list, nbond, bond_list, &
771 : nbend, bend_list, nub, ub_list, nimpr, impr_list, &
772 : nopbend, opbend_list, ntorsion, &
773 : torsion_list, fixd_list, ncolv, colv_list, ng3x3, &
774 : g3x3_list, ng4x6, nfixd, g4x6_list, nvsite, &
775 : vsite_list, ng3x3_restraint, ng4x6_restraint, &
776 : nfixd_restraint, nshell, shell_list, &
777 : nvsite_restraint, bond_kind_set, bend_kind_set, &
778 : ub_kind_set, torsion_kind_set, impr_kind_set, &
779 : opbend_kind_set, nelectron, nsgf, &
780 : molname_generated)
781 :
782 : TYPE(molecule_kind_type), INTENT(INOUT) :: molecule_kind
783 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name
784 : REAL(KIND=dp), OPTIONAL :: mass, charge
785 : INTEGER, INTENT(IN), OPTIONAL :: kind_number
786 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: molecule_list
787 : TYPE(atom_type), DIMENSION(:), OPTIONAL, POINTER :: atom_list
788 : INTEGER, INTENT(IN), OPTIONAL :: nbond
789 : TYPE(bond_type), DIMENSION(:), OPTIONAL, POINTER :: bond_list
790 : INTEGER, INTENT(IN), OPTIONAL :: nbend
791 : TYPE(bend_type), DIMENSION(:), OPTIONAL, POINTER :: bend_list
792 : INTEGER, INTENT(IN), OPTIONAL :: nub
793 : TYPE(ub_type), DIMENSION(:), OPTIONAL, POINTER :: ub_list
794 : INTEGER, INTENT(IN), OPTIONAL :: nimpr
795 : TYPE(impr_type), DIMENSION(:), OPTIONAL, POINTER :: impr_list
796 : INTEGER, INTENT(IN), OPTIONAL :: nopbend
797 : TYPE(opbend_type), DIMENSION(:), OPTIONAL, POINTER :: opbend_list
798 : INTEGER, INTENT(IN), OPTIONAL :: ntorsion
799 : TYPE(torsion_type), DIMENSION(:), OPTIONAL, &
800 : POINTER :: torsion_list
801 : TYPE(fixd_constraint_type), DIMENSION(:), &
802 : OPTIONAL, POINTER :: fixd_list
803 : TYPE(colvar_counters), INTENT(IN), OPTIONAL :: ncolv
804 : TYPE(colvar_constraint_type), DIMENSION(:), &
805 : OPTIONAL, POINTER :: colv_list
806 : INTEGER, INTENT(IN), OPTIONAL :: ng3x3
807 : TYPE(g3x3_constraint_type), DIMENSION(:), &
808 : OPTIONAL, POINTER :: g3x3_list
809 : INTEGER, INTENT(IN), OPTIONAL :: ng4x6, nfixd
810 : TYPE(g4x6_constraint_type), DIMENSION(:), &
811 : OPTIONAL, POINTER :: g4x6_list
812 : INTEGER, INTENT(IN), OPTIONAL :: nvsite
813 : TYPE(vsite_constraint_type), DIMENSION(:), &
814 : OPTIONAL, POINTER :: vsite_list
815 : INTEGER, INTENT(IN), OPTIONAL :: ng3x3_restraint, ng4x6_restraint, &
816 : nfixd_restraint, nshell
817 : TYPE(shell_type), DIMENSION(:), OPTIONAL, POINTER :: shell_list
818 : INTEGER, INTENT(IN), OPTIONAL :: nvsite_restraint
819 : TYPE(bond_kind_type), DIMENSION(:), OPTIONAL, &
820 : POINTER :: bond_kind_set
821 : TYPE(bend_kind_type), DIMENSION(:), OPTIONAL, &
822 : POINTER :: bend_kind_set
823 : TYPE(ub_kind_type), DIMENSION(:), OPTIONAL, &
824 : POINTER :: ub_kind_set
825 : TYPE(torsion_kind_type), DIMENSION(:), OPTIONAL, &
826 : POINTER :: torsion_kind_set
827 : TYPE(impr_kind_type), DIMENSION(:), OPTIONAL, &
828 : POINTER :: impr_kind_set
829 : TYPE(opbend_kind_type), DIMENSION(:), OPTIONAL, &
830 : POINTER :: opbend_kind_set
831 : INTEGER, INTENT(IN), OPTIONAL :: nelectron, nsgf
832 : LOGICAL, INTENT(IN), OPTIONAL :: molname_generated
833 :
834 : INTEGER :: n
835 :
836 2055811 : IF (PRESENT(atom_list)) THEN
837 272866 : n = SIZE(atom_list)
838 272866 : molecule_kind%natom = n
839 272866 : molecule_kind%atom_list => atom_list
840 : END IF
841 2055811 : IF (PRESENT(molname_generated)) molecule_kind%molname_generated = molname_generated
842 2055811 : IF (PRESENT(name)) molecule_kind%name = name
843 2055811 : IF (PRESENT(mass)) molecule_kind%mass = mass
844 2055811 : IF (PRESENT(charge)) molecule_kind%charge = charge
845 2055811 : IF (PRESENT(kind_number)) molecule_kind%kind_number = kind_number
846 2055811 : IF (PRESENT(nbond)) molecule_kind%nbond = nbond
847 2055811 : IF (PRESENT(bond_list)) molecule_kind%bond_list => bond_list
848 2055811 : IF (PRESENT(nbend)) molecule_kind%nbend = nbend
849 2055811 : IF (PRESENT(nelectron)) molecule_kind%nelectron = nelectron
850 2055811 : IF (PRESENT(nsgf)) molecule_kind%nsgf = nsgf
851 2055811 : IF (PRESENT(bend_list)) molecule_kind%bend_list => bend_list
852 2055811 : IF (PRESENT(nub)) molecule_kind%nub = nub
853 2055811 : IF (PRESENT(ub_list)) molecule_kind%ub_list => ub_list
854 2055811 : IF (PRESENT(ntorsion)) molecule_kind%ntorsion = ntorsion
855 2055811 : IF (PRESENT(torsion_list)) molecule_kind%torsion_list => torsion_list
856 2055811 : IF (PRESENT(nimpr)) molecule_kind%nimpr = nimpr
857 2055811 : IF (PRESENT(impr_list)) molecule_kind%impr_list => impr_list
858 2055811 : IF (PRESENT(nopbend)) molecule_kind%nopbend = nopbend
859 2055811 : IF (PRESENT(opbend_list)) molecule_kind%opbend_list => opbend_list
860 2055811 : IF (PRESENT(ncolv)) molecule_kind%ncolv = ncolv
861 2055811 : IF (PRESENT(colv_list)) molecule_kind%colv_list => colv_list
862 2055811 : IF (PRESENT(ng3x3)) molecule_kind%ng3x3 = ng3x3
863 2055811 : IF (PRESENT(g3x3_list)) molecule_kind%g3x3_list => g3x3_list
864 2055811 : IF (PRESENT(ng4x6)) molecule_kind%ng4x6 = ng4x6
865 2055811 : IF (PRESENT(nvsite)) molecule_kind%nvsite = nvsite
866 2055811 : IF (PRESENT(nfixd)) molecule_kind%nfixd = nfixd
867 2055811 : IF (PRESENT(nfixd_restraint)) molecule_kind%nfixd_restraint = nfixd_restraint
868 2055811 : IF (PRESENT(ng3x3_restraint)) molecule_kind%ng3x3_restraint = ng3x3_restraint
869 2055811 : IF (PRESENT(ng4x6_restraint)) molecule_kind%ng4x6_restraint = ng4x6_restraint
870 2055811 : IF (PRESENT(nvsite_restraint)) molecule_kind%nvsite_restraint = nvsite_restraint
871 2055811 : IF (PRESENT(g4x6_list)) molecule_kind%g4x6_list => g4x6_list
872 2055811 : IF (PRESENT(vsite_list)) molecule_kind%vsite_list => vsite_list
873 2055811 : IF (PRESENT(fixd_list)) molecule_kind%fixd_list => fixd_list
874 2055811 : IF (PRESENT(bond_kind_set)) molecule_kind%bond_kind_set => bond_kind_set
875 2055811 : IF (PRESENT(bend_kind_set)) molecule_kind%bend_kind_set => bend_kind_set
876 2055811 : IF (PRESENT(ub_kind_set)) molecule_kind%ub_kind_set => ub_kind_set
877 2055811 : IF (PRESENT(torsion_kind_set)) molecule_kind%torsion_kind_set => torsion_kind_set
878 2055811 : IF (PRESENT(impr_kind_set)) molecule_kind%impr_kind_set => impr_kind_set
879 2055811 : IF (PRESENT(opbend_kind_set)) molecule_kind%opbend_kind_set => opbend_kind_set
880 2055811 : IF (PRESENT(nshell)) molecule_kind%nshell = nshell
881 2055811 : IF (PRESENT(shell_list)) molecule_kind%shell_list => shell_list
882 2055811 : IF (PRESENT(molecule_list)) THEN
883 136433 : n = SIZE(molecule_list)
884 136433 : molecule_kind%nmolecule = n
885 136433 : molecule_kind%molecule_list => molecule_list
886 : END IF
887 2055811 : END SUBROUTINE set_molecule_kind
888 :
889 : ! **************************************************************************************************
890 : !> \brief Write a molecule kind data set to the output unit.
891 : !> \param molecule_kind ...
892 : !> \param output_unit ...
893 : !> \date 24.09.2003
894 : !> \author Matthias Krack
895 : !> \version 1.0
896 : ! **************************************************************************************************
897 2232 : SUBROUTINE write_molecule_kind(molecule_kind, output_unit)
898 : TYPE(molecule_kind_type), INTENT(IN) :: molecule_kind
899 : INTEGER, INTENT(in) :: output_unit
900 :
901 : CHARACTER(LEN=default_string_length) :: name
902 : INTEGER :: iatom, imolecule, natom, nmolecule
903 : TYPE(atomic_kind_type), POINTER :: atomic_kind
904 :
905 2232 : IF (output_unit > 0) THEN
906 2232 : natom = SIZE(molecule_kind%atom_list)
907 2232 : nmolecule = SIZE(molecule_kind%molecule_list)
908 :
909 2232 : IF (natom == 1) THEN
910 211 : atomic_kind => molecule_kind%atom_list(1)%atomic_kind
911 211 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
912 : WRITE (UNIT=output_unit, FMT="(/,T2,I5,A,T36,A,A,T64,A)") &
913 211 : molecule_kind%kind_number, &
914 211 : ". Molecule kind: "//TRIM(molecule_kind%name), &
915 422 : "Atomic kind name: ", TRIM(name)
916 : WRITE (UNIT=output_unit, FMT="(T9,A,L1,T55,A,T75,I6)") &
917 211 : "Automatic name: ", molecule_kind%molname_generated, &
918 422 : "Number of molecules:", nmolecule
919 : ELSE
920 : WRITE (UNIT=output_unit, FMT="(/,T2,I5,A,T50,A,T75,I6,/,T22,A)") &
921 2021 : molecule_kind%kind_number, &
922 2021 : ". Molecule kind: "//TRIM(molecule_kind%name), &
923 2021 : "Number of atoms: ", natom, &
924 4042 : "Atom Atomic kind name"
925 17014 : DO iatom = 1, natom
926 14993 : atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind
927 14993 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
928 : WRITE (UNIT=output_unit, FMT="(T20,I6,(7X,A18))") &
929 17014 : iatom, TRIM(name)
930 : END DO
931 : WRITE (UNIT=output_unit, FMT="(/,T9,A,L1)") &
932 2021 : "The name was automatically generated: ", &
933 4042 : molecule_kind%molname_generated
934 : WRITE (UNIT=output_unit, FMT="(T9,A,I6,/,T9,A,(T30,5I10))") &
935 2021 : "Number of molecules: ", nmolecule, "Molecule list:", &
936 33831 : (molecule_kind%molecule_list(imolecule), imolecule=1, nmolecule)
937 2021 : IF (molecule_kind%nbond > 0) &
938 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
939 1745 : "Number of bonds: ", molecule_kind%nbond
940 2021 : IF (molecule_kind%nbend > 0) &
941 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
942 1619 : "Number of bends: ", molecule_kind%nbend
943 2021 : IF (molecule_kind%nub > 0) &
944 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
945 266 : "Number of Urey-Bradley:", molecule_kind%nub
946 2021 : IF (molecule_kind%ntorsion > 0) &
947 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
948 1122 : "Number of torsions: ", molecule_kind%ntorsion
949 2021 : IF (molecule_kind%nimpr > 0) &
950 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
951 179 : "Number of improper: ", molecule_kind%nimpr
952 2021 : IF (molecule_kind%nopbend > 0) &
953 : WRITE (UNIT=output_unit, FMT="(1X,A30,I6)") &
954 4 : "Number of out opbends: ", molecule_kind%nopbend
955 : END IF
956 : END IF
957 2232 : END SUBROUTINE write_molecule_kind
958 :
959 : ! **************************************************************************************************
960 : !> \brief Write a moleculeatomic kind set data set to the output unit.
961 : !> \param molecule_kind_set ...
962 : !> \param subsys_section ...
963 : !> \date 24.09.2003
964 : !> \author Matthias Krack
965 : !> \version 1.0
966 : ! **************************************************************************************************
967 10255 : SUBROUTINE write_molecule_kind_set(molecule_kind_set, subsys_section)
968 : TYPE(molecule_kind_type), DIMENSION(:), INTENT(IN) :: molecule_kind_set
969 : TYPE(section_vals_type), INTENT(IN) :: subsys_section
970 :
971 : CHARACTER(len=*), PARAMETER :: routineN = 'write_molecule_kind_set'
972 :
973 : INTEGER :: handle, imolecule_kind, natom, nbend, &
974 : nbond, nimpr, nmolecule, &
975 : nmolecule_kind, nopbend, ntors, &
976 : ntotal, nub, output_unit
977 : LOGICAL :: all_single_atoms
978 : TYPE(cp_logger_type), POINTER :: logger
979 :
980 10255 : CALL timeset(routineN, handle)
981 :
982 10255 : NULLIFY (logger)
983 10255 : logger => cp_get_default_logger()
984 : output_unit = cp_print_key_unit_nr(logger, subsys_section, &
985 10255 : "PRINT%MOLECULES", extension=".Log")
986 10255 : IF (output_unit > 0) THEN
987 2591 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "MOLECULE KIND INFORMATION"
988 :
989 2591 : nmolecule_kind = SIZE(molecule_kind_set)
990 :
991 2591 : all_single_atoms = .TRUE.
992 31918 : DO imolecule_kind = 1, nmolecule_kind
993 29327 : natom = SIZE(molecule_kind_set(imolecule_kind)%atom_list)
994 29327 : nmolecule = SIZE(molecule_kind_set(imolecule_kind)%molecule_list)
995 31918 : IF (natom*nmolecule > 1) all_single_atoms = .FALSE.
996 : END DO
997 :
998 2591 : IF (all_single_atoms) THEN
999 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
1000 1901 : "All atoms are their own molecule, skipping detailed information"
1001 : ELSE
1002 2922 : DO imolecule_kind = 1, nmolecule_kind
1003 2922 : CALL write_molecule_kind(molecule_kind_set(imolecule_kind), output_unit)
1004 : END DO
1005 : END IF
1006 :
1007 : CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
1008 : nbond=nbond, &
1009 : nbend=nbend, &
1010 : nub=nub, &
1011 : ntorsion=ntors, &
1012 : nimpr=nimpr, &
1013 2591 : nopbend=nopbend)
1014 2591 : ntotal = nbond + nbend + nub + ntors + nimpr + nopbend
1015 2591 : IF (ntotal > 0) THEN
1016 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A,T45,A30,I6)") &
1017 588 : "MOLECULE KIND SET INFORMATION", &
1018 1176 : "Total Number of bonds: ", nbond
1019 : WRITE (UNIT=output_unit, FMT="(T45,A30,I6)") &
1020 588 : "Total Number of bends: ", nbend
1021 : WRITE (UNIT=output_unit, FMT="(T45,A30,I6)") &
1022 588 : "Total Number of Urey-Bradley:", nub
1023 : WRITE (UNIT=output_unit, FMT="(T45,A30,I6)") &
1024 588 : "Total Number of torsions: ", ntors
1025 : WRITE (UNIT=output_unit, FMT="(T45,A30,I6)") &
1026 588 : "Total Number of improper: ", nimpr
1027 : WRITE (UNIT=output_unit, FMT="(T45,A30,I6)") &
1028 588 : "Total Number of opbends: ", nopbend
1029 : END IF
1030 : END IF
1031 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
1032 10255 : "PRINT%MOLECULES")
1033 :
1034 10255 : CALL timestop(handle)
1035 :
1036 10255 : END SUBROUTINE write_molecule_kind_set
1037 :
1038 : ! **************************************************************************************************
1039 : !> \brief Write collective variable constraint information to output unit
1040 : !> \param colvar_constraint Data set of the collective variable constraint
1041 : !> \param icolv Collective variable number (index)
1042 : !> \param iw Logical unit number of the output unit
1043 : !> \author Matthias Krack (25.11.2025)
1044 : ! **************************************************************************************************
1045 26 : SUBROUTINE write_colvar_constraint(colvar_constraint, icolv, iw)
1046 :
1047 : TYPE(colvar_constraint_type), INTENT(IN), POINTER :: colvar_constraint
1048 : INTEGER, INTENT(IN) :: icolv, iw
1049 :
1050 : CHARACTER(LEN=30) :: type_string
1051 :
1052 26 : IF (iw > 0) THEN
1053 26 : CPASSERT(ASSOCIATED(colvar_constraint))
1054 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1055 26 : "COLVAR| Number", icolv
1056 26 : SELECT CASE (colvar_constraint%type_id)
1057 : CASE (no_colvar_id)
1058 0 : type_string = "Undefined"
1059 : CASE (dist_colvar_id)
1060 19 : type_string = "Distance"
1061 : CASE (coord_colvar_id)
1062 1 : type_string = "Coordination number"
1063 : CASE (torsion_colvar_id)
1064 0 : type_string = "Torsion"
1065 : CASE (angle_colvar_id)
1066 2 : type_string = "Angle"
1067 : CASE (plane_distance_colvar_id)
1068 0 : type_string = "Plane distance"
1069 : CASE (rotation_colvar_id)
1070 0 : type_string = "Rotation"
1071 : CASE (dfunct_colvar_id)
1072 2 : type_string = "Distance function"
1073 : CASE (qparm_colvar_id)
1074 0 : type_string = "Q parameter"
1075 : CASE (hydronium_shell_colvar_id)
1076 0 : type_string = "Hydronium shell"
1077 : CASE (reaction_path_colvar_id)
1078 0 : type_string = "Reaction path"
1079 : CASE (combine_colvar_id)
1080 0 : type_string = "Combine"
1081 : CASE (population_colvar_id)
1082 0 : type_string = "Population"
1083 : CASE (plane_plane_angle_colvar_id)
1084 2 : type_string = "Angle plane-plane"
1085 : CASE (gyration_colvar_id)
1086 0 : type_string = "Gyration radius"
1087 : CASE (rmsd_colvar_id)
1088 0 : type_string = "RMSD"
1089 : CASE (distance_from_path_colvar_id)
1090 0 : type_string = "Distance from path"
1091 : CASE (xyz_diag_colvar_id)
1092 0 : type_string = "XYZ diag"
1093 : CASE (xyz_outerdiag_colvar_id)
1094 0 : type_string = "XYZ outerdiag"
1095 : CASE (u_colvar_id)
1096 0 : type_string = "U"
1097 : CASE (Wc_colvar_id)
1098 0 : type_string = "WC"
1099 : CASE (HBP_colvar_id)
1100 0 : type_string = "HBP"
1101 : CASE (ring_puckering_colvar_id)
1102 0 : type_string = "Ring puckering"
1103 : CASE (mindist_colvar_id)
1104 0 : type_string = "Distance point-plane"
1105 : CASE (acid_hyd_dist_colvar_id)
1106 0 : type_string = "Acid hydronium distance"
1107 : CASE (acid_hyd_shell_colvar_id)
1108 0 : type_string = "Acid hydronium shell"
1109 : CASE (hydronium_dist_colvar_id)
1110 0 : type_string = "Hydronium distance"
1111 : CASE DEFAULT
1112 26 : CPABORT("Invalid collective variable ID specified. Check the code!")
1113 : END SELECT
1114 26 : IF (colvar_constraint%restraint%active) THEN
1115 : WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
1116 7 : "COLVAR| Restraint type", ADJUSTR(TRIM(type_string))
1117 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1118 7 : "COLVAR| Restraint constant k [a.u.]", colvar_constraint%restraint%k0
1119 : ELSE
1120 : WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
1121 19 : "COLVAR| Constraint type", ADJUSTR(TRIM(type_string))
1122 : END IF
1123 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1124 26 : "COLVAR| Target value", colvar_constraint%expected_value, &
1125 52 : "COLVAR| Target value growth speed", colvar_constraint%expected_value_growth_speed
1126 26 : IF (colvar_constraint%use_points) THEN
1127 4 : WRITE (UNIT=iw, FMT="(T2,A,T78,A3)") "COLVAR| Use points", "Yes"
1128 : ELSE
1129 22 : WRITE (UNIT=iw, FMT="(T2,A,T79,A2)") "COLVAR| Use points", "No"
1130 : END IF
1131 : END IF
1132 :
1133 26 : END SUBROUTINE write_colvar_constraint
1134 :
1135 : ! **************************************************************************************************
1136 : !> \brief Write fix atom constraint information to output unit
1137 : !> \param fixd_constraint Data set of the fix atom constraint
1138 : !> \param ifixd Fix atom constraint/restraint number (index)
1139 : !> \param iw Logical unit number of the output unit
1140 : !> \author Matthias Krack (26.11.2025)
1141 : ! **************************************************************************************************
1142 2 : SUBROUTINE write_fixd_constraint(fixd_constraint, ifixd, iw)
1143 :
1144 : TYPE(fixd_constraint_type), INTENT(IN), POINTER :: fixd_constraint
1145 : INTEGER, INTENT(IN) :: ifixd, iw
1146 :
1147 2 : IF (iw > 0) THEN
1148 2 : CPASSERT(ASSOCIATED(fixd_constraint))
1149 2 : IF (fixd_constraint%restraint%active) THEN
1150 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1151 2 : "FIX_ATOM| Number (restraint)", ifixd
1152 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1153 2 : "FIX_ATOM| Restraint constant k [a.u.]", fixd_constraint%restraint%k0
1154 : ELSE
1155 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1156 0 : "FIX_ATOM| Number (constraint)", ifixd
1157 : END IF
1158 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
1159 2 : "FIX_ATOM| Atom index", fixd_constraint%fixd
1160 : WRITE (UNIT=iw, FMT="(T2,A,T78,A3)") &
1161 2 : "FIX_ATOM| Fixed Cartesian components", periodicity_string(fixd_constraint%itype)
1162 2 : IF (INDEX(periodicity_string(fixd_constraint%itype), "X") > 0) THEN
1163 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1164 2 : "FIX_ATOM| X coordinate [Angstrom]", cp_unit_from_cp2k(fixd_constraint%coord(1), "Angstrom")
1165 : END IF
1166 2 : IF (INDEX(periodicity_string(fixd_constraint%itype), "Y") > 0) THEN
1167 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1168 2 : "FIX_ATOM| Y coordinate [Angstrom]", cp_unit_from_cp2k(fixd_constraint%coord(2), "Angstrom")
1169 : END IF
1170 2 : IF (INDEX(periodicity_string(fixd_constraint%itype), "Z") > 0) THEN
1171 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1172 2 : "FIX_ATOM| Z coordinate [Angstrom]", cp_unit_from_cp2k(fixd_constraint%coord(3), "Angstrom")
1173 : END IF
1174 : END IF
1175 :
1176 2 : END SUBROUTINE write_fixd_constraint
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief Write G3x3 constraint information to output unit
1180 : !> \param g3x3_constraint Data set of the g3x3 constraint
1181 : !> \param ig3x3 G3x3 constraint/restraint number (index)
1182 : !> \param iw Logical unit number of the output unit
1183 : !> \author Matthias Krack (26.11.2025)
1184 : ! **************************************************************************************************
1185 2 : SUBROUTINE write_g3x3_constraint(g3x3_constraint, ig3x3, iw)
1186 :
1187 : TYPE(g3x3_constraint_type), INTENT(IN), POINTER :: g3x3_constraint
1188 : INTEGER, INTENT(IN) :: ig3x3, iw
1189 :
1190 2 : IF (iw > 0) THEN
1191 2 : CPASSERT(ASSOCIATED(g3x3_constraint))
1192 2 : IF (g3x3_constraint%restraint%active) THEN
1193 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1194 1 : "G3X3| Number (restraint)", ig3x3
1195 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1196 1 : "G3X3| Restraint constant k [a.u.]", g3x3_constraint%restraint%k0
1197 : ELSE
1198 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1199 1 : "G3X3| Number (constraint)", ig3x3
1200 : END IF
1201 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
1202 2 : "G3X3| Atom index a", g3x3_constraint%a, &
1203 2 : "G3X3| Atom index b", g3x3_constraint%b, &
1204 4 : "G3X3| Atom index c", g3x3_constraint%c
1205 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1206 2 : "G3X3| Distance (a,b) [Angstrom]", cp_unit_from_cp2k(g3x3_constraint%dab, "Angstrom"), &
1207 2 : "G3X3| Distance (a,c) [Angstrom]", cp_unit_from_cp2k(g3x3_constraint%dac, "Angstrom"), &
1208 4 : "G3X3| Distance (b,c) [Angstrom]", cp_unit_from_cp2k(g3x3_constraint%dbc, "Angstrom")
1209 : END IF
1210 :
1211 2 : END SUBROUTINE write_g3x3_constraint
1212 :
1213 : ! **************************************************************************************************
1214 : !> \brief Write G4x6 constraint information to output unit
1215 : !> \param g4x6_constraint Data set of the g4x6 constraint
1216 : !> \param ig4x6 G4x6 constraint/restraint number (index)
1217 : !> \param iw Logical unit number of the output unit
1218 : !> \author Matthias Krack (26.11.2025)
1219 : ! **************************************************************************************************
1220 2 : SUBROUTINE write_g4x6_constraint(g4x6_constraint, ig4x6, iw)
1221 :
1222 : TYPE(g4x6_constraint_type), INTENT(IN), POINTER :: g4x6_constraint
1223 : INTEGER, INTENT(IN) :: ig4x6, iw
1224 :
1225 2 : IF (iw > 0) THEN
1226 2 : CPASSERT(ASSOCIATED(g4x6_constraint))
1227 2 : IF (g4x6_constraint%restraint%active) THEN
1228 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1229 1 : "G4X6| Number (restraint)", ig4x6
1230 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1231 1 : "G4X6| Restraint constant k [a.u.]", g4x6_constraint%restraint%k0
1232 : ELSE
1233 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1234 1 : "G4X6| Number (constraint)", ig4x6
1235 : END IF
1236 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
1237 2 : "G4X6| Atom index a", g4x6_constraint%a, &
1238 2 : "G4X6| Atom index b", g4x6_constraint%b, &
1239 2 : "G4X6| Atom index c", g4x6_constraint%c, &
1240 4 : "G4X6| Atom index d", g4x6_constraint%d
1241 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1242 2 : "G4X6| Distance (a,b) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dab, "Angstrom"), &
1243 2 : "G4X6| Distance (a,c) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dac, "Angstrom"), &
1244 2 : "G4X6| Distance (a,d) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dad, "Angstrom"), &
1245 2 : "G4X6| Distance (b,c) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dbc, "Angstrom"), &
1246 2 : "G4X6| Distance (b,d) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dbd, "Angstrom"), &
1247 4 : "G4X6| Distance (c,d) [Angstrom]", cp_unit_from_cp2k(g4x6_constraint%dcd, "Angstrom")
1248 : END IF
1249 :
1250 2 : END SUBROUTINE write_g4x6_constraint
1251 :
1252 : ! **************************************************************************************************
1253 : !> \brief Write virtual site constraint information to output unit
1254 : !> \param vsite_constraint Data set of the vsite constraint
1255 : !> \param ivsite Virtual site constraint/restraint number (index)
1256 : !> \param iw Logical unit number of the output unit
1257 : !> \author Matthias Krack (01.12.2025)
1258 : ! **************************************************************************************************
1259 0 : SUBROUTINE write_vsite_constraint(vsite_constraint, ivsite, iw)
1260 :
1261 : TYPE(vsite_constraint_type), INTENT(IN), POINTER :: vsite_constraint
1262 : INTEGER, INTENT(IN) :: ivsite, iw
1263 :
1264 0 : IF (iw > 0) THEN
1265 0 : CPASSERT(ASSOCIATED(vsite_constraint))
1266 0 : IF (vsite_constraint%restraint%active) THEN
1267 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1268 0 : "VSITE| Number (restraint)", ivsite
1269 : WRITE (UNIT=iw, FMT="(T2,A,T66,ES15.6)") &
1270 0 : "VSITE| Restraint constant k [a.u.]", vsite_constraint%restraint%k0
1271 : ELSE
1272 : WRITE (UNIT=iw, FMT="(/,T2,A,T71,I10)") &
1273 0 : "VSITE| Number (constraint)", ivsite
1274 : END IF
1275 : WRITE (UNIT=iw, FMT="(T2,A,T71,I10)") &
1276 0 : "VSITE| Atom index of virtual site", vsite_constraint%a, &
1277 0 : "VSITE| Atom index b", vsite_constraint%b, &
1278 0 : "VSITE| Atom index c", vsite_constraint%c, &
1279 0 : "VSITE| Atom index d", vsite_constraint%d
1280 : WRITE (UNIT=iw, FMT="(T2,A,T66,F15.8)") &
1281 0 : "VSITE| Distance (b,c) [Angstrom]", cp_unit_from_cp2k(vsite_constraint%wbc, "Angstrom"), &
1282 0 : "VSITE| Distance (d,c) [Angstrom]", cp_unit_from_cp2k(vsite_constraint%wdc, "Angstrom")
1283 : END IF
1284 :
1285 0 : END SUBROUTINE write_vsite_constraint
1286 :
1287 0 : END MODULE molecule_kind_types
|
