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