(git:e7e05ae)
topology_gromos.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
11  cp_logger_type
15  parser_get_object,&
17  USE cp_parser_types, ONLY: cp_parser_type,&
20  USE cp_units, ONLY: cp_unit_to_cp2k
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,&
34 #include "./base/base_uses.f90"
35 
36  IMPLICIT NONE
37 
38  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_gromos'
39 
40  PRIVATE
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  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  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  INTEGER, POINTER :: ba(:), bb(:), na(:)
69  LOGICAL :: found
70  REAL(kind=dp) :: ftemp
71  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  NULLIFY (logger)
78  NULLIFY (namearray1, namearray2)
79  logger => cp_get_default_logger()
80  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GTOP_INFO", &
81  extension=".subsysLog")
82  CALL timeset(routinen, handle)
83 
84  avail_section(1) = "TITLE"
85  avail_section(2) = "TOPPHYSCON"
86  avail_section(3) = "TOPVERSION"
87  avail_section(4) = "ATOMTYPENAME"
88  avail_section(5) = "RESNAME"
89  avail_section(6) = "SOLUTEATOM"
90  avail_section(7) = "BONDTYPE"
91  avail_section(8) = "BONDH"
92  avail_section(9) = "BOND"
93  avail_section(10) = "BONDANGLETYPE"
94  avail_section(11) = "BONDANGLEH"
95  avail_section(12) = "BONDANGLE"
96  avail_section(13) = "IMPDIHEDRALTYPE"
97  avail_section(14) = "IMPDIHEDRALH"
98  avail_section(15) = "IMPDIHEDRAL"
99  avail_section(16) = "DIHEDRALTYPE"
100  avail_section(17) = "DIHEDRALH"
101  avail_section(18) = "DIHEDRAL"
102  avail_section(19) = "LJPARAMETERS"
103  avail_section(20) = "SOLVENTATOM"
104  avail_section(21) = "SOLVENTCONSTR"
105 
106  atom_info => topology%atom_info
107  conn_info => topology%conn_info
108 
109  natom_prev = 0
110  IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
111  ! TITLE SECTION
112  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TITLE section'
113  CALL parser_create(parser, file_name, para_env=para_env)
114  label = trim(avail_section(1))
115  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
116  IF (found) THEN
117  DO
118  CALL parser_get_next_line(parser, 1)
119  CALL parser_get_object(parser, string, string_length=80)
120  IF (string == trim("END")) EXIT
121  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(string)
122  END DO
123  END IF
124  CALL parser_release(parser)
125  ! TOPPHYSCON SECTION
126  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPPHYSCON section'
127  CALL parser_create(parser, file_name, para_env=para_env)
128  label = trim(avail_section(2))
129  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
130  IF (found) THEN
131  DO
132  CALL parser_get_next_line(parser, 1)
133  CALL parser_get_object(parser, string)
134  IF (string == trim("END")) EXIT
135  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(string)
136  END DO
137  END IF
138  CALL parser_release(parser)
139  ! TOPVERSION SECTION
140  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPVERSION section'
141  CALL parser_create(parser, file_name, para_env=para_env)
142  label = trim(avail_section(3))
143  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
144  IF (found) THEN
145  DO
146  CALL parser_get_next_line(parser, 1)
147  CALL parser_get_object(parser, string)
148  IF (string == trim("END")) EXIT
149  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(string)
150  END DO
151  END IF
152  CALL parser_release(parser)
153  ! ATOMTYPENAME SECTION
154  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section'
155  CALL parser_create(parser, file_name, para_env=para_env)
156  label = trim(avail_section(4))
157  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
158  IF (found) THEN
159  CALL parser_get_next_line(parser, 1)
160  CALL parser_get_object(parser, ntype)
161  CALL reallocate(namearray1, 1, ntype)
162  DO itype = 1, ntype
163  CALL parser_get_next_line(parser, 1)
164  CALL parser_get_object(parser, namearray1(itype))
165  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(namearray1(itype))
166  END DO
167  END IF
168  CALL parser_release(parser)
169  ! RESNAME SECTION
170  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the RESNAME section'
171  CALL parser_create(parser, file_name, para_env=para_env)
172  label = trim(avail_section(5))
173  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
174  IF (found) THEN
175  CALL parser_get_next_line(parser, 1)
176  CALL parser_get_object(parser, ntype)
177  CALL reallocate(namearray2, 1, ntype)
178  DO itype = 1, ntype
179  CALL parser_get_next_line(parser, 1)
180  CALL parser_get_object(parser, namearray2(itype))
181  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", trim(namearray2(itype))
182  END DO
183  END IF
184  CALL parser_release(parser)
185  ! SOLUTEATOM SECTION
186  iresid = 1
187  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLUTEATOM section'
188  CALL parser_create(parser, file_name, para_env=para_env)
189  label = trim(avail_section(6))
190  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
191  IF (found) THEN
192  CALL parser_get_next_line(parser, 1)
193  CALL parser_get_object(parser, natom)
194  CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
195  CALL reallocate(atom_info%resid, 1, natom_prev + natom)
196  CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
197  CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
198  CALL reallocate(atom_info%id_element, 1, natom_prev + natom)
199  CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
200  CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
201  CALL parser_get_next_line(parser, 1)
202  DO iatom = 1, natom
203  index_now = iatom + natom_prev
204  CALL parser_get_object(parser, itemp)
205  CALL parser_get_object(parser, itemp)
206  atom_info%resid(index_now) = itemp
207  atom_info%id_molname(index_now) = str2id(s2s(namearray2(itemp)))
208  atom_info%id_resname(index_now) = str2id(s2s(namearray2(itemp)))
209  CALL parser_get_object(parser, string)
210  CALL parser_get_object(parser, itemp)
211  atom_info%id_atmname(index_now) = str2id(s2s(namearray1(itemp)))
212  atom_info%id_element(index_now) = str2id(s2s(namearray1(itemp)))
213  CALL parser_get_object(parser, atom_info%atm_mass(index_now))
214  CALL parser_get_object(parser, atom_info%atm_charge(index_now))
215  CALL parser_get_object(parser, itemp)
216  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLUTEATOM INFO HERE!!!!"
217  CALL parser_get_object(parser, ntype)
218  DO i = 1, 50
219  ii(i) = -1
220  END DO
221  stat = -1
222  begin = 1
223  IF (ntype /= 0) THEN
224  DO
225  IF (begin .EQ. 1) THEN
226  READ (parser%input_line, iostat=stat, fmt=*) itemp, itemp, ctemp, itemp, ftemp, ftemp, &
227  itemp, itemp, (ii(i), i=begin, ntype)
228  ELSE IF (begin .GT. 1) THEN
229  CALL parser_get_next_line(parser, 1)
230  READ (parser%input_line, iostat=stat, fmt=*) (ii(i), i=begin, ntype)
231  END IF
232  DO i = ntype, 1, -1
233  IF (ii(i) .LT. 0) THEN
234  begin = i
235  END IF
236  END DO
237  IF (stat .EQ. 0) EXIT
238  END DO
239  END IF
240  ! 1-4 list
241  CALL parser_get_next_line(parser, 1)
242  CALL parser_get_object(parser, ntype)
243  IF (ntype /= 0) THEN
244  itemp = (itemp - 1)/6 + 1
245  offset = 0
246  IF (ASSOCIATED(conn_info%onfo_a)) offset = SIZE(conn_info%onfo_a)
247  CALL reallocate(conn_info%onfo_a, 1, offset + ntype)
248  CALL reallocate(conn_info%onfo_b, 1, offset + ntype)
249  conn_info%onfo_a(offset + 1:offset + ntype) = index_now
250  DO i = 1, 50
251  ii(i) = -1
252  END DO
253  stat = -1
254  begin = 1
255  IF (ntype /= 0) THEN
256  DO
257  IF (begin .EQ. 1) THEN
258  READ (parser%input_line, iostat=stat, fmt=*) itemp, (ii(i), i=begin, ntype)
259  ELSE IF (begin .GT. 1) THEN
260  CALL parser_get_next_line(parser, 1)
261  READ (parser%input_line, iostat=stat, fmt=*) (ii(i), i=begin, ntype)
262  END IF
263  DO i = ntype, 1, -1
264  IF (ii(i) .LT. 0) THEN
265  begin = i
266  END IF
267  END DO
268  IF (stat .EQ. 0) EXIT
269  END DO
270  DO i = 1, ntype
271  conn_info%onfo_b(offset + i) = ii(i)
272  END DO
273  END IF
274  CALL parser_get_next_line(parser, 1)
275  ELSE
276  CALL parser_get_next_line(parser, 1)
277  END IF
278  END DO
279  END IF
280  nsolute = natom
281  CALL parser_release(parser)
282 
283  CALL parser_create(parser, file_name, para_env=para_env)
284  ! BONDH SECTION
285  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDH section'
286  label = trim(avail_section(8))
287  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
288  IF (found) THEN
289  CALL parser_get_next_line(parser, 1)
290  CALL parser_get_object(parser, ntype)
291  offset = 0
292  IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
293  CALL reallocate(conn_info%bond_a, 1, offset + ntype)
294  CALL reallocate(conn_info%bond_b, 1, offset + ntype)
295  CALL reallocate(conn_info%bond_type, 1, offset + ntype)
296  DO itype = 1, ntype
297  CALL parser_get_next_line(parser, 1)
298  CALL parser_get_object(parser, conn_info%bond_a(offset + itype))
299  CALL parser_get_object(parser, conn_info%bond_b(offset + itype))
300  CALL parser_get_object(parser, itemp)
301  conn_info%bond_type(offset + itype) = itemp
302  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDH INFO HERE!!!!"
303  END DO
304  conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev
305  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BOND section'
309  label = trim(avail_section(9))
310  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
311  IF (found) THEN
312  CALL parser_get_next_line(parser, 1)
313  CALL parser_get_object(parser, ntype)
314  offset = 0
315  IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
316  CALL reallocate(conn_info%bond_a, 1, offset + ntype)
317  CALL reallocate(conn_info%bond_b, 1, offset + ntype)
318  CALL reallocate(conn_info%bond_type, 1, offset + ntype)
319  DO itype = 1, ntype
320  CALL parser_get_next_line(parser, 1)
321  CALL parser_get_object(parser, conn_info%bond_a(offset + itype))
322  CALL parser_get_object(parser, conn_info%bond_b(offset + itype))
323  CALL parser_get_object(parser, itemp)
324  conn_info%bond_type(offset + itype) = itemp
325  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BOND INFO HERE!!!!"
326  END DO
327  conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev
328  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLEH section'
332  label = trim(avail_section(11))
333  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
334  IF (found) THEN
335  CALL parser_get_next_line(parser, 1)
336  CALL parser_get_object(parser, ntype)
337  offset = 0
338  IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a)
339  CALL reallocate(conn_info%theta_a, 1, offset + ntype)
340  CALL reallocate(conn_info%theta_b, 1, offset + ntype)
341  CALL reallocate(conn_info%theta_c, 1, offset + ntype)
342  CALL reallocate(conn_info%theta_type, 1, offset + ntype)
343  DO itype = 1, ntype
344  CALL parser_get_next_line(parser, 1)
345  CALL parser_get_object(parser, conn_info%theta_a(offset + itype))
346  CALL parser_get_object(parser, conn_info%theta_b(offset + itype))
347  CALL parser_get_object(parser, conn_info%theta_c(offset + itype))
348  CALL parser_get_object(parser, itemp)
349  conn_info%theta_type(offset + itype) = itemp
350  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLEH INFO HERE!!!!"
351  END DO
352  conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev
353  conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev
354  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLE section'
358  label = trim(avail_section(12))
359  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
360  IF (found) THEN
361  CALL parser_get_next_line(parser, 1)
362  CALL parser_get_object(parser, ntype)
363  offset = 0
364  IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a)
365  CALL reallocate(conn_info%theta_a, 1, offset + ntype)
366  CALL reallocate(conn_info%theta_b, 1, offset + ntype)
367  CALL reallocate(conn_info%theta_c, 1, offset + ntype)
368  CALL reallocate(conn_info%theta_type, 1, offset + ntype)
369  DO itype = 1, ntype
370  CALL parser_get_next_line(parser, 1)
371  CALL parser_get_object(parser, conn_info%theta_a(offset + itype))
372  CALL parser_get_object(parser, conn_info%theta_b(offset + itype))
373  CALL parser_get_object(parser, conn_info%theta_c(offset + itype))
374  CALL parser_get_object(parser, itemp)
375  conn_info%theta_type(offset + itype) = itemp
376  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLE INFO HERE!!!!"
377  END DO
378  conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev
379  conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev
380  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRALH section'
384  label = trim(avail_section(14))
385  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
386  IF (found) THEN
387  CALL parser_get_next_line(parser, 1)
388  CALL parser_get_object(parser, ntype)
389  offset = 0
390  IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a)
391  CALL reallocate(conn_info%impr_a, 1, offset + ntype)
392  CALL reallocate(conn_info%impr_b, 1, offset + ntype)
393  CALL reallocate(conn_info%impr_c, 1, offset + ntype)
394  CALL reallocate(conn_info%impr_d, 1, offset + ntype)
395  CALL reallocate(conn_info%impr_type, 1, offset + ntype)
396  DO itype = 1, ntype
397  CALL parser_get_next_line(parser, 1)
398  CALL parser_get_object(parser, conn_info%impr_a(offset + itype))
399  CALL parser_get_object(parser, conn_info%impr_b(offset + itype))
400  CALL parser_get_object(parser, conn_info%impr_c(offset + itype))
401  CALL parser_get_object(parser, conn_info%impr_d(offset + itype))
402  CALL parser_get_object(parser, itemp)
403  conn_info%impr_type(offset + itype) = itemp
404  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRALH INFO HERE!!!!"
405  END DO
406  conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev
407  conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev
408  conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev
409  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRAL section'
413  label = trim(avail_section(15))
414  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
415  IF (found) THEN
416  CALL parser_get_next_line(parser, 1)
417  CALL parser_get_object(parser, ntype)
418  offset = 0
419  IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a)
420  CALL reallocate(conn_info%impr_a, 1, offset + ntype)
421  CALL reallocate(conn_info%impr_b, 1, offset + ntype)
422  CALL reallocate(conn_info%impr_c, 1, offset + ntype)
423  CALL reallocate(conn_info%impr_d, 1, offset + ntype)
424  CALL reallocate(conn_info%impr_type, 1, offset + ntype)
425  DO itype = 1, ntype
426  CALL parser_get_next_line(parser, 1)
427  CALL parser_get_object(parser, conn_info%impr_a(offset + itype))
428  CALL parser_get_object(parser, conn_info%impr_b(offset + itype))
429  CALL parser_get_object(parser, conn_info%impr_c(offset + itype))
430  CALL parser_get_object(parser, conn_info%impr_d(offset + itype))
431  CALL parser_get_object(parser, itemp)
432  conn_info%impr_type(offset + itype) = itemp
433  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRAL INFO HERE!!!!"
434  END DO
435  conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev
436  conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev
437  conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev
438  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRALH section'
442  label = trim(avail_section(17))
443  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
444  IF (found) THEN
445  CALL parser_get_next_line(parser, 1)
446  CALL parser_get_object(parser, ntype)
447  offset = 0
448  IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a)
449  CALL reallocate(conn_info%phi_a, 1, offset + ntype)
450  CALL reallocate(conn_info%phi_b, 1, offset + ntype)
451  CALL reallocate(conn_info%phi_c, 1, offset + ntype)
452  CALL reallocate(conn_info%phi_d, 1, offset + ntype)
453  CALL reallocate(conn_info%phi_type, 1, offset + ntype)
454  DO itype = 1, ntype
455  CALL parser_get_next_line(parser, 1)
456  CALL parser_get_object(parser, conn_info%phi_a(offset + itype))
457  CALL parser_get_object(parser, conn_info%phi_b(offset + itype))
458  CALL parser_get_object(parser, conn_info%phi_c(offset + itype))
459  CALL parser_get_object(parser, conn_info%phi_d(offset + itype))
460  CALL parser_get_object(parser, itemp)
461  conn_info%phi_type(offset + itype) = itemp
462  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRALH INFO HERE!!!!"
463  END DO
464  conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev
465  conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev
466  conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev
467  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  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRAL section'
471  label = trim(avail_section(18))
472  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
473  IF (found) THEN
474  CALL parser_get_next_line(parser, 1)
475  CALL parser_get_object(parser, ntype)
476  offset = 0
477  IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a)
478  CALL reallocate(conn_info%phi_a, 1, offset + ntype)
479  CALL reallocate(conn_info%phi_b, 1, offset + ntype)
480  CALL reallocate(conn_info%phi_c, 1, offset + ntype)
481  CALL reallocate(conn_info%phi_d, 1, offset + ntype)
482  CALL reallocate(conn_info%phi_type, 1, offset + ntype)
483  DO itype = 1, ntype
484  CALL parser_get_next_line(parser, 1)
485  CALL parser_get_object(parser, conn_info%phi_a(offset + itype))
486  CALL parser_get_object(parser, conn_info%phi_b(offset + itype))
487  CALL parser_get_object(parser, conn_info%phi_c(offset + itype))
488  CALL parser_get_object(parser, conn_info%phi_d(offset + itype))
489  CALL parser_get_object(parser, itemp)
490  conn_info%phi_type(offset + itype) = itemp
491  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRAL INFO HERE!!!!"
492  END DO
493  conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev
494  conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev
495  conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev
496  conn_info%phi_d(offset + 1:offset + ntype) = conn_info%phi_d(offset + 1:offset + ntype) + natom_prev
497  END IF
498  CALL parser_release(parser)
499 
500  ! SOLVENTATOM and SOLVENTCONSTR SECTION
501  IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLVENTATOM section'
502  nsolvent = (SIZE(atom_info%r(1, :)) - nsolute)/3
503 
504  NULLIFY (na, am, ac, ba, bb)
505  CALL parser_create(parser, file_name, para_env=para_env)
506  label = trim(avail_section(20))
507  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
508  IF (found) THEN
509  CALL parser_get_next_line(parser, 1)
510  CALL parser_get_object(parser, natom)
511  CALL reallocate(na, 1, natom)
512  CALL reallocate(am, 1, natom)
513  CALL reallocate(ac, 1, natom)
514  DO iatom = 1, natom
515  CALL parser_get_next_line(parser, 1)
516  CALL parser_get_object(parser, itemp)
517  CALL parser_get_object(parser, string)
518  CALL parser_get_object(parser, na(iatom))
519  CALL parser_get_object(parser, am(iatom))
520  CALL parser_get_object(parser, ac(iatom))
521  IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLVENTATOM INFO HERE!!!!"
522  END DO
523  END IF
524  label = trim(avail_section(21))
525  ncon = 0
526  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
527  IF (found) THEN
528  CALL parser_get_next_line(parser, 1)
529  CALL parser_get_object(parser, ncon)
530  CALL reallocate(ba, 1, ncon)
531  CALL reallocate(bb, 1, ncon)
532  DO icon = 1, ncon
533  CALL parser_get_next_line(parser, 1)
534  CALL parser_get_object(parser, ba(icon))
535  CALL parser_get_object(parser, bb(icon))
536  END DO
537  END IF
538  CALL parser_release(parser)
539 
540  offset = 0
541  IF (ASSOCIATED(atom_info%id_molname)) offset = SIZE(atom_info%id_molname)
542  CALL reallocate(atom_info%id_molname, 1, offset + nsolvent*natom)
543  CALL reallocate(atom_info%resid, 1, offset + nsolvent*natom)
544  CALL reallocate(atom_info%id_resname, 1, offset + nsolvent*natom)
545  CALL reallocate(atom_info%id_atmname, 1, offset + nsolvent*natom)
546  CALL reallocate(atom_info%id_element, 1, offset + nsolvent*natom)
547  CALL reallocate(atom_info%atm_charge, 1, offset + nsolvent*natom)
548  CALL reallocate(atom_info%atm_mass, 1, offset + nsolvent*natom)
549  DO isolvent = 1, nsolvent
550  offset = nsolute + natom*isolvent - natom
551  DO iatom = 1, natom
552  index_now = offset + iatom
553  atom_info%id_atmname(index_now) = str2id(s2s(namearray1(na(iatom))))
554  atom_info%id_element(index_now) = str2id(s2s(namearray1(na(iatom))))
555  atom_info%id_molname(index_now) = str2id(s2s("SOL"))
556  atom_info%id_resname(index_now) = str2id(s2s("SOL"))
557  atom_info%resid(index_now) = isolvent
558  atom_info%atm_mass(index_now) = am(iatom)
559  atom_info%atm_charge(index_now) = ac(iatom)
560  END DO
561  END DO
562 
563  offset = 0
564  IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a)
565  offset2 = maxval(conn_info%bond_type(:))
566  CALL reallocate(conn_info%bond_a, 1, offset + ncon*nsolvent)
567  CALL reallocate(conn_info%bond_b, 1, offset + ncon*nsolvent)
568  CALL reallocate(conn_info%bond_type, 1, offset + ncon*nsolvent)
569  offset = offset - ncon
570  DO isolvent = 1, nsolvent
571  offset = offset + ncon
572  DO icon = 1, ncon
573  conn_info%bond_a(offset + icon) = nsolute + isolvent*ncon - ncon + ba(icon)
574  conn_info%bond_b(offset + icon) = nsolute + isolvent*ncon - ncon + bb(icon)
575  conn_info%bond_type(offset + icon) = offset2 + isolvent*ncon - ncon + icon
576  END DO
577  END DO
578  ! PARA_RES structure
579  i = 0
580  nbond_prev = 0
581  IF (ASSOCIATED(conn_info%c_bond_a)) i = SIZE(conn_info%c_bond_a)
582  nbond = SIZE(conn_info%bond_a)
583  DO ibond = 1 + nbond_prev, nbond + nbond_prev
584  iatom = conn_info%bond_a(ibond)
585  jatom = conn_info%bond_b(ibond)
586  IF (topology%para_res) THEN
587  IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
588  (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
589  (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN
590  IF (iw > 0) WRITE (iw, '(T2,A,2I3)') "GTOP_INFO| PARA_RES, bond between molecules atom ", &
591  iatom, jatom
592  i = i + 1
593  CALL reallocate(conn_info%c_bond_a, 1, i)
594  CALL reallocate(conn_info%c_bond_b, 1, i)
595  CALL reallocate(conn_info%c_bond_type, 1, i)
596  conn_info%c_bond_a(i) = iatom
597  conn_info%c_bond_b(i) = jatom
598  conn_info%c_bond_type(i) = conn_info%bond_type(ibond)
599  END IF
600  ELSE
601  IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN
602  cpabort("")
603  END IF
604  END IF
605  END DO
606 
607  DEALLOCATE (namearray1)
608  DEALLOCATE (namearray2)
609  DEALLOCATE (na)
610  DEALLOCATE (am)
611  DEALLOCATE (ac)
612  IF (ASSOCIATED(ba)) &
613  DEALLOCATE (ba)
614  IF (ASSOCIATED(bb)) &
615  DEALLOCATE (bb)
616 
617  CALL timestop(handle)
618  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
619  "PRINT%TOPOLOGY_INFO/GTOP_INFO")
620 
621  END SUBROUTINE read_topology_gromos
622 
623 ! **************************************************************************************************
624 !> \brief ...
625 !> \param topology ...
626 !> \param para_env ...
627 !> \param subsys_section ...
628 ! **************************************************************************************************
629  SUBROUTINE read_coordinate_g96(topology, para_env, subsys_section)
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  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  NULLIFY (logger, velocity)
649  logger => cp_get_default_logger()
650  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/G96_INFO", &
651  extension=".subsysLog")
652  CALL timeset(routinen, handle)
653 
654  pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
655  atom_info => topology%atom_info
656  IF (iw > 0) WRITE (iw, *) " Reading in G96 file ", trim(topology%coord_file_name)
657  avail_section(1) = "TITLE"
658  avail_section(2) = "TIMESTEP"
659  avail_section(3) = "POSITION"
660  avail_section(4) = "VELOCITY"
661  avail_section(5) = "BOX"
662 
663  ! Element is assigned on the basis of the atm_name
664  topology%aa_element = .true.
665 
666  CALL reallocate(atom_info%id_molname, 1, nblock)
667  CALL reallocate(atom_info%id_resname, 1, nblock)
668  CALL reallocate(atom_info%resid, 1, nblock)
669  CALL reallocate(atom_info%id_atmname, 1, nblock)
670  CALL reallocate(atom_info%id_element, 1, nblock)
671  CALL reallocate(atom_info%r, 1, 3, 1, nblock)
672  CALL reallocate(atom_info%atm_mass, 1, nblock)
673  CALL reallocate(atom_info%atm_charge, 1, nblock)
674  CALL reallocate(atom_info%occup, 1, nblock)
675  CALL reallocate(atom_info%beta, 1, nblock)
676  ! TITLE SECTION
677  IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the TITLE section'
678  CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
679  label = trim(avail_section(1))
680  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
681  IF (found) THEN
682  DO
683  CALL parser_get_next_line(parser, 1)
684  CALL parser_get_object(parser, string, string_length=default_string_length)
685  IF (string == trim("END")) EXIT
686  IF (iw > 0) WRITE (iw, *) "G96_INFO| ", trim(string)
687  END DO
688  END IF
689  CALL parser_release(parser)
690  ! POSITION SECTION
691  IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the POSITION section'
692  CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
693  label = trim(avail_section(3))
694  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
695  IF (found) THEN
696  natom = 0
697  CALL parser_get_next_line(parser, 1)
698  CALL parser_get_object(parser, string, string_length=default_string_length)
699  DO
700  IF (string == trim("END")) EXIT
701  natom = natom + 1
702  IF (natom > SIZE(atom_info%id_molname)) THEN
703  newsize = int(pfactor*natom)
704  CALL reallocate(atom_info%id_molname, 1, newsize)
705  CALL reallocate(atom_info%id_resname, 1, newsize)
706  CALL reallocate(atom_info%resid, 1, newsize)
707  CALL reallocate(atom_info%id_atmname, 1, newsize)
708  CALL reallocate(atom_info%id_element, 1, newsize)
709  CALL reallocate(atom_info%r, 1, 3, 1, newsize)
710  CALL reallocate(atom_info%atm_mass, 1, newsize)
711  CALL reallocate(atom_info%atm_charge, 1, newsize)
712  CALL reallocate(atom_info%occup, 1, newsize)
713  CALL reallocate(atom_info%beta, 1, newsize)
714  END IF
715  READ (string, *) &
716  atom_info%resid(natom), strtmp, strtmp2, &
717  itemp, atom_info%r(1, natom), atom_info%r(2, natom), atom_info%r(3, natom)
718  atom_info%id_resname(natom) = str2id(s2s(strtmp))
719  atom_info%id_atmname(natom) = str2id(s2s(strtmp2))
720  atom_info%id_molname(natom) = atom_info%id_resname(natom)
721  atom_info%id_element(natom) = atom_info%id_atmname(natom)
722  atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "nm")
723  atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "nm")
724  atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "nm")
725  IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT POSITION INFO HERE!!!!"
726  CALL parser_get_next_line(parser, 1)
727  CALL parser_get_object(parser, string, string_length=default_string_length)
728  END DO
729  END IF
730  CALL parser_release(parser)
731  CALL reallocate(velocity, 1, 3, 1, natom)
732  ! VELOCITY SECTION
733  IF (topology%use_g96_velocity) THEN
734  IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the VELOCITY section'
735  CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
736  label = trim(avail_section(4))
737  CALL parser_search_string(parser, label, .true., found, begin_line=.true.)
738  IF (found) THEN
739  natom = 0
740  CALL parser_get_next_line(parser, 1)
741  CALL parser_get_object(parser, string, string_length=default_string_length)
742  DO
743  IF (string == trim("END")) EXIT
744  natom = natom + 1
745  READ (string, *) &
746  atom_info%resid(natom), strtmp, strtmp2, &
747  itemp, velocity(1, natom), velocity(2, natom), velocity(3, natom)
748  atom_info%id_resname(natom) = str2id(strtmp)
749  atom_info%id_atmname(natom) = str2id(strtmp2)
750  atom_info%id_molname(natom) = atom_info%id_resname(natom)
751  atom_info%id_element(natom) = atom_info%id_atmname(natom)
752  velocity(1, natom) = cp_unit_to_cp2k(velocity(1, natom), "nm*ps^-1")
753  velocity(2, natom) = cp_unit_to_cp2k(velocity(2, natom), "nm*ps^-1")
754  velocity(3, natom) = cp_unit_to_cp2k(velocity(3, natom), "nm*ps^-1")
755  IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT VELOCITY INFO HERE!!!!"
756  CALL parser_get_next_line(parser, 1)
757  CALL parser_get_object(parser, string, string_length=default_string_length)
758  END DO
759  CALL parser_release(parser)
760  velocity_section => section_vals_get_subs_vals(subsys_section, "VELOCITY")
761  CALL section_velocity_val_set(velocity_section, velocity=velocity, &
762  conv_factor=1.0_dp)
763  END IF
764  END IF
765  DEALLOCATE (velocity)
766 
767  CALL reallocate(atom_info%id_molname, 1, natom)
768  CALL reallocate(atom_info%id_resname, 1, natom)
769  CALL reallocate(atom_info%resid, 1, natom)
770  CALL reallocate(atom_info%id_atmname, 1, natom)
771  CALL reallocate(atom_info%id_element, 1, natom)
772  CALL reallocate(atom_info%r, 1, 3, 1, natom)
773  CALL reallocate(atom_info%atm_mass, 1, natom)
774  CALL reallocate(atom_info%atm_charge, 1, natom)
775  CALL reallocate(atom_info%occup, 1, natom)
776  CALL reallocate(atom_info%beta, 1, natom)
777 
778  IF (.NOT. topology%para_res) atom_info%resid(:) = 1
779 
780  topology%natoms = natom
781  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
782  "PRINT%TOPOLOGY_INFO/G96_INFO")
783  CALL timestop(handle)
784 
785  END SUBROUTINE read_coordinate_g96
786 
787 END MODULE topology_gromos
788 
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition: cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
subroutine, public section_velocity_val_set(velocity_section, particles, velocity, conv_factor)
routine to dump velocities.. fast implementation
objects that represent the structure of input sections and the data contained in an input section
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Utility routines for the memory handling.
Interface to the message passing library MPI.
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition: string_table.F:72
subroutine, public read_coordinate_g96(topology, para_env, subsys_section)
...
subroutine, public read_topology_gromos(file_name, topology, para_env, subsys_section)
Read GROMOS topology file.
Control for reading in different topologies and coordinates.
Definition: topology.F:13