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 : MODULE topology_gromos
10 : USE cp_log_handling, ONLY: cp_get_default_logger,&
11 : cp_logger_type
12 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
13 : cp_print_key_unit_nr
14 : USE cp_parser_methods, ONLY: parser_get_next_line,&
15 : parser_get_object,&
16 : parser_search_string
17 : USE cp_parser_types, ONLY: cp_parser_type,&
18 : parser_create,&
19 : parser_release
20 : USE cp_units, ONLY: cp_unit_to_cp2k
21 : USE input_cp2k_restarts_util, ONLY: section_velocity_val_set
22 : USE input_section_types, ONLY: section_get_rval,&
23 : section_vals_get_subs_vals,&
24 : section_vals_type
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE memory_utilities, ONLY: reallocate
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE string_table, ONLY: s2s,&
30 : str2id
31 : USE topology_types, ONLY: atom_info_type,&
32 : connectivity_info_type,&
33 : topology_parameters_type
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_gromos'
39 :
40 : PRIVATE
41 : PUBLIC :: read_topology_gromos, read_coordinate_g96
42 :
43 : CONTAINS
44 :
45 : ! **************************************************************************************************
46 : !> \brief Read GROMOS topology file
47 : !> \param file_name ...
48 : !> \param topology ...
49 : !> \param para_env ...
50 : !> \param subsys_section ...
51 : ! **************************************************************************************************
52 140 : SUBROUTINE read_topology_gromos(file_name, topology, para_env, subsys_section)
53 : CHARACTER(LEN=*), INTENT(IN) :: file_name
54 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
55 : TYPE(mp_para_env_type), POINTER :: para_env
56 : TYPE(section_vals_type), POINTER :: subsys_section
57 :
58 : CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_gromos'
59 :
60 : CHARACTER :: ctemp
61 : CHARACTER(LEN=default_string_length) :: label, string
62 : CHARACTER(LEN=default_string_length), &
63 : DIMENSION(21) :: avail_section
64 14 : CHARACTER(LEN=default_string_length), POINTER :: namearray1(:), namearray2(:)
65 : INTEGER :: begin, handle, i, iatom, ibond, icon, ii(50), index_now, iresid, isolvent, itemp, &
66 : itype, iw, jatom, natom, natom_prev, nbond, nbond_prev, ncon, nsolute, nsolvent, ntype, &
67 : offset, offset2, stat
68 14 : INTEGER, POINTER :: ba(:), bb(:), na(:)
69 : LOGICAL :: found
70 : REAL(KIND=dp) :: ftemp
71 14 : REAL(KIND=dp), POINTER :: ac(:), am(:)
72 : TYPE(atom_info_type), POINTER :: atom_info
73 : TYPE(connectivity_info_type), POINTER :: conn_info
74 : TYPE(cp_logger_type), POINTER :: logger
75 : TYPE(cp_parser_type) :: parser
76 :
77 14 : NULLIFY (logger)
78 14 : NULLIFY (namearray1, namearray2)
79 28 : logger => cp_get_default_logger()
80 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GTOP_INFO", &
81 14 : extension=".subsysLog")
82 14 : CALL timeset(routineN, handle)
83 :
84 14 : avail_section(1) = "TITLE"
85 14 : avail_section(2) = "TOPPHYSCON"
86 14 : avail_section(3) = "TOPVERSION"
87 14 : avail_section(4) = "ATOMTYPENAME"
88 14 : avail_section(5) = "RESNAME"
89 14 : avail_section(6) = "SOLUTEATOM"
90 14 : avail_section(7) = "BONDTYPE"
91 14 : avail_section(8) = "BONDH"
92 14 : avail_section(9) = "BOND"
93 14 : avail_section(10) = "BONDANGLETYPE"
94 14 : avail_section(11) = "BONDANGLEH"
95 14 : avail_section(12) = "BONDANGLE"
96 14 : avail_section(13) = "IMPDIHEDRALTYPE"
97 14 : avail_section(14) = "IMPDIHEDRALH"
98 14 : avail_section(15) = "IMPDIHEDRAL"
99 14 : avail_section(16) = "DIHEDRALTYPE"
100 14 : avail_section(17) = "DIHEDRALH"
101 14 : avail_section(18) = "DIHEDRAL"
102 14 : avail_section(19) = "LJPARAMETERS"
103 14 : avail_section(20) = "SOLVENTATOM"
104 14 : avail_section(21) = "SOLVENTCONSTR"
105 :
106 14 : atom_info => topology%atom_info
107 14 : conn_info => topology%conn_info
108 :
109 14 : natom_prev = 0
110 14 : IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
111 : ! TITLE SECTION
112 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TITLE section'
113 14 : CALL parser_create(parser, file_name, para_env=para_env)
114 14 : label = TRIM(avail_section(1))
115 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
116 14 : IF (found) THEN
117 : DO
118 56 : CALL parser_get_next_line(parser, 1)
119 56 : CALL parser_get_object(parser, string, string_length=80)
120 56 : IF (string == TRIM("END")) EXIT
121 56 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string)
122 : END DO
123 : END IF
124 14 : CALL parser_release(parser)
125 : ! TOPPHYSCON SECTION
126 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPPHYSCON section'
127 14 : CALL parser_create(parser, file_name, para_env=para_env)
128 14 : label = TRIM(avail_section(2))
129 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
130 14 : IF (found) THEN
131 : DO
132 42 : CALL parser_get_next_line(parser, 1)
133 42 : CALL parser_get_object(parser, string)
134 42 : IF (string == TRIM("END")) EXIT
135 42 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string)
136 : END DO
137 : END IF
138 14 : CALL parser_release(parser)
139 : ! TOPVERSION SECTION
140 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPVERSION section'
141 14 : CALL parser_create(parser, file_name, para_env=para_env)
142 14 : label = TRIM(avail_section(3))
143 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
144 14 : IF (found) THEN
145 : DO
146 28 : CALL parser_get_next_line(parser, 1)
147 28 : CALL parser_get_object(parser, string)
148 28 : IF (string == TRIM("END")) EXIT
149 28 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string)
150 : END DO
151 : END IF
152 14 : CALL parser_release(parser)
153 : ! ATOMTYPENAME SECTION
154 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section'
155 14 : CALL parser_create(parser, file_name, para_env=para_env)
156 14 : label = TRIM(avail_section(4))
157 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
158 14 : IF (found) THEN
159 14 : CALL parser_get_next_line(parser, 1)
160 14 : CALL parser_get_object(parser, ntype)
161 14 : CALL reallocate(namearray1, 1, ntype)
162 104 : DO itype = 1, ntype
163 90 : CALL parser_get_next_line(parser, 1)
164 90 : CALL parser_get_object(parser, namearray1(itype))
165 104 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray1(itype))
166 : END DO
167 : END IF
168 14 : CALL parser_release(parser)
169 : ! RESNAME SECTION
170 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the RESNAME section'
171 14 : CALL parser_create(parser, file_name, para_env=para_env)
172 14 : label = TRIM(avail_section(5))
173 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
174 14 : IF (found) THEN
175 14 : CALL parser_get_next_line(parser, 1)
176 14 : CALL parser_get_object(parser, ntype)
177 14 : CALL reallocate(namearray2, 1, ntype)
178 28 : DO itype = 1, ntype
179 14 : CALL parser_get_next_line(parser, 1)
180 14 : CALL parser_get_object(parser, namearray2(itype))
181 28 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray2(itype))
182 : END DO
183 : END IF
184 14 : CALL parser_release(parser)
185 : ! SOLUTEATOM SECTION
186 14 : iresid = 1
187 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLUTEATOM section'
188 14 : CALL parser_create(parser, file_name, para_env=para_env)
189 14 : label = TRIM(avail_section(6))
190 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
191 14 : IF (found) THEN
192 14 : CALL parser_get_next_line(parser, 1)
193 14 : CALL parser_get_object(parser, natom)
194 14 : CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
195 14 : CALL reallocate(atom_info%resid, 1, natom_prev + natom)
196 14 : CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
197 14 : CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
198 14 : CALL reallocate(atom_info%id_element, 1, natom_prev + natom)
199 14 : CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
200 14 : CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
201 14 : CALL parser_get_next_line(parser, 1)
202 136 : DO iatom = 1, natom
203 122 : index_now = iatom + natom_prev
204 122 : CALL parser_get_object(parser, itemp)
205 122 : CALL parser_get_object(parser, itemp)
206 122 : atom_info%resid(index_now) = itemp
207 122 : atom_info%id_molname(index_now) = str2id(s2s(namearray2(itemp)))
208 122 : atom_info%id_resname(index_now) = str2id(s2s(namearray2(itemp)))
209 122 : CALL parser_get_object(parser, string)
210 122 : CALL parser_get_object(parser, itemp)
211 122 : atom_info%id_atmname(index_now) = str2id(s2s(namearray1(itemp)))
212 122 : atom_info%id_element(index_now) = str2id(s2s(namearray1(itemp)))
213 122 : CALL parser_get_object(parser, atom_info%atm_mass(index_now))
214 122 : CALL parser_get_object(parser, atom_info%atm_charge(index_now))
215 122 : CALL parser_get_object(parser, itemp)
216 122 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLUTEATOM INFO HERE!"
217 122 : CALL parser_get_object(parser, ntype)
218 6222 : DO i = 1, 50
219 6222 : ii(i) = -1
220 : END DO
221 122 : stat = -1
222 122 : begin = 1
223 122 : IF (ntype /= 0) THEN
224 : DO
225 108 : IF (begin .EQ. 1) THEN
226 100 : READ (parser%input_line, iostat=stat, FMT=*) itemp, itemp, ctemp, itemp, ftemp, ftemp, &
227 200 : itemp, itemp, (ii(i), i=begin, ntype)
228 8 : ELSE IF (begin .GT. 1) THEN
229 8 : CALL parser_get_next_line(parser, 1)
230 8 : READ (parser%input_line, iostat=stat, FMT=*) (ii(i), i=begin, ntype)
231 : END IF
232 438 : DO i = ntype, 1, -1
233 438 : IF (ii(i) .LT. 0) THEN
234 16 : begin = i
235 : END IF
236 : END DO
237 108 : IF (stat .EQ. 0) EXIT
238 : END DO
239 : END IF
240 : ! 1-4 list
241 122 : CALL parser_get_next_line(parser, 1)
242 122 : CALL parser_get_object(parser, ntype)
243 624 : IF (ntype /= 0) THEN
244 80 : itemp = (itemp - 1)/6 + 1
245 80 : offset = 0
246 80 : IF (ASSOCIATED(conn_info%onfo_a)) offset = SIZE(conn_info%onfo_a)
247 80 : CALL reallocate(conn_info%onfo_a, 1, offset + ntype)
248 80 : CALL reallocate(conn_info%onfo_b, 1, offset + ntype)
249 264 : conn_info%onfo_a(offset + 1:offset + ntype) = index_now
250 4080 : DO i = 1, 50
251 4080 : ii(i) = -1
252 : END DO
253 80 : stat = -1
254 80 : begin = 1
255 : IF (ntype /= 0) THEN
256 : DO
257 80 : IF (begin .EQ. 1) THEN
258 80 : READ (parser%input_line, iostat=stat, FMT=*) itemp, (ii(i), i=begin, ntype)
259 0 : ELSE IF (begin .GT. 1) THEN
260 0 : CALL parser_get_next_line(parser, 1)
261 0 : READ (parser%input_line, iostat=stat, FMT=*) (ii(i), i=begin, ntype)
262 : END IF
263 264 : DO i = ntype, 1, -1
264 264 : IF (ii(i) .LT. 0) THEN
265 0 : begin = i
266 : END IF
267 : END DO
268 80 : IF (stat .EQ. 0) EXIT
269 : END DO
270 264 : DO i = 1, ntype
271 264 : conn_info%onfo_b(offset + i) = ii(i)
272 : END DO
273 : END IF
274 80 : CALL parser_get_next_line(parser, 1)
275 : ELSE
276 42 : CALL parser_get_next_line(parser, 1)
277 : END IF
278 : END DO
279 : END IF
280 14 : nsolute = natom
281 14 : CALL parser_release(parser)
282 :
283 14 : CALL parser_create(parser, file_name, para_env=para_env)
284 : ! BONDH SECTION
285 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDH section'
286 14 : label = TRIM(avail_section(8))
287 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
288 14 : IF (found) THEN
289 14 : CALL parser_get_next_line(parser, 1)
290 14 : CALL parser_get_object(parser, ntype)
291 14 : offset = 0
292 14 : IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
293 14 : CALL reallocate(conn_info%bond_a, 1, offset + ntype)
294 14 : CALL reallocate(conn_info%bond_b, 1, offset + ntype)
295 14 : CALL reallocate(conn_info%bond_type, 1, offset + ntype)
296 82 : DO itype = 1, ntype
297 68 : CALL parser_get_next_line(parser, 1)
298 68 : CALL parser_get_object(parser, conn_info%bond_a(offset + itype))
299 68 : CALL parser_get_object(parser, conn_info%bond_b(offset + itype))
300 68 : CALL parser_get_object(parser, itemp)
301 68 : conn_info%bond_type(offset + itype) = itemp
302 82 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDH INFO HERE!"
303 : END DO
304 82 : conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev
305 82 : conn_info%bond_b(offset + 1:offset + ntype) = conn_info%bond_b(offset + 1:offset + ntype) + natom_prev
306 : END IF
307 : ! BOND SECTION
308 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BOND section'
309 14 : label = TRIM(avail_section(9))
310 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
311 14 : IF (found) THEN
312 14 : CALL parser_get_next_line(parser, 1)
313 14 : CALL parser_get_object(parser, ntype)
314 14 : offset = 0
315 14 : IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
316 14 : CALL reallocate(conn_info%bond_a, 1, offset + ntype)
317 14 : CALL reallocate(conn_info%bond_b, 1, offset + ntype)
318 14 : CALL reallocate(conn_info%bond_type, 1, offset + ntype)
319 54 : DO itype = 1, ntype
320 40 : CALL parser_get_next_line(parser, 1)
321 40 : CALL parser_get_object(parser, conn_info%bond_a(offset + itype))
322 40 : CALL parser_get_object(parser, conn_info%bond_b(offset + itype))
323 40 : CALL parser_get_object(parser, itemp)
324 40 : conn_info%bond_type(offset + itype) = itemp
325 54 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BOND INFO HERE!"
326 : END DO
327 54 : conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev
328 54 : conn_info%bond_b(offset + 1:offset + ntype) = conn_info%bond_b(offset + 1:offset + ntype) + natom_prev
329 : END IF
330 : ! BONDANGLEH SECTION
331 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLEH section'
332 14 : label = TRIM(avail_section(11))
333 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
334 14 : IF (found) THEN
335 14 : CALL parser_get_next_line(parser, 1)
336 14 : CALL parser_get_object(parser, ntype)
337 14 : offset = 0
338 14 : IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a)
339 14 : CALL reallocate(conn_info%theta_a, 1, offset + ntype)
340 14 : CALL reallocate(conn_info%theta_b, 1, offset + ntype)
341 14 : CALL reallocate(conn_info%theta_c, 1, offset + ntype)
342 14 : CALL reallocate(conn_info%theta_type, 1, offset + ntype)
343 124 : DO itype = 1, ntype
344 110 : CALL parser_get_next_line(parser, 1)
345 110 : CALL parser_get_object(parser, conn_info%theta_a(offset + itype))
346 110 : CALL parser_get_object(parser, conn_info%theta_b(offset + itype))
347 110 : CALL parser_get_object(parser, conn_info%theta_c(offset + itype))
348 110 : CALL parser_get_object(parser, itemp)
349 110 : conn_info%theta_type(offset + itype) = itemp
350 124 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLEH INFO HERE!"
351 : END DO
352 124 : conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev
353 124 : conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev
354 124 : conn_info%theta_c(offset + 1:offset + ntype) = conn_info%theta_c(offset + 1:offset + ntype) + natom_prev
355 : END IF
356 : ! BONDANGLE SECTION
357 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLE section'
358 14 : label = TRIM(avail_section(12))
359 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
360 14 : IF (found) THEN
361 14 : CALL parser_get_next_line(parser, 1)
362 14 : CALL parser_get_object(parser, ntype)
363 14 : offset = 0
364 14 : IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a)
365 14 : CALL reallocate(conn_info%theta_a, 1, offset + ntype)
366 14 : CALL reallocate(conn_info%theta_b, 1, offset + ntype)
367 14 : CALL reallocate(conn_info%theta_c, 1, offset + ntype)
368 14 : CALL reallocate(conn_info%theta_type, 1, offset + ntype)
369 62 : DO itype = 1, ntype
370 48 : CALL parser_get_next_line(parser, 1)
371 48 : CALL parser_get_object(parser, conn_info%theta_a(offset + itype))
372 48 : CALL parser_get_object(parser, conn_info%theta_b(offset + itype))
373 48 : CALL parser_get_object(parser, conn_info%theta_c(offset + itype))
374 48 : CALL parser_get_object(parser, itemp)
375 48 : conn_info%theta_type(offset + itype) = itemp
376 62 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLE INFO HERE!"
377 : END DO
378 62 : conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev
379 62 : conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev
380 62 : conn_info%theta_c(offset + 1:offset + ntype) = conn_info%theta_c(offset + 1:offset + ntype) + natom_prev
381 : END IF
382 : ! IMPDIHEDRALH SECTION
383 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRALH section'
384 14 : label = TRIM(avail_section(14))
385 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
386 14 : IF (found) THEN
387 14 : CALL parser_get_next_line(parser, 1)
388 14 : CALL parser_get_object(parser, ntype)
389 14 : offset = 0
390 14 : IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a)
391 14 : CALL reallocate(conn_info%impr_a, 1, offset + ntype)
392 14 : CALL reallocate(conn_info%impr_b, 1, offset + ntype)
393 14 : CALL reallocate(conn_info%impr_c, 1, offset + ntype)
394 14 : CALL reallocate(conn_info%impr_d, 1, offset + ntype)
395 14 : CALL reallocate(conn_info%impr_type, 1, offset + ntype)
396 14 : DO itype = 1, ntype
397 0 : CALL parser_get_next_line(parser, 1)
398 0 : CALL parser_get_object(parser, conn_info%impr_a(offset + itype))
399 0 : CALL parser_get_object(parser, conn_info%impr_b(offset + itype))
400 0 : CALL parser_get_object(parser, conn_info%impr_c(offset + itype))
401 0 : CALL parser_get_object(parser, conn_info%impr_d(offset + itype))
402 0 : CALL parser_get_object(parser, itemp)
403 0 : conn_info%impr_type(offset + itype) = itemp
404 14 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRALH INFO HERE!"
405 : END DO
406 14 : conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev
407 14 : conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev
408 14 : conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev
409 14 : conn_info%impr_d(offset + 1:offset + ntype) = conn_info%impr_d(offset + 1:offset + ntype) + natom_prev
410 : END IF
411 : ! IMPDIHEDRAL SECTION
412 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRAL section'
413 14 : label = TRIM(avail_section(15))
414 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
415 14 : IF (found) THEN
416 14 : CALL parser_get_next_line(parser, 1)
417 14 : CALL parser_get_object(parser, ntype)
418 14 : offset = 0
419 14 : IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a)
420 14 : CALL reallocate(conn_info%impr_a, 1, offset + ntype)
421 14 : CALL reallocate(conn_info%impr_b, 1, offset + ntype)
422 14 : CALL reallocate(conn_info%impr_c, 1, offset + ntype)
423 14 : CALL reallocate(conn_info%impr_d, 1, offset + ntype)
424 14 : CALL reallocate(conn_info%impr_type, 1, offset + ntype)
425 14 : DO itype = 1, ntype
426 0 : CALL parser_get_next_line(parser, 1)
427 0 : CALL parser_get_object(parser, conn_info%impr_a(offset + itype))
428 0 : CALL parser_get_object(parser, conn_info%impr_b(offset + itype))
429 0 : CALL parser_get_object(parser, conn_info%impr_c(offset + itype))
430 0 : CALL parser_get_object(parser, conn_info%impr_d(offset + itype))
431 0 : CALL parser_get_object(parser, itemp)
432 0 : conn_info%impr_type(offset + itype) = itemp
433 14 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRAL INFO HERE!"
434 : END DO
435 14 : conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev
436 14 : conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev
437 14 : conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev
438 14 : conn_info%impr_d(offset + 1:offset + ntype) = conn_info%impr_d(offset + 1:offset + ntype) + natom_prev
439 : END IF
440 : ! DIHEDRALH SECTION
441 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRALH section'
442 14 : label = TRIM(avail_section(17))
443 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
444 14 : IF (found) THEN
445 14 : CALL parser_get_next_line(parser, 1)
446 14 : CALL parser_get_object(parser, ntype)
447 14 : offset = 0
448 14 : IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a)
449 14 : CALL reallocate(conn_info%phi_a, 1, offset + ntype)
450 14 : CALL reallocate(conn_info%phi_b, 1, offset + ntype)
451 14 : CALL reallocate(conn_info%phi_c, 1, offset + ntype)
452 14 : CALL reallocate(conn_info%phi_d, 1, offset + ntype)
453 14 : CALL reallocate(conn_info%phi_type, 1, offset + ntype)
454 270 : DO itype = 1, ntype
455 256 : CALL parser_get_next_line(parser, 1)
456 256 : CALL parser_get_object(parser, conn_info%phi_a(offset + itype))
457 256 : CALL parser_get_object(parser, conn_info%phi_b(offset + itype))
458 256 : CALL parser_get_object(parser, conn_info%phi_c(offset + itype))
459 256 : CALL parser_get_object(parser, conn_info%phi_d(offset + itype))
460 256 : CALL parser_get_object(parser, itemp)
461 256 : conn_info%phi_type(offset + itype) = itemp
462 270 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRALH INFO HERE!"
463 : END DO
464 270 : conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev
465 270 : conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev
466 270 : conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev
467 270 : conn_info%phi_d(offset + 1:offset + ntype) = conn_info%phi_d(offset + 1:offset + ntype) + natom_prev
468 : END IF
469 : ! DIHEDRAL SECTION
470 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRAL section'
471 14 : label = TRIM(avail_section(18))
472 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
473 14 : IF (found) THEN
474 14 : CALL parser_get_next_line(parser, 1)
475 14 : CALL parser_get_object(parser, ntype)
476 14 : offset = 0
477 14 : IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a)
478 14 : CALL reallocate(conn_info%phi_a, 1, offset + ntype)
479 14 : CALL reallocate(conn_info%phi_b, 1, offset + ntype)
480 14 : CALL reallocate(conn_info%phi_c, 1, offset + ntype)
481 14 : CALL reallocate(conn_info%phi_d, 1, offset + ntype)
482 14 : CALL reallocate(conn_info%phi_type, 1, offset + ntype)
483 70 : DO itype = 1, ntype
484 56 : CALL parser_get_next_line(parser, 1)
485 56 : CALL parser_get_object(parser, conn_info%phi_a(offset + itype))
486 56 : CALL parser_get_object(parser, conn_info%phi_b(offset + itype))
487 56 : CALL parser_get_object(parser, conn_info%phi_c(offset + itype))
488 56 : CALL parser_get_object(parser, conn_info%phi_d(offset + itype))
489 56 : CALL parser_get_object(parser, itemp)
490 56 : conn_info%phi_type(offset + itype) = itemp
491 70 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRAL INFO HERE!"
492 : END DO
493 70 : conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev
494 70 : conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev
495 70 : conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev
496 70 : conn_info%phi_d(offset + 1:offset + ntype) = conn_info%phi_d(offset + 1:offset + ntype) + natom_prev
497 : END IF
498 14 : CALL parser_release(parser)
499 :
500 : ! SOLVENTATOM and SOLVENTCONSTR SECTION
501 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLVENTATOM section'
502 14 : nsolvent = (SIZE(atom_info%r(1, :)) - nsolute)/3
503 :
504 14 : NULLIFY (na, am, ac, ba, bb)
505 14 : CALL parser_create(parser, file_name, para_env=para_env)
506 14 : label = TRIM(avail_section(20))
507 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
508 14 : IF (found) THEN
509 14 : CALL parser_get_next_line(parser, 1)
510 14 : CALL parser_get_object(parser, natom)
511 14 : CALL reallocate(na, 1, natom)
512 14 : CALL reallocate(am, 1, natom)
513 14 : CALL reallocate(ac, 1, natom)
514 38 : DO iatom = 1, natom
515 24 : CALL parser_get_next_line(parser, 1)
516 24 : CALL parser_get_object(parser, itemp)
517 24 : CALL parser_get_object(parser, string)
518 24 : CALL parser_get_object(parser, na(iatom))
519 24 : CALL parser_get_object(parser, am(iatom))
520 24 : CALL parser_get_object(parser, ac(iatom))
521 62 : IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLVENTATOM INFO HERE!"
522 : END DO
523 : END IF
524 14 : label = TRIM(avail_section(21))
525 14 : ncon = 0
526 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
527 14 : IF (found) THEN
528 8 : CALL parser_get_next_line(parser, 1)
529 8 : CALL parser_get_object(parser, ncon)
530 8 : CALL reallocate(ba, 1, ncon)
531 8 : CALL reallocate(bb, 1, ncon)
532 32 : DO icon = 1, ncon
533 24 : CALL parser_get_next_line(parser, 1)
534 24 : CALL parser_get_object(parser, ba(icon))
535 32 : CALL parser_get_object(parser, bb(icon))
536 : END DO
537 : END IF
538 14 : CALL parser_release(parser)
539 :
540 14 : offset = 0
541 14 : IF (ASSOCIATED(atom_info%id_molname)) offset = SIZE(atom_info%id_molname)
542 14 : CALL reallocate(atom_info%id_molname, 1, offset + nsolvent*natom)
543 14 : CALL reallocate(atom_info%resid, 1, offset + nsolvent*natom)
544 14 : CALL reallocate(atom_info%id_resname, 1, offset + nsolvent*natom)
545 14 : CALL reallocate(atom_info%id_atmname, 1, offset + nsolvent*natom)
546 14 : CALL reallocate(atom_info%id_element, 1, offset + nsolvent*natom)
547 14 : CALL reallocate(atom_info%atm_charge, 1, offset + nsolvent*natom)
548 14 : CALL reallocate(atom_info%atm_mass, 1, offset + nsolvent*natom)
549 870 : DO isolvent = 1, nsolvent
550 856 : offset = nsolute + natom*isolvent - natom
551 3438 : DO iatom = 1, natom
552 2568 : index_now = offset + iatom
553 2568 : atom_info%id_atmname(index_now) = str2id(s2s(namearray1(na(iatom))))
554 2568 : atom_info%id_element(index_now) = str2id(s2s(namearray1(na(iatom))))
555 2568 : atom_info%id_molname(index_now) = str2id(s2s("SOL"))
556 2568 : atom_info%id_resname(index_now) = str2id(s2s("SOL"))
557 2568 : atom_info%resid(index_now) = isolvent
558 2568 : atom_info%atm_mass(index_now) = am(iatom)
559 3424 : atom_info%atm_charge(index_now) = ac(iatom)
560 : END DO
561 : END DO
562 :
563 14 : offset = 0
564 14 : IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
565 122 : offset2 = MAXVAL(conn_info%bond_type(:))
566 14 : CALL reallocate(conn_info%bond_a, 1, offset + ncon*nsolvent)
567 14 : CALL reallocate(conn_info%bond_b, 1, offset + ncon*nsolvent)
568 14 : CALL reallocate(conn_info%bond_type, 1, offset + ncon*nsolvent)
569 14 : offset = offset - ncon
570 870 : DO isolvent = 1, nsolvent
571 856 : offset = offset + ncon
572 3438 : DO icon = 1, ncon
573 2568 : conn_info%bond_a(offset + icon) = nsolute + isolvent*ncon - ncon + ba(icon)
574 2568 : conn_info%bond_b(offset + icon) = nsolute + isolvent*ncon - ncon + bb(icon)
575 3424 : conn_info%bond_type(offset + icon) = offset2 + isolvent*ncon - ncon + icon
576 : END DO
577 : END DO
578 : ! PARA_RES structure
579 14 : i = 0
580 14 : nbond_prev = 0
581 14 : IF (ASSOCIATED(conn_info%c_bond_a)) i = SIZE(conn_info%c_bond_a)
582 14 : nbond = SIZE(conn_info%bond_a)
583 2690 : DO ibond = 1 + nbond_prev, nbond + nbond_prev
584 2676 : iatom = conn_info%bond_a(ibond)
585 2676 : jatom = conn_info%bond_b(ibond)
586 2690 : IF (topology%para_res) THEN
587 : IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
588 2676 : (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
589 : (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN
590 0 : IF (iw > 0) WRITE (iw, '(T2,A,2I3)') "GTOP_INFO| PARA_RES, bond between molecules atom ", &
591 0 : iatom, jatom
592 0 : i = i + 1
593 0 : CALL reallocate(conn_info%c_bond_a, 1, i)
594 0 : CALL reallocate(conn_info%c_bond_b, 1, i)
595 0 : CALL reallocate(conn_info%c_bond_type, 1, i)
596 0 : conn_info%c_bond_a(i) = iatom
597 0 : conn_info%c_bond_b(i) = jatom
598 0 : conn_info%c_bond_type(i) = conn_info%bond_type(ibond)
599 : END IF
600 : ELSE
601 0 : IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN
602 0 : CPABORT("")
603 : END IF
604 : END IF
605 : END DO
606 :
607 14 : DEALLOCATE (namearray1)
608 14 : DEALLOCATE (namearray2)
609 14 : DEALLOCATE (na)
610 14 : DEALLOCATE (am)
611 14 : DEALLOCATE (ac)
612 14 : IF (ASSOCIATED(ba)) &
613 8 : DEALLOCATE (ba)
614 14 : IF (ASSOCIATED(bb)) &
615 8 : DEALLOCATE (bb)
616 :
617 14 : CALL timestop(handle)
618 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
619 14 : "PRINT%TOPOLOGY_INFO/GTOP_INFO")
620 :
621 42 : END SUBROUTINE read_topology_gromos
622 :
623 : ! **************************************************************************************************
624 : !> \brief ...
625 : !> \param topology ...
626 : !> \param para_env ...
627 : !> \param subsys_section ...
628 : ! **************************************************************************************************
629 56 : SUBROUTINE read_coordinate_g96(topology, para_env, subsys_section)
630 : TYPE(topology_parameters_type) :: topology
631 : TYPE(mp_para_env_type), POINTER :: para_env
632 : TYPE(section_vals_type), POINTER :: subsys_section
633 :
634 : CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_g96'
635 : INTEGER, PARAMETER :: nblock = 1000
636 :
637 : CHARACTER(LEN=default_string_length) :: label, string, strtmp, strtmp2
638 : CHARACTER(LEN=default_string_length), DIMENSION(5) :: avail_section
639 : INTEGER :: handle, itemp, iw, natom, newsize
640 : LOGICAL :: found
641 : REAL(KIND=dp) :: pfactor
642 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: velocity
643 : TYPE(atom_info_type), POINTER :: atom_info
644 : TYPE(cp_logger_type), POINTER :: logger
645 : TYPE(cp_parser_type) :: parser
646 : TYPE(section_vals_type), POINTER :: velocity_section
647 :
648 14 : NULLIFY (logger, velocity)
649 28 : logger => cp_get_default_logger()
650 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/G96_INFO", &
651 14 : extension=".subsysLog")
652 14 : CALL timeset(routineN, handle)
653 :
654 14 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
655 14 : atom_info => topology%atom_info
656 14 : IF (iw > 0) WRITE (iw, *) " Reading in G96 file ", TRIM(topology%coord_file_name)
657 14 : avail_section(1) = "TITLE"
658 14 : avail_section(2) = "TIMESTEP"
659 14 : avail_section(3) = "POSITION"
660 14 : avail_section(4) = "VELOCITY"
661 14 : avail_section(5) = "BOX"
662 :
663 : ! Element is assigned on the basis of the atm_name
664 14 : topology%aa_element = .TRUE.
665 :
666 14 : CALL reallocate(atom_info%id_molname, 1, nblock)
667 14 : CALL reallocate(atom_info%id_resname, 1, nblock)
668 14 : CALL reallocate(atom_info%resid, 1, nblock)
669 14 : CALL reallocate(atom_info%id_atmname, 1, nblock)
670 14 : CALL reallocate(atom_info%id_element, 1, nblock)
671 14 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
672 14 : CALL reallocate(atom_info%atm_mass, 1, nblock)
673 14 : CALL reallocate(atom_info%atm_charge, 1, nblock)
674 14 : CALL reallocate(atom_info%occup, 1, nblock)
675 14 : CALL reallocate(atom_info%beta, 1, nblock)
676 : ! TITLE SECTION
677 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the TITLE section'
678 14 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
679 14 : label = TRIM(avail_section(1))
680 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
681 14 : IF (found) THEN
682 : DO
683 28 : CALL parser_get_next_line(parser, 1)
684 28 : CALL parser_get_object(parser, string, string_length=default_string_length)
685 28 : IF (string == TRIM("END")) EXIT
686 28 : IF (iw > 0) WRITE (iw, *) "G96_INFO| ", TRIM(string)
687 : END DO
688 : END IF
689 14 : CALL parser_release(parser)
690 : ! POSITION SECTION
691 14 : IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the POSITION section'
692 14 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
693 14 : label = TRIM(avail_section(3))
694 14 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
695 14 : IF (found) THEN
696 14 : natom = 0
697 14 : CALL parser_get_next_line(parser, 1)
698 14 : CALL parser_get_object(parser, string, string_length=default_string_length)
699 2690 : DO
700 2704 : IF (string == TRIM("END")) EXIT
701 2690 : natom = natom + 1
702 2690 : IF (natom > SIZE(atom_info%id_molname)) THEN
703 0 : newsize = INT(pfactor*natom)
704 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
705 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
706 0 : CALL reallocate(atom_info%resid, 1, newsize)
707 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
708 0 : CALL reallocate(atom_info%id_element, 1, newsize)
709 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
710 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
711 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
712 0 : CALL reallocate(atom_info%occup, 1, newsize)
713 0 : CALL reallocate(atom_info%beta, 1, newsize)
714 : END IF
715 : READ (string, *) &
716 2690 : atom_info%resid(natom), strtmp, strtmp2, &
717 5380 : itemp, atom_info%r(1, natom), atom_info%r(2, natom), atom_info%r(3, natom)
718 2690 : atom_info%id_resname(natom) = str2id(s2s(strtmp))
719 2690 : atom_info%id_atmname(natom) = str2id(s2s(strtmp2))
720 2690 : atom_info%id_molname(natom) = atom_info%id_resname(natom)
721 2690 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
722 2690 : atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "nm")
723 2690 : atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "nm")
724 2690 : atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "nm")
725 2690 : IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT POSITION INFO HERE!"
726 2690 : CALL parser_get_next_line(parser, 1)
727 2704 : CALL parser_get_object(parser, string, string_length=default_string_length)
728 : END DO
729 : END IF
730 14 : CALL parser_release(parser)
731 14 : CALL reallocate(velocity, 1, 3, 1, natom)
732 : ! VELOCITY SECTION
733 14 : IF (topology%use_g96_velocity) THEN
734 2 : IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the VELOCITY section'
735 2 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
736 2 : label = TRIM(avail_section(4))
737 2 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
738 2 : IF (found) THEN
739 2 : natom = 0
740 2 : CALL parser_get_next_line(parser, 1)
741 2 : CALL parser_get_object(parser, string, string_length=default_string_length)
742 6 : DO
743 8 : IF (string == TRIM("END")) EXIT
744 6 : natom = natom + 1
745 : READ (string, *) &
746 6 : atom_info%resid(natom), strtmp, strtmp2, &
747 12 : itemp, velocity(1, natom), velocity(2, natom), velocity(3, natom)
748 6 : atom_info%id_resname(natom) = str2id(strtmp)
749 6 : atom_info%id_atmname(natom) = str2id(strtmp2)
750 6 : atom_info%id_molname(natom) = atom_info%id_resname(natom)
751 6 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
752 6 : velocity(1, natom) = cp_unit_to_cp2k(velocity(1, natom), "nm*ps^-1")
753 6 : velocity(2, natom) = cp_unit_to_cp2k(velocity(2, natom), "nm*ps^-1")
754 6 : velocity(3, natom) = cp_unit_to_cp2k(velocity(3, natom), "nm*ps^-1")
755 6 : IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT VELOCITY INFO HERE!"
756 6 : CALL parser_get_next_line(parser, 1)
757 8 : CALL parser_get_object(parser, string, string_length=default_string_length)
758 : END DO
759 2 : CALL parser_release(parser)
760 2 : velocity_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
761 : CALL section_velocity_val_set(velocity_section, velocity=velocity, &
762 2 : conv_factor=1.0_dp)
763 : END IF
764 : END IF
765 14 : DEALLOCATE (velocity)
766 :
767 14 : CALL reallocate(atom_info%id_molname, 1, natom)
768 14 : CALL reallocate(atom_info%id_resname, 1, natom)
769 14 : CALL reallocate(atom_info%resid, 1, natom)
770 14 : CALL reallocate(atom_info%id_atmname, 1, natom)
771 14 : CALL reallocate(atom_info%id_element, 1, natom)
772 14 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
773 14 : CALL reallocate(atom_info%atm_mass, 1, natom)
774 14 : CALL reallocate(atom_info%atm_charge, 1, natom)
775 14 : CALL reallocate(atom_info%occup, 1, natom)
776 14 : CALL reallocate(atom_info%beta, 1, natom)
777 :
778 14 : IF (.NOT. topology%para_res) atom_info%resid(:) = 1
779 :
780 14 : topology%natoms = natom
781 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
782 14 : "PRINT%TOPOLOGY_INFO/G96_INFO")
783 14 : CALL timestop(handle)
784 :
785 42 : END SUBROUTINE read_coordinate_g96
786 :
787 : END MODULE topology_gromos
788 :
|
