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