Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Functionality to read in PSF topologies and convert it into local
10 : !> data structures
11 : !> \author ikuo
12 : !> tlaino 10.2006
13 : ! **************************************************************************************************
14 : MODULE topology_psf
15 : USE cp_log_handling, ONLY: cp_get_default_logger,&
16 : cp_logger_get_default_io_unit,&
17 : cp_logger_type,&
18 : cp_to_string
19 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
20 : cp_print_key_generate_filename,&
21 : cp_print_key_unit_nr
22 : USE cp_parser_methods, ONLY: parser_get_next_line,&
23 : parser_get_object,&
24 : parser_search_string,&
25 : parser_test_next_token
26 : USE cp_parser_types, ONLY: cp_parser_type,&
27 : parser_create,&
28 : parser_release
29 : USE force_fields_input, ONLY: read_chrg_section
30 : USE input_constants, ONLY: do_conn_psf,&
31 : do_conn_psf_u
32 : USE input_section_types, ONLY: section_vals_get,&
33 : section_vals_get_subs_vals,&
34 : section_vals_type,&
35 : section_vals_val_get
36 : USE kinds, ONLY: default_path_length,&
37 : default_string_length,&
38 : dp
39 : USE memory_utilities, ONLY: reallocate
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
42 : USE string_table, ONLY: id2str,&
43 : s2s,&
44 : str2id
45 : USE string_utilities, ONLY: uppercase
46 : USE topology_types, ONLY: atom_info_type,&
47 : connectivity_info_type,&
48 : topology_parameters_type
49 : USE topology_util, ONLY: array1_list_type,&
50 : reorder_structure,&
51 : tag_molecule
52 : USE util, ONLY: sort
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_psf'
58 :
59 : PRIVATE
60 : PUBLIC :: read_topology_psf, &
61 : write_topology_psf, &
62 : psf_post_process, &
63 : idm_psf
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Read PSF topology file
69 : !> Teodoro Laino - Introduced CHARMM31 EXT PSF standard format
70 : !> \param filename ...
71 : !> \param topology ...
72 : !> \param para_env ...
73 : !> \param subsys_section ...
74 : !> \param psf_type ...
75 : !> \par History
76 : !> 04-2007 Teodoro Laino - Zurich University [tlaino]
77 : !> This routine should contain only information read from the PSF format
78 : !> and all post_process should be performef in the psf_post_process
79 : ! **************************************************************************************************
80 24960 : SUBROUTINE read_topology_psf(filename, topology, para_env, subsys_section, psf_type)
81 : CHARACTER(LEN=*), INTENT(IN) :: filename
82 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
83 : TYPE(mp_para_env_type), POINTER :: para_env
84 : TYPE(section_vals_type), POINTER :: subsys_section
85 : INTEGER, INTENT(IN) :: psf_type
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_psf'
88 :
89 : CHARACTER(LEN=2*default_string_length) :: psf_format
90 : CHARACTER(LEN=3) :: c_int
91 : CHARACTER(LEN=default_string_length) :: dummy_field, field, label, strtmp1, &
92 : strtmp2, strtmp3
93 : INTEGER :: handle, i, iatom, ibond, idum, index_now, iphi, itheta, iw, natom, natom_prev, &
94 : nbond, nbond_prev, nphi, nphi_prev, ntheta, ntheta_prev, output_unit
95 : LOGICAL :: found
96 : TYPE(atom_info_type), POINTER :: atom_info
97 : TYPE(connectivity_info_type), POINTER :: conn_info
98 : TYPE(cp_logger_type), POINTER :: logger
99 : TYPE(cp_parser_type) :: parser
100 :
101 8320 : NULLIFY (logger)
102 16640 : logger => cp_get_default_logger()
103 8320 : output_unit = cp_logger_get_default_io_unit(logger)
104 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", &
105 8320 : extension=".subsysLog")
106 8320 : CALL timeset(routineN, handle)
107 8320 : CALL parser_create(parser, filename, para_env=para_env)
108 :
109 8320 : atom_info => topology%atom_info
110 8320 : conn_info => topology%conn_info
111 8320 : natom_prev = 0
112 8320 : IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname)
113 8320 : c_int = 'I8'
114 8320 : label = 'PSF'
115 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
116 8320 : IF (.NOT. found) &
117 : CALL cp_abort(__LOCATION__, &
118 0 : "Missing PSF specification line in <"//TRIM(filename)//">")
119 18632 : DO WHILE (parser_test_next_token(parser) /= "EOL")
120 10312 : CALL parser_get_object(parser, field)
121 8320 : SELECT CASE (field(1:3))
122 : CASE ("PSF")
123 8320 : IF (psf_type == do_conn_psf) THEN
124 : ! X-PLOR PSF format "similar" to the plain CHARMM PSF format
125 7974 : psf_format = '(I8,1X,A4,I5,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8)'
126 : END IF
127 : CASE ("EXT")
128 1992 : IF (psf_type == do_conn_psf) THEN
129 : ! EXTEnded CHARMM31 format
130 1992 : psf_format = '(I10,T12,A7,T21,I8,T30,A7,T39,A6,T47,A6,T53,F10.6,T69,F8.3,T88,I1)'
131 1992 : c_int = 'I10'
132 : ELSE
133 0 : CPABORT("PSF_INFO| "//field(1:3)//" :: not available for UPSF format!")
134 : END IF
135 : CASE DEFAULT
136 10312 : CPABORT("PSF_INFO| "//field(1:3)//" :: Unimplemented keyword in CP2K PSF/UPSF format!")
137 : END SELECT
138 : END DO
139 8320 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NATOM section'
140 : !
141 : ! ATOM section
142 : !
143 8320 : label = '!NATOM'
144 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
145 8320 : IF (.NOT. found) THEN
146 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NATOM section '
147 0 : natom = 0
148 : ELSE
149 8320 : CALL parser_get_object(parser, natom)
150 8320 : IF (natom_prev + natom > topology%natoms) &
151 : CALL cp_abort(__LOCATION__, &
152 : "Number of atoms in connectivity control is larger than the "// &
153 : "number of atoms in coordinate control. check coordinates and "// &
154 0 : "connectivity. ")
155 8320 : IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NATOM = ', natom
156 : !malloc the memory that we need
157 8320 : CALL reallocate(atom_info%id_molname, 1, natom_prev + natom)
158 8320 : CALL reallocate(atom_info%resid, 1, natom_prev + natom)
159 8320 : CALL reallocate(atom_info%id_resname, 1, natom_prev + natom)
160 8320 : CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom)
161 8320 : CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom)
162 8320 : CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom)
163 : !Read in the atom info
164 8320 : IF (psf_type == do_conn_psf_u) THEN
165 194588 : DO iatom = 1, natom
166 194242 : index_now = iatom + natom_prev
167 194242 : CALL parser_get_next_line(parser, 1)
168 194242 : READ (parser%input_line, FMT=*, ERR=9) i, &
169 194242 : strtmp1, &
170 194242 : atom_info%resid(index_now), &
171 194242 : strtmp2, &
172 194242 : dummy_field, &
173 194242 : strtmp3, &
174 194242 : atom_info%atm_charge(index_now), &
175 388484 : atom_info%atm_mass(index_now)
176 194242 : atom_info%id_molname(index_now) = str2id(s2s(strtmp1))
177 194242 : atom_info%id_resname(index_now) = str2id(s2s(strtmp2))
178 194588 : atom_info%id_atmname(index_now) = str2id(s2s(strtmp3))
179 : END DO
180 : ELSE
181 185051 : DO iatom = 1, natom
182 177077 : index_now = iatom + natom_prev
183 177077 : CALL parser_get_next_line(parser, 1)
184 : READ (parser%input_line, FMT=psf_format) &
185 177077 : i, &
186 177077 : strtmp1, &
187 177077 : atom_info%resid(index_now), &
188 177077 : strtmp2, &
189 177077 : dummy_field, &
190 177077 : strtmp3, &
191 177077 : atom_info%atm_charge(index_now), &
192 177077 : atom_info%atm_mass(index_now), &
193 354154 : idum
194 177077 : atom_info%id_molname(index_now) = str2id(s2s(strtmp1))
195 177077 : atom_info%id_resname(index_now) = str2id(s2s(strtmp2))
196 185051 : atom_info%id_atmname(index_now) = str2id(s2s(ADJUSTL(strtmp3)))
197 : END DO
198 : END IF
199 : END IF
200 :
201 : !
202 : ! BOND section
203 : !
204 8320 : nbond_prev = 0
205 8320 : IF (ASSOCIATED(conn_info%bond_a)) nbond_prev = SIZE(conn_info%bond_a)
206 :
207 8320 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NBOND section'
208 8320 : IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated BOND: ', nbond_prev
209 8320 : label = '!NBOND'
210 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
211 8320 : IF (.NOT. found) THEN
212 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NBOND section '
213 0 : nbond = 0
214 : ELSE
215 8320 : CALL parser_get_object(parser, nbond)
216 8320 : IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NBOND = ', nbond
217 : !malloc the memory that we need
218 8320 : CALL reallocate(conn_info%bond_a, 1, nbond_prev + nbond)
219 8320 : CALL reallocate(conn_info%bond_b, 1, nbond_prev + nbond)
220 : !Read in the bond info
221 8320 : IF (psf_type == do_conn_psf_u) THEN
222 44250 : DO ibond = 1, nbond, 4
223 43904 : CALL parser_get_next_line(parser, 1)
224 43904 : index_now = nbond_prev + ibond - 1
225 175468 : READ (parser%input_line, FMT=*, ERR=9) (conn_info%bond_a(index_now + i), &
226 219372 : conn_info%bond_b(index_now + i), &
227 263622 : i=1, MIN(4, (nbond - ibond + 1)))
228 : END DO
229 : ELSE
230 7974 : DO ibond = 1, nbond, 4
231 40818 : CALL parser_get_next_line(parser, 1)
232 40818 : index_now = nbond_prev + ibond - 1
233 : READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') &
234 152917 : (conn_info%bond_a(index_now + i), &
235 193735 : conn_info%bond_b(index_now + i), &
236 237925 : i=1, MIN(4, (nbond - ibond + 1)))
237 : END DO
238 : END IF
239 : IF (ANY(conn_info%bond_a(nbond_prev + 1:) <= 0) .OR. &
240 : ANY(conn_info%bond_a(nbond_prev + 1:) > natom) .OR. &
241 1330180 : ANY(conn_info%bond_b(nbond_prev + 1:) <= 0) .OR. &
242 : ANY(conn_info%bond_b(nbond_prev + 1:) > natom)) THEN
243 0 : CPABORT("topology_read, invalid bond")
244 : END IF
245 336705 : conn_info%bond_a(nbond_prev + 1:) = conn_info%bond_a(nbond_prev + 1:) + natom_prev
246 336705 : conn_info%bond_b(nbond_prev + 1:) = conn_info%bond_b(nbond_prev + 1:) + natom_prev
247 : END IF
248 : !
249 : ! THETA section
250 : !
251 8320 : ntheta_prev = 0
252 8320 : IF (ASSOCIATED(conn_info%theta_a)) ntheta_prev = SIZE(conn_info%theta_a)
253 :
254 8320 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NTHETA section'
255 8320 : IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated THETA: ', ntheta_prev
256 8320 : label = '!NTHETA'
257 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
258 8320 : IF (.NOT. found) THEN
259 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NTHETA section '
260 0 : ntheta = 0
261 : ELSE
262 8320 : CALL parser_get_object(parser, ntheta)
263 8320 : IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NTHETA = ', ntheta
264 : !malloc the memory that we need
265 8320 : CALL reallocate(conn_info%theta_a, 1, ntheta_prev + ntheta)
266 8320 : CALL reallocate(conn_info%theta_b, 1, ntheta_prev + ntheta)
267 8320 : CALL reallocate(conn_info%theta_c, 1, ntheta_prev + ntheta)
268 : !Read in the bend info
269 8320 : IF (psf_type == do_conn_psf_u) THEN
270 28292 : DO itheta = 1, ntheta, 3
271 27946 : CALL parser_get_next_line(parser, 1)
272 27946 : index_now = ntheta_prev + itheta - 1
273 83778 : READ (parser%input_line, FMT=*, ERR=9) (conn_info%theta_a(index_now + i), &
274 83778 : conn_info%theta_b(index_now + i), &
275 111724 : conn_info%theta_c(index_now + i), &
276 140016 : i=1, MIN(3, (ntheta - itheta + 1)))
277 : END DO
278 : ELSE
279 7974 : DO itheta = 1, ntheta, 3
280 17045 : CALL parser_get_next_line(parser, 1)
281 17045 : index_now = ntheta_prev + itheta - 1
282 : READ (parser%input_line, FMT='(9'//TRIM(c_int)//')') &
283 42711 : (conn_info%theta_a(index_now + i), &
284 42711 : conn_info%theta_b(index_now + i), &
285 59756 : conn_info%theta_c(index_now + i), &
286 80195 : i=1, MIN(3, (ntheta - itheta + 1)))
287 : END DO
288 : END IF
289 134809 : conn_info%theta_a(ntheta_prev + 1:) = conn_info%theta_a(ntheta_prev + 1:) + natom_prev
290 134809 : conn_info%theta_b(ntheta_prev + 1:) = conn_info%theta_b(ntheta_prev + 1:) + natom_prev
291 134809 : conn_info%theta_c(ntheta_prev + 1:) = conn_info%theta_c(ntheta_prev + 1:) + natom_prev
292 : END IF
293 : !
294 : ! PHI section
295 : !
296 8320 : nphi_prev = 0
297 8320 : IF (ASSOCIATED(conn_info%phi_a)) nphi_prev = SIZE(conn_info%phi_a)
298 :
299 8320 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NPHI section'
300 8320 : IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated PHI: ', nphi_prev
301 8320 : label = '!NPHI'
302 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
303 8320 : IF (.NOT. found) THEN
304 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NPHI section '
305 : nphi = 0
306 : ELSE
307 8320 : CALL parser_get_object(parser, nphi)
308 8320 : IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NPHI = ', nphi
309 : !malloc the memory that we need
310 8320 : CALL reallocate(conn_info%phi_a, 1, nphi_prev + nphi)
311 8320 : CALL reallocate(conn_info%phi_b, 1, nphi_prev + nphi)
312 8320 : CALL reallocate(conn_info%phi_c, 1, nphi_prev + nphi)
313 8320 : CALL reallocate(conn_info%phi_d, 1, nphi_prev + nphi)
314 : !Read in the torsion info
315 8320 : IF (psf_type == do_conn_psf_u) THEN
316 23270 : DO iphi = 1, nphi, 2
317 22924 : CALL parser_get_next_line(parser, 1)
318 22924 : index_now = nphi_prev + iphi - 1
319 45832 : READ (parser%input_line, FMT=*, ERR=9) (conn_info%phi_a(index_now + i), &
320 45832 : conn_info%phi_b(index_now + i), &
321 45832 : conn_info%phi_c(index_now + i), &
322 68756 : conn_info%phi_d(index_now + i), &
323 92026 : i=1, MIN(2, (nphi - iphi + 1)))
324 : END DO
325 : ELSE
326 7974 : DO iphi = 1, nphi, 2
327 11187 : CALL parser_get_next_line(parser, 1)
328 11187 : index_now = nphi_prev + iphi - 1
329 : READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') &
330 21209 : (conn_info%phi_a(index_now + i), &
331 21209 : conn_info%phi_b(index_now + i), &
332 21209 : conn_info%phi_c(index_now + i), &
333 32396 : conn_info%phi_d(index_now + i), &
334 50059 : i=1, MIN(2, (nphi - iphi + 1)))
335 : END DO
336 : END IF
337 75361 : conn_info%phi_a(nphi_prev + 1:) = conn_info%phi_a(nphi_prev + 1:) + natom_prev
338 75361 : conn_info%phi_b(nphi_prev + 1:) = conn_info%phi_b(nphi_prev + 1:) + natom_prev
339 75361 : conn_info%phi_c(nphi_prev + 1:) = conn_info%phi_c(nphi_prev + 1:) + natom_prev
340 75361 : conn_info%phi_d(nphi_prev + 1:) = conn_info%phi_d(nphi_prev + 1:) + natom_prev
341 : END IF
342 : !
343 : ! IMPHI section
344 : !
345 8320 : nphi_prev = 0
346 8320 : IF (ASSOCIATED(conn_info%impr_a)) nphi_prev = SIZE(conn_info%impr_a)
347 :
348 8320 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NIMPHI section'
349 8320 : IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated IMPHI: ', nphi_prev
350 8320 : label = '!NIMPHI'
351 8320 : CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.)
352 8320 : IF (.NOT. found) THEN
353 0 : IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NIMPHI section '
354 : nphi = 0
355 : ELSE
356 8320 : CALL parser_get_object(parser, nphi)
357 8320 : IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NIMPR = ', nphi
358 : !malloc the memory that we need
359 8320 : CALL reallocate(conn_info%impr_a, 1, nphi_prev + nphi)
360 8320 : CALL reallocate(conn_info%impr_b, 1, nphi_prev + nphi)
361 8320 : CALL reallocate(conn_info%impr_c, 1, nphi_prev + nphi)
362 8320 : CALL reallocate(conn_info%impr_d, 1, nphi_prev + nphi)
363 : !Read in the improper torsion info
364 8320 : IF (psf_type == do_conn_psf_u) THEN
365 1674 : DO iphi = 1, nphi, 2
366 1328 : CALL parser_get_next_line(parser, 1)
367 1328 : index_now = nphi_prev + iphi - 1
368 2654 : READ (parser%input_line, FMT=*, ERR=9) (conn_info%impr_a(index_now + i), &
369 2654 : conn_info%impr_b(index_now + i), &
370 2654 : conn_info%impr_c(index_now + i), &
371 3982 : conn_info%impr_d(index_now + i), &
372 5656 : i=1, MIN(2, (nphi - iphi + 1)))
373 : END DO
374 : ELSE
375 7974 : DO iphi = 1, nphi, 2
376 430 : CALL parser_get_next_line(parser, 1)
377 430 : index_now = nphi_prev + iphi - 1
378 : READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') &
379 850 : (conn_info%impr_a(index_now + i), &
380 850 : conn_info%impr_b(index_now + i), &
381 850 : conn_info%impr_c(index_now + i), &
382 1280 : conn_info%impr_d(index_now + i), &
383 9644 : i=1, MIN(2, (nphi - iphi + 1)))
384 : END DO
385 : END IF
386 11824 : conn_info%impr_a(nphi_prev + 1:) = conn_info%impr_a(nphi_prev + 1:) + natom_prev
387 11824 : conn_info%impr_b(nphi_prev + 1:) = conn_info%impr_b(nphi_prev + 1:) + natom_prev
388 11824 : conn_info%impr_c(nphi_prev + 1:) = conn_info%impr_c(nphi_prev + 1:) + natom_prev
389 11824 : conn_info%impr_d(nphi_prev + 1:) = conn_info%impr_d(nphi_prev + 1:) + natom_prev
390 : END IF
391 :
392 8320 : CALL parser_release(parser)
393 8320 : CALL timestop(handle)
394 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
395 8320 : "PRINT%TOPOLOGY_INFO/PSF_INFO")
396 8320 : RETURN
397 : 9 CONTINUE
398 : ! Print error and exit
399 0 : IF (output_unit > 0) THEN
400 : WRITE (output_unit, '(T2,A)') &
401 0 : "PSF_INFO| Error while reading PSF using the unformatted PSF reading option!", &
402 0 : "PSF_INFO| Try using PSF instead of UPSF."
403 : END IF
404 :
405 0 : CPABORT("Error while reading PSF data!")
406 :
407 24960 : END SUBROUTINE read_topology_psf
408 :
409 : ! **************************************************************************************************
410 : !> \brief Post processing of PSF informations
411 : !> \param topology ...
412 : !> \param subsys_section ...
413 : ! **************************************************************************************************
414 783 : SUBROUTINE psf_post_process(topology, subsys_section)
415 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
416 : TYPE(section_vals_type), POINTER :: subsys_section
417 :
418 : CHARACTER(len=*), PARAMETER :: routineN = 'psf_post_process'
419 :
420 : INTEGER :: handle, i, iatom, ibond, ionfo, iw, &
421 : jatom, N, natom, nbond, nonfo, nphi, &
422 : ntheta
423 783 : TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bend_list, ex_bond_list
424 : TYPE(atom_info_type), POINTER :: atom_info
425 : TYPE(connectivity_info_type), POINTER :: conn_info
426 : TYPE(cp_logger_type), POINTER :: logger
427 :
428 783 : NULLIFY (logger)
429 1566 : logger => cp_get_default_logger()
430 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", &
431 783 : extension=".subsysLog")
432 783 : CALL timeset(routineN, handle)
433 783 : atom_info => topology%atom_info
434 783 : conn_info => topology%conn_info
435 : !
436 : ! PARA_RES structure
437 : !
438 783 : natom = 0
439 783 : nbond = 0
440 783 : i = 0
441 783 : IF (ASSOCIATED(atom_info%id_molname)) natom = SIZE(atom_info%id_molname)
442 783 : IF (ASSOCIATED(conn_info%bond_a)) nbond = SIZE(conn_info%bond_a)
443 783 : IF (ASSOCIATED(conn_info%c_bond_a)) i = SIZE(conn_info%c_bond_a)
444 386352 : DO ibond = 1, nbond
445 385569 : iatom = conn_info%bond_a(ibond)
446 385569 : jatom = conn_info%bond_b(ibond)
447 386352 : IF (topology%para_res) THEN
448 : IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
449 385441 : (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
450 : (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN
451 2397 : IF (iw > 0) WRITE (iw, '(T2,A,2I6)') "PSF_INFO| PARA_RES, bond between molecules atom ", &
452 306 : iatom, jatom
453 2244 : i = i + 1
454 2244 : CALL reallocate(conn_info%c_bond_a, 1, i)
455 2244 : CALL reallocate(conn_info%c_bond_b, 1, i)
456 2244 : conn_info%c_bond_a(i) = iatom
457 2244 : conn_info%c_bond_b(i) = jatom
458 : END IF
459 : ELSE
460 128 : CPASSERT(atom_info%id_molname(iatom) == atom_info%id_molname(jatom))
461 : END IF
462 : END DO
463 : !
464 : ! UB structure
465 : !
466 783 : ntheta = 0
467 783 : IF (ASSOCIATED(conn_info%theta_a)) ntheta = SIZE(conn_info%theta_a)
468 783 : CALL reallocate(conn_info%ub_a, 1, ntheta)
469 783 : CALL reallocate(conn_info%ub_b, 1, ntheta)
470 783 : CALL reallocate(conn_info%ub_c, 1, ntheta)
471 173292 : conn_info%ub_a(:) = conn_info%theta_a(:)
472 173292 : conn_info%ub_b(:) = conn_info%theta_b(:)
473 173292 : conn_info%ub_c(:) = conn_info%theta_c(:)
474 : !
475 : ! ONFO structure
476 : !
477 783 : nphi = 0
478 783 : nonfo = 0
479 783 : IF (ASSOCIATED(conn_info%phi_a)) nphi = SIZE(conn_info%phi_a)
480 783 : CALL reallocate(conn_info%onfo_a, 1, nphi)
481 783 : CALL reallocate(conn_info%onfo_b, 1, nphi)
482 105502 : conn_info%onfo_a(1:) = conn_info%phi_a(1:)
483 105502 : conn_info%onfo_b(1:) = conn_info%phi_d(1:)
484 : ! Reorder bonds
485 452506 : ALLOCATE (ex_bond_list(natom))
486 450940 : DO I = 1, natom
487 450940 : ALLOCATE (ex_bond_list(I)%array1(0))
488 : END DO
489 783 : N = 0
490 783 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
491 783 : CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, N)
492 : ! Reorder bends
493 451723 : ALLOCATE (ex_bend_list(natom))
494 450940 : DO I = 1, natom
495 450940 : ALLOCATE (ex_bend_list(I)%array1(0))
496 : END DO
497 783 : N = 0
498 783 : IF (ASSOCIATED(conn_info%theta_a)) N = SIZE(conn_info%theta_a)
499 783 : CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, N)
500 105502 : DO ionfo = 1, nphi
501 : ! Check if the torsion is not shared between angles or bonds
502 742483 : IF (ANY(ex_bond_list(conn_info%onfo_a(ionfo))%array1 == conn_info%onfo_b(ionfo)) .OR. &
503 : ANY(ex_bend_list(conn_info%onfo_a(ionfo))%array1 == conn_info%onfo_b(ionfo))) CYCLE
504 99201 : nonfo = nonfo + 1
505 99201 : conn_info%onfo_a(nonfo) = conn_info%onfo_a(ionfo)
506 105502 : conn_info%onfo_b(nonfo) = conn_info%onfo_b(ionfo)
507 : END DO
508 : ! deallocate bends
509 450940 : DO I = 1, natom
510 450940 : DEALLOCATE (ex_bend_list(I)%array1)
511 : END DO
512 783 : DEALLOCATE (ex_bend_list)
513 : ! deallocate bonds
514 450940 : DO I = 1, natom
515 450940 : DEALLOCATE (ex_bond_list(I)%array1)
516 : END DO
517 783 : DEALLOCATE (ex_bond_list)
518 : ! Get unique onfo
519 451723 : ALLOCATE (ex_bond_list(natom))
520 450940 : DO I = 1, natom
521 450940 : ALLOCATE (ex_bond_list(I)%array1(0))
522 : END DO
523 783 : N = 0
524 783 : IF (ASSOCIATED(conn_info%onfo_a)) N = nonfo
525 783 : CALL reorder_structure(ex_bond_list, conn_info%onfo_a, conn_info%onfo_b, N)
526 783 : nonfo = 0
527 450940 : DO I = 1, natom
528 649342 : DO ionfo = 1, SIZE(ex_bond_list(I)%array1)
529 1759715 : IF (COUNT(ex_bond_list(I)%array1 == ex_bond_list(I)%array1(ionfo)) /= 1) THEN
530 3240 : ex_bond_list(I)%array1(ionfo) = 0
531 : ELSE
532 195162 : IF (ex_bond_list(I)%array1(ionfo) <= I) CYCLE
533 97581 : nonfo = nonfo + 1
534 97581 : conn_info%onfo_a(nonfo) = I
535 97581 : conn_info%onfo_b(nonfo) = ex_bond_list(I)%array1(ionfo)
536 : END IF
537 : END DO
538 : END DO
539 450940 : DO I = 1, natom
540 450940 : DEALLOCATE (ex_bond_list(I)%array1)
541 : END DO
542 783 : DEALLOCATE (ex_bond_list)
543 783 : CALL reallocate(conn_info%onfo_a, 1, nonfo)
544 783 : CALL reallocate(conn_info%onfo_b, 1, nonfo)
545 :
546 783 : CALL timestop(handle)
547 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
548 783 : "PRINT%TOPOLOGY_INFO/PSF_INFO")
549 1566 : END SUBROUTINE psf_post_process
550 :
551 : ! **************************************************************************************************
552 : !> \brief Input driven modification (IDM) of PSF defined structures
553 : !> \param topology ...
554 : !> \param section ...
555 : !> \param subsys_section ...
556 : !> \author Teodoro Laino - Zurich University 04.2007
557 : ! **************************************************************************************************
558 819 : SUBROUTINE idm_psf(topology, section, subsys_section)
559 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
560 : TYPE(section_vals_type), POINTER :: section, subsys_section
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'idm_psf'
563 :
564 : INTEGER :: handle, i, iend, iend1, istart, istart1, &
565 : item, iw, j, mol_id, n_rep, natom, &
566 : nbond, nimpr, noe, nphi, ntheta
567 273 : INTEGER, DIMENSION(:), POINTER :: tag_mols, tmp, wrk
568 : LOGICAL :: explicit
569 273 : TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bond_list
570 : TYPE(atom_info_type), POINTER :: atom_info
571 : TYPE(connectivity_info_type), POINTER :: conn_info
572 : TYPE(cp_logger_type), POINTER :: logger
573 : TYPE(section_vals_type), POINTER :: subsection
574 :
575 273 : NULLIFY (logger)
576 546 : logger => cp_get_default_logger()
577 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", &
578 273 : extension=".subsysLog")
579 273 : CALL timeset(routineN, handle)
580 273 : CALL section_vals_get(section, explicit=explicit)
581 273 : IF (explicit) THEN
582 2 : atom_info => topology%atom_info
583 2 : conn_info => topology%conn_info
584 2 : natom = 0
585 2 : IF (ASSOCIATED(atom_info%id_molname)) natom = SIZE(atom_info%id_molname)
586 2 : nbond = 0
587 2 : IF (ASSOCIATED(conn_info%bond_a)) nbond = SIZE(conn_info%bond_a)
588 2 : ntheta = 0
589 2 : IF (ASSOCIATED(conn_info%theta_a)) ntheta = SIZE(conn_info%theta_a)
590 2 : nphi = 0
591 2 : IF (ASSOCIATED(conn_info%phi_a)) nphi = SIZE(conn_info%phi_a)
592 2 : nimpr = 0
593 2 : IF (ASSOCIATED(conn_info%impr_a)) nimpr = SIZE(conn_info%impr_a)
594 : ! Any new defined bond
595 2 : subsection => section_vals_get_subs_vals(section, "BONDS")
596 2 : CALL section_vals_get(subsection, explicit=explicit)
597 2 : IF (explicit) THEN
598 2 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
599 2 : CALL reallocate(conn_info%bond_a, 1, n_rep + nbond)
600 2 : CALL reallocate(conn_info%bond_b, 1, n_rep + nbond)
601 8 : DO i = 1, n_rep
602 6 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp)
603 6 : conn_info%bond_a(nbond + i) = tmp(1)
604 8 : conn_info%bond_b(nbond + i) = tmp(2)
605 : END DO
606 : ! And now modify the molecule name if two molecules have been bridged
607 6712 : ALLOCATE (ex_bond_list(natom))
608 6 : ALLOCATE (tag_mols(natom))
609 4 : ALLOCATE (wrk(natom))
610 6708 : DO j = 1, natom
611 6708 : ALLOCATE (ex_bond_list(j)%array1(0))
612 : END DO
613 2 : CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, nbond + n_rep)
614 : ! Loop over atoms to possiblyt change molecule name
615 6708 : tag_mols = -1
616 2 : mol_id = 1
617 6708 : DO i = 1, natom
618 6706 : IF (tag_mols(i) /= -1) CYCLE
619 1110 : CALL tag_molecule(tag_mols, ex_bond_list, i, mol_id)
620 6708 : mol_id = mol_id + 1
621 : END DO
622 2 : mol_id = mol_id - 1
623 2 : IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Number of molecules detected after merging: ', mol_id
624 : ! Now simply check about the contiguousness of molecule definition
625 2 : CALL sort(tag_mols, natom, wrk)
626 2 : item = tag_mols(1)
627 2 : istart = 1
628 6706 : DO i = 2, natom
629 6704 : IF (tag_mols(i) == item) CYCLE
630 1108 : iend = i - 1
631 1108 : noe = iend - istart + 1
632 7802 : istart1 = MINVAL(wrk(istart:iend))
633 7802 : iend1 = MAXVAL(wrk(istart:iend))
634 1108 : CPASSERT(iend1 - istart1 + 1 == noe)
635 7802 : atom_info%id_molname(istart1:iend1) = str2id(s2s("MOL"//cp_to_string(item)))
636 1108 : item = tag_mols(i)
637 6706 : istart = i
638 : END DO
639 2 : iend = i - 1
640 2 : noe = iend - istart + 1
641 14 : istart1 = MINVAL(wrk(istart:iend))
642 14 : iend1 = MAXVAL(wrk(istart:iend))
643 2 : CPASSERT(iend1 - istart1 + 1 == noe)
644 14 : atom_info%id_molname(istart1:iend1) = str2id(s2s("MOL"//cp_to_string(item)))
645 : ! Deallocate bonds
646 6708 : DO i = 1, natom
647 6708 : DEALLOCATE (ex_bond_list(i)%array1)
648 : END DO
649 2 : DEALLOCATE (ex_bond_list)
650 2 : DEALLOCATE (tag_mols)
651 4 : DEALLOCATE (wrk)
652 : END IF
653 : ! Any new defined angle
654 2 : subsection => section_vals_get_subs_vals(section, "ANGLES")
655 2 : CALL section_vals_get(subsection, explicit=explicit)
656 2 : IF (explicit) THEN
657 2 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
658 2 : CALL reallocate(conn_info%theta_a, 1, n_rep + ntheta)
659 2 : CALL reallocate(conn_info%theta_b, 1, n_rep + ntheta)
660 2 : CALL reallocate(conn_info%theta_c, 1, n_rep + ntheta)
661 26 : DO i = 1, n_rep
662 24 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp)
663 24 : conn_info%theta_a(ntheta + i) = tmp(1)
664 24 : conn_info%theta_b(ntheta + i) = tmp(2)
665 26 : conn_info%theta_c(ntheta + i) = tmp(3)
666 : END DO
667 : END IF
668 : ! Any new defined torsion
669 2 : subsection => section_vals_get_subs_vals(section, "TORSIONS")
670 2 : CALL section_vals_get(subsection, explicit=explicit)
671 2 : IF (explicit) THEN
672 2 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
673 2 : CALL reallocate(conn_info%phi_a, 1, n_rep + nphi)
674 2 : CALL reallocate(conn_info%phi_b, 1, n_rep + nphi)
675 2 : CALL reallocate(conn_info%phi_c, 1, n_rep + nphi)
676 2 : CALL reallocate(conn_info%phi_d, 1, n_rep + nphi)
677 74 : DO i = 1, n_rep
678 72 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp)
679 72 : conn_info%phi_a(nphi + i) = tmp(1)
680 72 : conn_info%phi_b(nphi + i) = tmp(2)
681 72 : conn_info%phi_c(nphi + i) = tmp(3)
682 74 : conn_info%phi_d(nphi + i) = tmp(4)
683 : END DO
684 : END IF
685 : ! Any new defined improper
686 2 : subsection => section_vals_get_subs_vals(section, "IMPROPERS")
687 2 : CALL section_vals_get(subsection, explicit=explicit)
688 2 : IF (explicit) THEN
689 0 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep)
690 0 : CALL reallocate(conn_info%impr_a, 1, n_rep + nimpr)
691 0 : CALL reallocate(conn_info%impr_b, 1, n_rep + nimpr)
692 0 : CALL reallocate(conn_info%impr_c, 1, n_rep + nimpr)
693 0 : CALL reallocate(conn_info%impr_d, 1, n_rep + nimpr)
694 0 : DO i = 1, n_rep
695 0 : CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp)
696 0 : conn_info%impr_a(nimpr + i) = tmp(1)
697 0 : conn_info%impr_b(nimpr + i) = tmp(2)
698 0 : conn_info%impr_c(nimpr + i) = tmp(3)
699 0 : conn_info%impr_d(nimpr + i) = tmp(4)
700 : END DO
701 : END IF
702 : END IF
703 273 : CALL timestop(handle)
704 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
705 273 : "PRINT%TOPOLOGY_INFO/PSF_INFO")
706 :
707 273 : END SUBROUTINE idm_psf
708 :
709 : ! **************************************************************************************************
710 : !> \brief Teodoro Laino - 01.2006
711 : !> Write PSF topology file in the CHARMM31 EXT standard format
712 : !> \param file_unit ...
713 : !> \param topology ...
714 : !> \param subsys_section ...
715 : !> \param force_env_section ...
716 : ! **************************************************************************************************
717 56 : SUBROUTINE write_topology_psf(file_unit, topology, subsys_section, force_env_section)
718 : INTEGER, INTENT(IN) :: file_unit
719 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
720 : TYPE(section_vals_type), POINTER :: subsys_section, force_env_section
721 :
722 : CHARACTER(len=*), PARAMETER :: routineN = 'write_topology_psf'
723 :
724 : CHARACTER(LEN=2*default_string_length) :: psf_format
725 : CHARACTER(LEN=default_path_length) :: record
726 : CHARACTER(LEN=default_string_length) :: c_int, my_tag1, my_tag2, my_tag3
727 : CHARACTER(LEN=default_string_length), &
728 56 : DIMENSION(:), POINTER :: charge_atm
729 : INTEGER :: handle, i, iw, j, my_index, nchg
730 : LOGICAL :: explicit, ldum
731 56 : REAL(KIND=dp), DIMENSION(:), POINTER :: charge_inp, charges
732 : TYPE(atom_info_type), POINTER :: atom_info
733 : TYPE(connectivity_info_type), POINTER :: conn_info
734 : TYPE(cp_logger_type), POINTER :: logger
735 : TYPE(section_vals_type), POINTER :: print_key, tmp_section
736 :
737 56 : NULLIFY (logger)
738 112 : logger => cp_get_default_logger()
739 56 : print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PSF")
740 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", &
741 56 : extension=".subsysLog")
742 56 : CALL timeset(routineN, handle)
743 :
744 56 : atom_info => topology%atom_info
745 56 : conn_info => topology%conn_info
746 :
747 : ! Check for charges.. (need to dump them in the PSF..)
748 168 : ALLOCATE (charges(topology%natoms))
749 57858 : charges = atom_info%atm_charge
750 : ! Collect charges from Input file..
751 56 : NULLIFY (tmp_section)
752 56 : tmp_section => section_vals_get_subs_vals(force_env_section, "MM%FORCEFIELD%CHARGE")
753 56 : CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg)
754 56 : IF (explicit) THEN
755 63 : ALLOCATE (charge_atm(nchg))
756 63 : ALLOCATE (charge_inp(nchg))
757 21 : CALL read_chrg_section(charge_atm, charge_inp, section=tmp_section, start=0)
758 3672 : DO j = 1, topology%natoms
759 3651 : record = id2str(atom_info%id_atmname(j))
760 3651 : ldum = qmmm_ff_precond_only_qm(record)
761 3651 : CALL uppercase(record)
762 8882 : DO i = 1, nchg
763 8861 : IF (record == charge_atm(i)) THEN
764 3651 : charges(j) = charge_inp(i)
765 3651 : EXIT
766 : END IF
767 : END DO
768 : END DO
769 21 : DEALLOCATE (charge_atm)
770 21 : DEALLOCATE (charge_inp)
771 : END IF
772 : ! fixup for topology output
773 28957 : DO j = 1, topology%natoms
774 28957 : IF (charges(j) == -HUGE(0.0_dp)) charges(j) = -99.0_dp
775 : END DO
776 : record = cp_print_key_generate_filename(logger, print_key, &
777 56 : extension=".psf", my_local=.FALSE.)
778 : ! build the EXT format
779 56 : c_int = "I10"
780 56 : psf_format = '(I10,T12,A,T21,I0,T30,A,T39,A,T47,A,T53,F10.6,T69,F8.3,T88,I1)'
781 56 : IF (iw > 0) WRITE (iw, '(T2,A)') &
782 0 : "PSF_WRITE| Writing out PSF file with CHARMM31 EXTErnal format: ", TRIM(record)
783 :
784 56 : WRITE (file_unit, FMT='(A)') "PSF EXT"
785 56 : WRITE (file_unit, FMT='(A)') ""
786 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 1, " !NTITLE"
787 56 : WRITE (file_unit, FMT='(A)') " CP2K generated DUMP of connectivity"
788 56 : WRITE (file_unit, FMT='(A)') ""
789 :
790 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') topology%natoms, " !NATOM"
791 56 : my_index = 1
792 56 : i = 1
793 56 : my_tag1 = id2str(atom_info%id_molname(i))
794 56 : my_tag2 = id2str(atom_info%id_resname(i))
795 56 : my_tag3 = id2str(atom_info%id_atmname(i))
796 56 : ldum = qmmm_ff_precond_only_qm(my_tag1)
797 56 : ldum = qmmm_ff_precond_only_qm(my_tag2)
798 56 : ldum = qmmm_ff_precond_only_qm(my_tag3)
799 : WRITE (file_unit, FMT=psf_format) &
800 56 : i, &
801 56 : TRIM(my_tag1), &
802 56 : my_index, &
803 56 : TRIM(my_tag2), &
804 56 : TRIM(my_tag3), &
805 56 : TRIM(my_tag3), &
806 56 : charges(i), &
807 56 : atom_info%atm_mass(i), &
808 112 : 0
809 28901 : DO i = 2, topology%natoms
810 28845 : IF ((atom_info%map_mol_num(i) /= atom_info%map_mol_num(i - 1)) .OR. &
811 8107 : (atom_info%map_mol_res(i) /= atom_info%map_mol_res(i - 1))) my_index = my_index + 1
812 28845 : my_tag1 = id2str(atom_info%id_molname(i))
813 28845 : my_tag2 = id2str(atom_info%id_resname(i))
814 28845 : my_tag3 = id2str(atom_info%id_atmname(i))
815 28845 : ldum = qmmm_ff_precond_only_qm(my_tag1)
816 28845 : ldum = qmmm_ff_precond_only_qm(my_tag2)
817 28845 : ldum = qmmm_ff_precond_only_qm(my_tag3)
818 : WRITE (file_unit, FMT=psf_format) &
819 28845 : i, &
820 28845 : TRIM(my_tag1), &
821 28845 : my_index, &
822 28845 : TRIM(my_tag2), &
823 28845 : TRIM(my_tag3), &
824 28845 : TRIM(my_tag3), &
825 28845 : charges(i), &
826 28845 : atom_info%atm_mass(i), &
827 57746 : 0
828 : END DO
829 56 : WRITE (file_unit, FMT='(/)')
830 56 : DEALLOCATE (charges)
831 :
832 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%bond_a), " !NBOND"
833 5251 : DO i = 1, SIZE(conn_info%bond_a), 4
834 5195 : j = 0
835 25927 : DO WHILE ((j < 4) .AND. ((i + j) <= SIZE(conn_info%bond_a)))
836 : WRITE (file_unit, FMT='(2('//TRIM(c_int)//'))', ADVANCE="NO") &
837 20732 : conn_info%bond_a(i + j), conn_info%bond_b(i + j)
838 20754 : j = j + 1
839 : END DO
840 5251 : WRITE (file_unit, FMT='(/)', ADVANCE="NO")
841 : END DO
842 56 : WRITE (file_unit, FMT='(/)')
843 :
844 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%theta_a), " !NTHETA"
845 5504 : DO i = 1, SIZE(conn_info%theta_a), 3
846 5448 : j = 0
847 21757 : DO WHILE ((j < 3) .AND. ((i + j) <= SIZE(conn_info%theta_a)))
848 : WRITE (file_unit, FMT='(3('//TRIM(c_int)//'))', ADVANCE="NO") &
849 16309 : conn_info%theta_a(i + j), conn_info%theta_b(i + j), &
850 32618 : conn_info%theta_c(i + j)
851 16333 : j = j + 1
852 : END DO
853 5504 : WRITE (file_unit, FMT='(/)', ADVANCE="NO")
854 : END DO
855 56 : WRITE (file_unit, FMT='(/)')
856 :
857 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%phi_a), " !NPHI"
858 4992 : DO i = 1, SIZE(conn_info%phi_a), 2
859 4936 : j = 0
860 14796 : DO WHILE ((j < 2) .AND. ((i + j) <= SIZE(conn_info%phi_a)))
861 : WRITE (file_unit, FMT='(4('//TRIM(c_int)//'))', ADVANCE="NO") &
862 9860 : conn_info%phi_a(i + j), conn_info%phi_b(i + j), &
863 19720 : conn_info%phi_c(i + j), conn_info%phi_d(i + j)
864 9872 : j = j + 1
865 : END DO
866 4992 : WRITE (file_unit, FMT='(/)', ADVANCE="NO")
867 : END DO
868 56 : WRITE (file_unit, FMT='(/)')
869 :
870 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%impr_a), " !NIMPHI"
871 364 : DO i = 1, SIZE(conn_info%impr_a), 2
872 308 : j = 0
873 920 : DO WHILE ((j < 2) .AND. ((i + j) <= SIZE(conn_info%impr_a)))
874 : WRITE (file_unit, FMT='(4('//TRIM(c_int)//'))', ADVANCE="NO") &
875 612 : conn_info%impr_a(i + j), conn_info%impr_b(i + j), &
876 1224 : conn_info%impr_c(i + j), conn_info%impr_d(i + j)
877 616 : j = j + 1
878 : END DO
879 364 : WRITE (file_unit, FMT='(/)', ADVANCE="NO")
880 : END DO
881 56 : WRITE (file_unit, FMT='(/)')
882 :
883 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NDON"
884 56 : WRITE (file_unit, FMT='(/)')
885 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NACC"
886 56 : WRITE (file_unit, FMT='(/)')
887 56 : WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NNB"
888 56 : WRITE (file_unit, FMT='(/)')
889 :
890 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
891 56 : "PRINT%TOPOLOGY_INFO/PSF_INFO")
892 56 : CALL timestop(handle)
893 :
894 168 : END SUBROUTINE write_topology_psf
895 :
896 : END MODULE topology_psf
897 :
|
