Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Reads the input sections "topology"
10 : !> \par History
11 : !> JGH (26-01-2002) Added read_topology_section
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE topology_input
15 : USE colvar_types, ONLY: colvar_clone,&
16 : colvar_p_type
17 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
18 : cp_to_string
19 : USE input_constants, ONLY: do_conn_generate,&
20 : do_conn_mol_set,&
21 : do_conn_off,&
22 : do_conn_user,&
23 : do_constr_none,&
24 : do_coord_off
25 : USE input_section_types, ONLY: section_vals_get,&
26 : section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get,&
29 : section_vals_val_unset
30 : USE kinds, ONLY: default_string_length,&
31 : dp
32 : USE memory_utilities, ONLY: reallocate
33 : USE topology_types, ONLY: constraint_info_type,&
34 : topology_parameters_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_input'
40 :
41 : PRIVATE
42 : PUBLIC :: read_topology_section, read_constraints_section
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief reads the input section topology
48 : !> \param topology ...
49 : !> \param topology_section ...
50 : !> \par History
51 : !> none
52 : !> \author JGH (26-01-2002)
53 : ! **************************************************************************************************
54 10950 : SUBROUTINE read_topology_section(topology, topology_section)
55 : TYPE(topology_parameters_type) :: topology
56 : TYPE(section_vals_type), POINTER :: topology_section
57 :
58 : CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_section'
59 :
60 : INTEGER :: handle, ival
61 :
62 10950 : CALL timeset(routineN, handle)
63 10950 : CALL section_vals_val_get(topology_section, "CHARGE_OCCUP", l_val=topology%charge_occup)
64 10950 : CALL section_vals_val_get(topology_section, "CHARGE_BETA", l_val=topology%charge_beta)
65 10950 : CALL section_vals_val_get(topology_section, "CHARGE_EXTENDED", l_val=topology%charge_extended)
66 43800 : ival = COUNT([topology%charge_occup, topology%charge_beta, topology%charge_extended])
67 10950 : IF (ival > 1) THEN
68 0 : CPABORT("Only one between <CHARGE_OCCUP,CHARGE_BETA,CHARGE_EXTENDED> can be defined! ")
69 : END IF
70 10950 : CALL section_vals_val_get(topology_section, "PARA_RES", l_val=topology%para_res)
71 10950 : CALL section_vals_val_get(topology_section, "GENERATE%REORDER", l_val=topology%reorder_atom)
72 10950 : CALL section_vals_val_get(topology_section, "GENERATE%CREATE_MOLECULES", l_val=topology%create_molecules)
73 10950 : CALL section_vals_val_get(topology_section, "MOL_CHECK", l_val=topology%molecules_check)
74 10950 : CALL section_vals_val_get(topology_section, "USE_G96_VELOCITY", l_val=topology%use_g96_velocity)
75 10950 : CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=topology%coord_type)
76 12879 : SELECT CASE (topology%coord_type)
77 : CASE (do_coord_off)
78 : ! Do Nothing
79 : CASE DEFAULT
80 1929 : topology%coordinate = .TRUE.
81 10950 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=topology%coord_file_name)
82 : END SELECT
83 10950 : CALL section_vals_val_get(topology_section, "CONN_FILE_FORMAT", i_val=topology%conn_type)
84 11474 : SELECT CASE (topology%conn_type)
85 : CASE (do_conn_off, do_conn_generate, do_conn_mol_set, do_conn_user)
86 : ! Do Nothing
87 : CASE DEFAULT
88 10950 : CALL section_vals_val_get(topology_section, "CONN_FILE_NAME", c_val=topology%conn_file_name)
89 : END SELECT
90 10950 : CALL section_vals_val_get(topology_section, "EXCLUDE_VDW", i_val=topology%exclude_vdw)
91 10950 : CALL section_vals_val_get(topology_section, "EXCLUDE_EI", i_val=topology%exclude_ei)
92 10950 : CALL section_vals_val_get(topology_section, "GENERATE%BONDPARM", i_val=topology%bondparm_type)
93 10950 : CALL section_vals_val_get(topology_section, "GENERATE%BONDPARM_FACTOR", r_val=topology%bondparm_factor)
94 10950 : CALL timestop(handle)
95 10950 : END SUBROUTINE read_topology_section
96 :
97 : ! **************************************************************************************************
98 : !> \brief Read all the distance parameters. Put them in the
99 : !> constraint_distance array.
100 : !> \param topology ...
101 : !> \param colvar_p ...
102 : !> \param constraint_section ...
103 : !> \par History
104 : !> JGH (26-01-2002) Distance parameters are now stored in tables. The position
105 : !> within the table is used as handle for the topology
106 : !> teo Read the CONSTRAINT section within the new input style
107 : !> \author teo
108 : ! **************************************************************************************************
109 10950 : SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section)
110 :
111 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
112 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
113 : TYPE(section_vals_type), POINTER :: constraint_section
114 :
115 : CHARACTER(LEN=default_string_length), &
116 10950 : DIMENSION(:), POINTER :: tmpstringlist
117 : INTEGER :: icolvar, ig, isize, isize_old, itype, &
118 : jg, msize, msize_old, n_rep, ncons, &
119 : nrep
120 10950 : INTEGER, DIMENSION(:), POINTER :: ilist, tmplist
121 : LOGICAL :: explicit
122 10950 : REAL(KIND=dp), DIMENSION(:), POINTER :: rlist
123 : TYPE(constraint_info_type), POINTER :: cons_info
124 : TYPE(section_vals_type), POINTER :: collective_section, fix_atom_section, &
125 : g3x3_section, g4x6_section, &
126 : hbonds_section, vsite_section
127 :
128 10950 : cons_info => topology%cons_info
129 64025 : IF (ASSOCIATED(constraint_section)) THEN
130 10615 : hbonds_section => section_vals_get_subs_vals(constraint_section, "HBONDS")
131 10615 : g3x3_section => section_vals_get_subs_vals(constraint_section, "G3X3")
132 10615 : g4x6_section => section_vals_get_subs_vals(constraint_section, "G4X6")
133 10615 : vsite_section => section_vals_get_subs_vals(constraint_section, "VIRTUAL_SITE")
134 10615 : fix_atom_section => section_vals_get_subs_vals(constraint_section, "FIXED_ATOMS")
135 10615 : collective_section => section_vals_get_subs_vals(constraint_section, "COLLECTIVE")
136 : ! HBONDS
137 10615 : CALL section_vals_get(hbonds_section, explicit=topology%const_hydr)
138 : CALL check_restraint(hbonds_section, &
139 : is_restraint=cons_info%hbonds_restraint, &
140 : k0=cons_info%hbonds_k0, &
141 10615 : label="HBONDS")
142 : ! G3X3
143 10615 : CALL section_vals_get(g3x3_section, explicit=explicit, n_repetition=ncons)
144 10615 : IF (explicit) THEN
145 156 : topology%const_33 = .TRUE.
146 156 : cons_info%nconst_g33 = ncons
147 : !
148 468 : ALLOCATE (cons_info%const_g33_mol(ncons))
149 468 : ALLOCATE (cons_info%const_g33_molname(ncons))
150 312 : ALLOCATE (cons_info%const_g33_a(ncons))
151 312 : ALLOCATE (cons_info%const_g33_b(ncons))
152 312 : ALLOCATE (cons_info%const_g33_c(ncons))
153 468 : ALLOCATE (cons_info%const_g33_dab(ncons))
154 312 : ALLOCATE (cons_info%const_g33_dac(ncons))
155 312 : ALLOCATE (cons_info%const_g33_dbc(ncons))
156 312 : ALLOCATE (cons_info%g33_intermolecular(ncons))
157 312 : ALLOCATE (cons_info%g33_restraint(ncons))
158 312 : ALLOCATE (cons_info%g33_k0(ncons))
159 312 : ALLOCATE (cons_info%g33_exclude_qm(ncons))
160 312 : ALLOCATE (cons_info%g33_exclude_mm(ncons))
161 316 : DO ig = 1, ncons
162 : CALL check_restraint(g3x3_section, &
163 : is_restraint=cons_info%g33_restraint(ig), &
164 : k0=cons_info%g33_k0(ig), &
165 : i_rep_section=ig, &
166 160 : label="G3X3")
167 160 : cons_info%const_g33_mol(ig) = 0
168 160 : cons_info%const_g33_molname(ig) = "UNDEF"
169 : ! Exclude QM or MM
170 : CALL section_vals_val_get(g3x3_section, "EXCLUDE_QM", i_rep_section=ig, &
171 160 : l_val=cons_info%g33_exclude_qm(ig))
172 : CALL section_vals_val_get(g3x3_section, "EXCLUDE_MM", i_rep_section=ig, &
173 160 : l_val=cons_info%g33_exclude_mm(ig))
174 : ! Intramolecular restraint
175 : CALL section_vals_val_get(g3x3_section, "INTERMOLECULAR", i_rep_section=ig, &
176 160 : l_val=cons_info%g33_intermolecular(ig))
177 : ! If it is intramolecular let's unset (in case user did it)
178 : ! the molecule and molname field
179 160 : IF (cons_info%g33_intermolecular(ig)) THEN
180 4 : CALL section_vals_val_unset(g3x3_section, "MOLECULE", i_rep_section=ig)
181 4 : CALL section_vals_val_unset(g3x3_section, "MOLNAME", i_rep_section=ig)
182 : END IF
183 : ! Let's tag to which molecule we want to apply constraints
184 : CALL section_vals_val_get(g3x3_section, "MOLECULE", i_rep_section=ig, &
185 160 : n_rep_val=nrep)
186 160 : IF (nrep /= 0) THEN
187 : CALL section_vals_val_get(g3x3_section, "MOLECULE", i_rep_section=ig, &
188 120 : i_val=cons_info%const_g33_mol(ig))
189 : END IF
190 : CALL section_vals_val_get(g3x3_section, "MOLNAME", i_rep_section=ig, &
191 160 : n_rep_val=nrep)
192 160 : IF (nrep /= 0) THEN
193 : CALL section_vals_val_get(g3x3_section, "MOLNAME", i_rep_section=ig, &
194 36 : c_val=cons_info%const_g33_molname(ig))
195 : END IF
196 160 : IF ((cons_info%const_g33_mol(ig) /= 0) .AND. (cons_info%const_g33_molname(ig) /= "UNDEF")) THEN
197 : CALL cp_abort(__LOCATION__, &
198 : "Invalid G3X3 constraint section "//cp_to_string(ig)//": "// &
199 0 : "check MOLECULE and MOLNAME setup!")
200 : END IF
201 160 : IF ((cons_info%const_g33_mol(ig) == 0) .AND. (cons_info%const_g33_molname(ig) == "UNDEF") .AND. &
202 : (.NOT. cons_info%g33_intermolecular(ig))) THEN
203 : CALL cp_abort(__LOCATION__, &
204 : "Invalid G3X3 constraint section "//cp_to_string(ig)//": "// &
205 0 : "check MOLECULE and MOLNAME setup!")
206 : END IF
207 : CALL section_vals_val_get(g3x3_section, "ATOMS", i_rep_section=ig, &
208 160 : i_vals=ilist)
209 : CALL section_vals_val_get(g3x3_section, "DISTANCES", i_rep_section=ig, &
210 160 : r_vals=rlist)
211 160 : cons_info%const_g33_a(ig) = ilist(1)
212 160 : cons_info%const_g33_b(ig) = ilist(2)
213 160 : cons_info%const_g33_c(ig) = ilist(3)
214 :
215 160 : cons_info%const_g33_dab(ig) = rlist(1)
216 160 : cons_info%const_g33_dac(ig) = rlist(2)
217 636 : cons_info%const_g33_dbc(ig) = rlist(3)
218 : END DO
219 : END IF
220 : ! G4X6
221 10615 : CALL section_vals_get(g4x6_section, explicit=explicit, n_repetition=ncons)
222 10615 : IF (explicit) THEN
223 16 : topology%const_46 = .TRUE.
224 16 : cons_info%nconst_g46 = ncons
225 : !
226 48 : ALLOCATE (cons_info%const_g46_mol(ncons))
227 48 : ALLOCATE (cons_info%const_g46_molname(ncons))
228 32 : ALLOCATE (cons_info%const_g46_a(ncons))
229 32 : ALLOCATE (cons_info%const_g46_b(ncons))
230 32 : ALLOCATE (cons_info%const_g46_c(ncons))
231 32 : ALLOCATE (cons_info%const_g46_d(ncons))
232 48 : ALLOCATE (cons_info%const_g46_dab(ncons))
233 32 : ALLOCATE (cons_info%const_g46_dac(ncons))
234 32 : ALLOCATE (cons_info%const_g46_dbc(ncons))
235 32 : ALLOCATE (cons_info%const_g46_dad(ncons))
236 32 : ALLOCATE (cons_info%const_g46_dbd(ncons))
237 32 : ALLOCATE (cons_info%const_g46_dcd(ncons))
238 32 : ALLOCATE (cons_info%g46_intermolecular(ncons))
239 32 : ALLOCATE (cons_info%g46_restraint(ncons))
240 32 : ALLOCATE (cons_info%g46_k0(ncons))
241 32 : ALLOCATE (cons_info%g46_exclude_qm(ncons))
242 32 : ALLOCATE (cons_info%g46_exclude_mm(ncons))
243 32 : DO ig = 1, ncons
244 : CALL check_restraint(g4x6_section, &
245 : is_restraint=cons_info%g46_restraint(ig), &
246 : k0=cons_info%g46_k0(ig), &
247 : i_rep_section=ig, &
248 16 : label="G4X6")
249 16 : cons_info%const_g46_mol(ig) = 0
250 16 : cons_info%const_g46_molname(ig) = "UNDEF"
251 : ! Exclude QM or MM
252 : CALL section_vals_val_get(g4x6_section, "EXCLUDE_QM", i_rep_section=ig, &
253 16 : l_val=cons_info%g46_exclude_qm(ig))
254 : CALL section_vals_val_get(g4x6_section, "EXCLUDE_MM", i_rep_section=ig, &
255 16 : l_val=cons_info%g46_exclude_mm(ig))
256 : ! Intramolecular restraint
257 : CALL section_vals_val_get(g4x6_section, "INTERMOLECULAR", i_rep_section=ig, &
258 16 : l_val=cons_info%g46_intermolecular(ig))
259 : ! If it is intramolecular let's unset (in case user did it)
260 : ! the molecule and molname field
261 16 : IF (cons_info%g46_intermolecular(ig)) THEN
262 4 : CALL section_vals_val_unset(g4x6_section, "MOLECULE", i_rep_section=ig)
263 4 : CALL section_vals_val_unset(g4x6_section, "MOLNAME", i_rep_section=ig)
264 : END IF
265 : ! Let's tag to which molecule we want to apply constraints
266 : CALL section_vals_val_get(g4x6_section, "MOLECULE", i_rep_section=ig, &
267 16 : n_rep_val=nrep)
268 16 : IF (nrep /= 0) THEN
269 : CALL section_vals_val_get(g4x6_section, "MOLECULE", i_rep_section=ig, &
270 8 : i_val=cons_info%const_g46_mol(ig))
271 : END IF
272 : CALL section_vals_val_get(g4x6_section, "MOLNAME", i_rep_section=ig, &
273 16 : n_rep_val=nrep)
274 16 : IF (nrep /= 0) THEN
275 : CALL section_vals_val_get(g4x6_section, "MOLNAME", i_rep_section=ig, &
276 4 : c_val=cons_info%const_g46_molname(ig))
277 : END IF
278 16 : IF ((cons_info%const_g46_mol(ig) /= 0) .AND. (cons_info%const_g46_molname(ig) /= "UNDEF")) THEN
279 : CALL cp_abort(__LOCATION__, &
280 : "Invalid G4X6 constraint section "//cp_to_string(ig)//": "// &
281 0 : "check MOLECULE and MOLNAME setup!")
282 : END IF
283 16 : IF ((cons_info%const_g46_mol(ig) == 0) .AND. (cons_info%const_g46_molname(ig) == "UNDEF") .AND. &
284 : (.NOT. cons_info%g46_intermolecular(ig))) THEN
285 : CALL cp_abort(__LOCATION__, &
286 : "Invalid G4X6 constraint section "//cp_to_string(ig)//": "// &
287 0 : "check MOLECULE and MOLNAME setup!")
288 : END IF
289 : CALL section_vals_val_get(g4x6_section, "ATOMS", i_rep_section=ig, &
290 16 : i_vals=ilist)
291 : CALL section_vals_val_get(g4x6_section, "DISTANCES", i_rep_section=ig, &
292 16 : r_vals=rlist)
293 16 : cons_info%const_g46_a(ig) = ilist(1)
294 16 : cons_info%const_g46_b(ig) = ilist(2)
295 16 : cons_info%const_g46_c(ig) = ilist(3)
296 16 : cons_info%const_g46_d(ig) = ilist(4)
297 16 : cons_info%const_g46_dab(ig) = rlist(1)
298 16 : cons_info%const_g46_dac(ig) = rlist(2)
299 16 : cons_info%const_g46_dad(ig) = rlist(3)
300 16 : cons_info%const_g46_dbc(ig) = rlist(4)
301 16 : cons_info%const_g46_dbd(ig) = rlist(5)
302 64 : cons_info%const_g46_dcd(ig) = rlist(6)
303 : END DO
304 : END IF
305 : ! virtual
306 10615 : CALL section_vals_get(vsite_section, explicit=explicit, n_repetition=ncons)
307 10615 : IF (explicit) THEN
308 8 : topology%const_vsite = .TRUE.
309 8 : cons_info%nconst_vsite = ncons
310 : !
311 24 : ALLOCATE (cons_info%const_vsite_mol(ncons))
312 24 : ALLOCATE (cons_info%const_vsite_molname(ncons))
313 16 : ALLOCATE (cons_info%const_vsite_a(ncons))
314 16 : ALLOCATE (cons_info%const_vsite_b(ncons))
315 16 : ALLOCATE (cons_info%const_vsite_c(ncons))
316 16 : ALLOCATE (cons_info%const_vsite_d(ncons))
317 24 : ALLOCATE (cons_info%const_vsite_wbc(ncons))
318 16 : ALLOCATE (cons_info%const_vsite_wdc(ncons))
319 16 : ALLOCATE (cons_info%vsite_intermolecular(ncons))
320 16 : ALLOCATE (cons_info%vsite_restraint(ncons))
321 16 : ALLOCATE (cons_info%vsite_k0(ncons))
322 16 : ALLOCATE (cons_info%vsite_exclude_qm(ncons))
323 16 : ALLOCATE (cons_info%vsite_exclude_mm(ncons))
324 16 : DO ig = 1, ncons
325 : CALL check_restraint(vsite_section, &
326 : is_restraint=cons_info%vsite_restraint(ig), &
327 : k0=cons_info%vsite_k0(ig), &
328 : i_rep_section=ig, &
329 8 : label="Virtual_SITE")
330 8 : cons_info%const_vsite_mol(ig) = 0
331 8 : cons_info%const_vsite_molname(ig) = "UNDEF"
332 : ! Exclude QM or MM
333 : CALL section_vals_val_get(vsite_section, "EXCLUDE_QM", i_rep_section=ig, &
334 8 : l_val=cons_info%vsite_exclude_qm(ig))
335 : CALL section_vals_val_get(vsite_section, "EXCLUDE_MM", i_rep_section=ig, &
336 8 : l_val=cons_info%vsite_exclude_mm(ig))
337 : ! Intramolecular restraint
338 : CALL section_vals_val_get(vsite_section, "INTERMOLECULAR", i_rep_section=ig, &
339 8 : l_val=cons_info%vsite_intermolecular(ig))
340 : ! If it is intramolecular let's unset (in case user did it)
341 : ! the molecule and molname field
342 8 : IF (cons_info%vsite_intermolecular(ig)) THEN
343 0 : CALL section_vals_val_unset(vsite_section, "MOLECULE", i_rep_section=ig)
344 0 : CALL section_vals_val_unset(vsite_section, "MOLNAME", i_rep_section=ig)
345 : END IF
346 : ! Let's tag to which molecule we want to apply constraints
347 : CALL section_vals_val_get(vsite_section, "MOLECULE", i_rep_section=ig, &
348 8 : n_rep_val=nrep)
349 8 : IF (nrep /= 0) THEN
350 : CALL section_vals_val_get(vsite_section, "MOLECULE", i_rep_section=ig, &
351 8 : i_val=cons_info%const_vsite_mol(ig))
352 : END IF
353 : CALL section_vals_val_get(vsite_section, "MOLNAME", i_rep_section=ig, &
354 8 : n_rep_val=nrep)
355 8 : IF (nrep /= 0) THEN
356 : CALL section_vals_val_get(vsite_section, "MOLNAME", i_rep_section=ig, &
357 0 : c_val=cons_info%const_vsite_molname(ig))
358 : END IF
359 8 : IF ((cons_info%const_vsite_mol(ig) /= 0) .AND. (cons_info%const_vsite_molname(ig) /= "UNDEF")) THEN
360 : CALL cp_abort(__LOCATION__, &
361 : "Invalid VIRTUAL_SITE constraint section "//cp_to_string(ig)//": "// &
362 0 : "check MOLECULE and MOLNAME setup!")
363 : END IF
364 8 : IF ((cons_info%const_vsite_mol(ig) == 0) .AND. (cons_info%const_vsite_molname(ig) == "UNDEF") .AND. &
365 : (.NOT. cons_info%vsite_intermolecular(ig))) THEN
366 : CALL cp_abort(__LOCATION__, &
367 : "Invalid VIRTUAL_SITE constraint section "//cp_to_string(ig)//": "// &
368 0 : "check MOLECULE and MOLNAME setup!")
369 : END IF
370 : CALL section_vals_val_get(vsite_section, "ATOMS", i_rep_section=ig, &
371 8 : i_vals=ilist)
372 : CALL section_vals_val_get(vsite_section, "PARAMETERS", i_rep_section=ig, &
373 8 : r_vals=rlist)
374 8 : cons_info%const_vsite_a(ig) = ilist(1)
375 8 : cons_info%const_vsite_b(ig) = ilist(2)
376 8 : cons_info%const_vsite_c(ig) = ilist(3)
377 8 : cons_info%const_vsite_d(ig) = ilist(4)
378 8 : cons_info%const_vsite_wbc(ig) = rlist(1)
379 32 : cons_info%const_vsite_wdc(ig) = rlist(2)
380 : END DO
381 : END IF
382 : ! FIXED ATOMS
383 10615 : CALL section_vals_get(fix_atom_section, explicit=explicit, n_repetition=ncons)
384 10615 : IF (explicit) THEN
385 110 : NULLIFY (tmplist, tmpstringlist)
386 110 : isize = 0
387 110 : msize = 0
388 110 : ALLOCATE (cons_info%fixed_atoms(isize))
389 110 : ALLOCATE (cons_info%fixed_type(isize))
390 110 : ALLOCATE (cons_info%fixed_restraint(isize))
391 110 : ALLOCATE (cons_info%fixed_k0(isize))
392 110 : ALLOCATE (cons_info%fixed_molnames(msize))
393 110 : ALLOCATE (cons_info%fixed_mol_type(isize))
394 110 : ALLOCATE (cons_info%fixed_mol_restraint(msize))
395 110 : ALLOCATE (cons_info%fixed_mol_k0(msize))
396 330 : ALLOCATE (cons_info%fixed_exclude_qm(ncons))
397 220 : ALLOCATE (cons_info%fixed_exclude_mm(ncons))
398 246 : DO ig = 1, ncons
399 136 : isize_old = isize
400 136 : msize_old = msize
401 : CALL section_vals_val_get(fix_atom_section, "COMPONENTS_TO_FIX", i_rep_section=ig, &
402 136 : i_val=itype)
403 : CALL section_vals_val_get(fix_atom_section, "LIST", i_rep_section=ig, &
404 136 : n_rep_val=n_rep)
405 254 : DO jg = 1, n_rep
406 : CALL section_vals_val_get(fix_atom_section, "LIST", i_rep_section=ig, &
407 118 : i_rep_val=jg, i_vals=tmplist)
408 118 : CALL reallocate(cons_info%fixed_atoms, 1, isize + SIZE(tmplist))
409 20142 : cons_info%fixed_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
410 118 : CALL reallocate(cons_info%fixed_restraint, 1, isize + SIZE(tmplist))
411 118 : CALL reallocate(cons_info%fixed_k0, 1, isize + SIZE(tmplist))
412 118 : CALL reallocate(cons_info%fixed_type, 1, isize + SIZE(tmplist))
413 10130 : cons_info%fixed_type(isize + 1:isize + SIZE(tmplist)) = itype
414 254 : isize = SIZE(cons_info%fixed_atoms)
415 : END DO
416 : !Check for restraints
417 136 : IF ((isize - isize_old) > 0) THEN
418 : CALL check_restraint(fix_atom_section, &
419 : is_restraint=cons_info%fixed_restraint(isize_old + 1), &
420 : k0=cons_info%fixed_k0(isize_old + 1), &
421 : i_rep_section=ig, &
422 112 : label="FIXED ATOM")
423 10124 : cons_info%fixed_restraint(isize_old + 1:isize) = cons_info%fixed_restraint(isize_old + 1)
424 10124 : cons_info%fixed_k0(isize_old + 1:isize) = cons_info%fixed_k0(isize_old + 1)
425 : END IF
426 : CALL section_vals_val_get(fix_atom_section, "MOLNAME", i_rep_section=ig, &
427 136 : n_rep_val=n_rep)
428 136 : IF (n_rep /= 0) THEN
429 12 : DO jg = 1, n_rep
430 : CALL section_vals_val_get(fix_atom_section, "MOLNAME", i_rep_section=ig, &
431 6 : i_rep_val=jg, c_vals=tmpstringlist)
432 6 : CALL reallocate(cons_info%fixed_molnames, 1, msize + SIZE(tmpstringlist, 1))
433 6 : CALL reallocate(cons_info%fixed_mol_type, 1, msize + SIZE(tmpstringlist, 1))
434 6 : CALL reallocate(cons_info%fixed_mol_restraint, 1, msize + SIZE(tmpstringlist, 1))
435 6 : CALL reallocate(cons_info%fixed_mol_k0, 1, msize + SIZE(tmpstringlist, 1))
436 18 : cons_info%fixed_molnames(msize + 1:msize + SIZE(tmpstringlist, 1)) = tmpstringlist
437 12 : cons_info%fixed_mol_type(msize + 1:msize + SIZE(tmpstringlist, 1)) = itype
438 12 : msize = SIZE(cons_info%fixed_molnames)
439 : END DO
440 : ! Exclude QM or MM work only if defined MOLNAME
441 6 : CALL reallocate(cons_info%fixed_exclude_qm, 1, msize)
442 6 : CALL reallocate(cons_info%fixed_exclude_mm, 1, msize)
443 : CALL section_vals_val_get(fix_atom_section, "EXCLUDE_QM", i_rep_section=ig, &
444 6 : l_val=cons_info%fixed_exclude_qm(msize_old + 1))
445 : CALL section_vals_val_get(fix_atom_section, "EXCLUDE_MM", i_rep_section=ig, &
446 6 : l_val=cons_info%fixed_exclude_mm(msize_old + 1))
447 12 : cons_info%fixed_exclude_qm(msize_old + 1:msize) = cons_info%fixed_exclude_qm(msize_old + 1)
448 12 : cons_info%fixed_exclude_mm(msize_old + 1:msize) = cons_info%fixed_exclude_mm(msize_old + 1)
449 : END IF
450 : !Check for restraints
451 136 : IF (n_rep /= 0) THEN
452 : CALL check_restraint(fix_atom_section, &
453 : is_restraint=cons_info%fixed_mol_restraint(msize_old + 1), &
454 : k0=cons_info%fixed_mol_k0(msize_old + 1), &
455 : i_rep_section=ig, &
456 6 : label="FIXED ATOM")
457 12 : cons_info%fixed_mol_restraint(msize_old + 1:msize) = cons_info%fixed_mol_restraint(msize_old + 1)
458 12 : cons_info%fixed_mol_k0(msize_old + 1:msize) = cons_info%fixed_mol_k0(msize_old + 1)
459 : END IF
460 : CALL section_vals_val_get(fix_atom_section, "MM_SUBSYS", i_rep_section=ig, &
461 136 : n_rep_val=nrep, explicit=explicit)
462 136 : IF (nrep == 1 .AND. explicit) THEN
463 16 : CPASSERT(cons_info%freeze_mm == do_constr_none)
464 : CALL section_vals_val_get(fix_atom_section, "MM_SUBSYS", i_val=cons_info%freeze_mm, &
465 16 : i_rep_section=ig)
466 16 : cons_info%freeze_mm_type = itype
467 : END IF
468 : CALL section_vals_val_get(fix_atom_section, "QM_SUBSYS", i_rep_section=ig, &
469 136 : n_rep_val=nrep, explicit=explicit)
470 136 : IF (nrep == 1 .AND. explicit) THEN
471 2 : CPASSERT(cons_info%freeze_qm == do_constr_none)
472 : CALL section_vals_val_get(fix_atom_section, "QM_SUBSYS", i_val=cons_info%freeze_qm, &
473 2 : i_rep_section=ig)
474 2 : cons_info%freeze_qm_type = itype
475 : END IF
476 136 : IF (cons_info%freeze_mm /= do_constr_none) THEN
477 : CALL check_restraint(fix_atom_section, &
478 : is_restraint=cons_info%fixed_mm_restraint, &
479 : k0=cons_info%fixed_mm_k0, &
480 : i_rep_section=ig, &
481 28 : label="FIXED ATOM")
482 : END IF
483 790 : IF (cons_info%freeze_qm /= do_constr_none) THEN
484 : CALL check_restraint(fix_atom_section, &
485 : is_restraint=cons_info%fixed_qm_restraint, &
486 : k0=cons_info%fixed_qm_k0, &
487 : i_rep_section=ig, &
488 2 : label="FIXED ATOM")
489 : END IF
490 :
491 : END DO
492 : IF ((isize /= 0) .OR. (msize /= 0) .OR. &
493 110 : (cons_info%freeze_mm /= do_constr_none) .OR. &
494 : (cons_info%freeze_qm /= do_constr_none)) THEN
495 110 : topology%const_atom = .TRUE.
496 : END IF
497 : END IF
498 : ! Collective Constraints
499 10615 : CALL section_vals_get(collective_section, explicit=explicit, n_repetition=ncons)
500 10615 : IF (explicit) THEN
501 120 : topology%const_colv = .TRUE.
502 390 : DO ig = 1, ncons
503 270 : CALL section_vals_val_get(collective_section, "COLVAR", i_rep_section=ig, i_val=icolvar)
504 390 : IF (icolvar > SIZE(colvar_p)) THEN
505 0 : CPABORT("More collective constraints than collective variables specified.")
506 : END IF
507 : END DO
508 120 : cons_info%nconst_colv = ncons
509 360 : ALLOCATE (cons_info%const_colv_mol(ncons))
510 360 : ALLOCATE (cons_info%const_colv_molname(ncons))
511 360 : ALLOCATE (cons_info%const_colv_target(ncons))
512 240 : ALLOCATE (cons_info%const_colv_target_growth(ncons))
513 510 : ALLOCATE (cons_info%colvar_set(ncons))
514 240 : ALLOCATE (cons_info%colv_intermolecular(ncons))
515 240 : ALLOCATE (cons_info%colv_restraint(ncons))
516 240 : ALLOCATE (cons_info%colv_k0(ncons))
517 240 : ALLOCATE (cons_info%colv_exclude_qm(ncons))
518 240 : ALLOCATE (cons_info%colv_exclude_mm(ncons))
519 390 : DO ig = 1, ncons
520 : CALL check_restraint(collective_section, &
521 : is_restraint=cons_info%colv_restraint(ig), &
522 : k0=cons_info%colv_k0(ig), &
523 : i_rep_section=ig, &
524 270 : label="COLLECTIVE")
525 270 : cons_info%const_colv_mol(ig) = 0
526 270 : cons_info%const_colv_molname(ig) = "UNDEF"
527 : ! Exclude QM or MM
528 : CALL section_vals_val_get(collective_section, "EXCLUDE_QM", i_rep_section=ig, &
529 270 : l_val=cons_info%colv_exclude_qm(ig))
530 : CALL section_vals_val_get(collective_section, "EXCLUDE_MM", i_rep_section=ig, &
531 270 : l_val=cons_info%colv_exclude_mm(ig))
532 : ! Intramolecular restraint
533 : CALL section_vals_val_get(collective_section, "INTERMOLECULAR", i_rep_section=ig, &
534 270 : l_val=cons_info%colv_intermolecular(ig))
535 : ! If it is intramolecular let's unset (in case user did it)
536 : ! the molecule and molname field
537 270 : IF (cons_info%colv_intermolecular(ig)) THEN
538 66 : CALL section_vals_val_unset(collective_section, "MOLECULE", i_rep_section=ig)
539 66 : CALL section_vals_val_unset(collective_section, "MOLNAME", i_rep_section=ig)
540 : END IF
541 : ! Let's tag to which molecule we want to apply constraints
542 : CALL section_vals_val_get(collective_section, "MOLECULE", i_rep_section=ig, &
543 270 : n_rep_val=nrep)
544 270 : IF (nrep /= 0) THEN
545 : CALL section_vals_val_get(collective_section, "MOLECULE", i_rep_section=ig, &
546 160 : i_val=cons_info%const_colv_mol(ig))
547 : END IF
548 : CALL section_vals_val_get(collective_section, "MOLNAME", i_rep_section=ig, &
549 270 : n_rep_val=nrep)
550 270 : IF (nrep /= 0) THEN
551 : CALL section_vals_val_get(collective_section, "MOLNAME", i_rep_section=ig, &
552 44 : c_val=cons_info%const_colv_molname(ig))
553 : END IF
554 270 : IF (((cons_info%const_colv_mol(ig) /= 0) .AND. (cons_info%const_colv_molname(ig) /= "UNDEF"))) THEN
555 0 : CPABORT("Both MOLNAME and MOLECULE specified for CONSTRAINT section. ")
556 : END IF
557 270 : IF ((cons_info%const_colv_mol(ig) == 0) .AND. (cons_info%const_colv_molname(ig) == "UNDEF") .AND. &
558 : (.NOT. cons_info%colv_intermolecular(ig))) THEN
559 : CALL cp_abort(__LOCATION__, &
560 : "Constraint section error: you have to specify at least one of the "// &
561 0 : "following keywords: MOLECULE, MOLNAME or INTERMOLECULAR! ")
562 : END IF
563 270 : NULLIFY (cons_info%colvar_set(ig)%colvar)
564 : CALL section_vals_val_get(collective_section, "COLVAR", i_rep_section=ig, &
565 270 : i_val=icolvar)
566 : CALL colvar_clone(cons_info%colvar_set(ig)%colvar, &
567 270 : colvar_p(icolvar)%colvar)
568 : CALL section_vals_val_get(collective_section, "TARGET", &
569 270 : n_rep_val=n_rep, i_rep_section=ig)
570 270 : IF (n_rep /= 0) THEN
571 : CALL section_vals_val_get(collective_section, "TARGET", &
572 168 : r_val=cons_info%const_colv_target(ig), i_rep_section=ig)
573 : ELSE
574 102 : cons_info%const_colv_target(ig) = -HUGE(0.0_dp)
575 : END IF
576 : CALL section_vals_val_get(collective_section, "TARGET_GROWTH", &
577 1470 : r_val=cons_info%const_colv_target_growth(ig), i_rep_section=ig)
578 : END DO
579 : END IF
580 : END IF
581 :
582 10950 : END SUBROUTINE read_constraints_section
583 :
584 : ! **************************************************************************************************
585 : !> \brief Reads input and decides if apply restraints instead of constraints
586 : !> \param cons_section ...
587 : !> \param is_restraint ...
588 : !> \param k0 ...
589 : !> \param i_rep_section ...
590 : !> \param label ...
591 : !> \author teo
592 : ! **************************************************************************************************
593 22434 : SUBROUTINE check_restraint(cons_section, is_restraint, k0, i_rep_section, label)
594 : TYPE(section_vals_type), POINTER :: cons_section
595 : LOGICAL, INTENT(OUT) :: is_restraint
596 : REAL(KIND=dp), INTENT(OUT) :: k0
597 : INTEGER, INTENT(IN), OPTIONAL :: i_rep_section
598 : CHARACTER(LEN=*), INTENT(IN) :: label
599 :
600 : CHARACTER(LEN=default_string_length) :: nlabel
601 : INTEGER :: output_unit
602 : LOGICAL :: explicit
603 : TYPE(section_vals_type), POINTER :: restraint_section
604 :
605 11217 : is_restraint = .FALSE.
606 11217 : output_unit = cp_logger_get_default_io_unit()
607 11217 : CALL section_vals_get(cons_section, explicit=explicit)
608 11217 : IF (explicit) THEN
609 : restraint_section => section_vals_get_subs_vals(cons_section, "RESTRAINT", &
610 618 : i_rep_section=i_rep_section)
611 618 : CALL section_vals_get(restraint_section, explicit=is_restraint)
612 618 : IF (is_restraint) THEN
613 124 : CALL section_vals_val_get(restraint_section, "K", r_val=k0)
614 124 : IF (output_unit > 0) THEN
615 64 : nlabel = cp_to_string(i_rep_section)
616 : WRITE (output_unit, FMT='(T2,"RESTRAINT|",1X,A,F9.6)') &
617 : "Active restraint on "//label//" section Nr."// &
618 64 : TRIM(nlabel)//". K [a.u.]=", k0
619 : END IF
620 : END IF
621 : END IF
622 11217 : END SUBROUTINE check_restraint
623 :
624 : END MODULE topology_input
625 :
|
