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