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 Handles PDB files
10 : !>
11 : !> PDB Format Description Version 2.2 from http://www.rcsb.org
12 : !> COLUMNS DATA TYPE FIELD DEFINITION
13 : !>
14 : !> 1 - 6 Record name "ATOM "
15 : !> 7 - 11 Integer serial Atom serial number.
16 : !> 13 - 16 Atom name Atom name.
17 : !> 17 Character altLoc Alternate location indicator.
18 : !> 18 - 20 Residue name resName Residue name.
19 : !> 22 Character chainID Chain identifier.
20 : !> 23 - 26 Integer resSeq Residue sequence number.
21 : !> 27 AChar iCode Code for insertion of residues.
22 : !> 31 - 38 Real(8.3) x Orthogonal coordinates for X in
23 : !> Angstroms.
24 : !> 39 - 46 Real(8.3) y Orthogonal coordinates for Y in
25 : !> Angstroms.
26 : !> 47 - 54 Real(8.3) z Orthogonal coordinates for Z in
27 : !> Angstroms.
28 : !> 55 - 60 Real(6.2) occupancy Occupancy.
29 : !> 61 - 66 Real(6.2) tempFactor Temperature factor.
30 : !> 73 - 76 LString(4) segID Segment identifier, left-justified.
31 : !> 77 - 78 LString(2) element Element symbol, right-justified.
32 : !> 79 - 80 LString(2) charge Charge on the atom.
33 : !>
34 : !> 81 - Real(*) Charge Ext. This last field is an extenstion to
35 : !> standard PDB to provide a full charge
36 : !> without limitation of digits.
37 : !>
38 : !> 1 - 6 Record name "CRYST1"
39 : !> 7 - 15 Real(9.3) a (Angstroms)
40 : !> 16 - 24 Real(9.3) b (Angstroms)
41 : !> 25 - 33 Real(9.3) c (Angstroms)
42 : !> 34 - 40 Real(7.2) alpha (degrees)
43 : !> 41 - 47 Real(7.2) beta (degrees)
44 : !> 48 - 54 Real(7.2) gamma (degrees)
45 : !> 56 - 66 LString Space group
46 : !> 67 - 70 Integer Z value
47 : ! **************************************************************************************************
48 : MODULE topology_pdb
49 : USE cell_types, ONLY: get_cell
50 : USE cp2k_info, ONLY: compile_revision,&
51 : cp2k_version,&
52 : r_host_name,&
53 : r_user_name
54 : USE cp_log_handling, ONLY: cp_get_default_logger,&
55 : cp_logger_type
56 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
57 : cp_print_key_generate_filename,&
58 : cp_print_key_unit_nr
59 : USE cp_parser_methods, ONLY: parser_get_next_line
60 : USE cp_parser_types, ONLY: cp_parser_type,&
61 : parser_create,&
62 : parser_release
63 : USE cp_units, ONLY: cp_unit_to_cp2k
64 : USE input_constants, ONLY: do_conn_user
65 : USE input_section_types, ONLY: section_get_rval,&
66 : section_vals_get_subs_vals,&
67 : section_vals_type,&
68 : section_vals_val_get
69 : USE kinds, ONLY: default_path_length,&
70 : default_string_length,&
71 : dp
72 : USE machine, ONLY: m_timestamp,&
73 : timestamp_length
74 : USE memory_utilities, ONLY: reallocate
75 : USE message_passing, ONLY: mp_para_env_type
76 : USE physcon, ONLY: angstrom
77 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
78 : USE string_table, ONLY: id2str,&
79 : s2s,&
80 : str2id
81 : USE topology_types, ONLY: atom_info_type,&
82 : topology_parameters_type
83 : #include "./base/base_uses.f90"
84 :
85 : IMPLICIT NONE
86 :
87 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_pdb'
88 :
89 : PRIVATE
90 : PUBLIC :: read_coordinate_pdb, write_coordinate_pdb
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief ...
96 : !> \param topology ...
97 : !> \param para_env ...
98 : !> \param subsys_section ...
99 : !> \par History
100 : !> TLAINO 05.2004 - Added the TER option to use different non-bonded molecules
101 : ! **************************************************************************************************
102 2322 : SUBROUTINE read_coordinate_pdb(topology, para_env, subsys_section)
103 : TYPE(topology_parameters_type) :: topology
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 : TYPE(section_vals_type), POINTER :: subsys_section
106 :
107 : CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_pdb'
108 : INTEGER, PARAMETER :: nblock = 1000
109 :
110 : CHARACTER(LEN=default_path_length) :: line
111 : CHARACTER(LEN=default_string_length) :: record, root_mol_name, strtmp
112 : INTEGER :: handle, id0, inum_mol, istat, iw, natom, &
113 : newsize
114 : LOGICAL :: my_end
115 : REAL(KIND=dp) :: pfactor
116 : TYPE(atom_info_type), POINTER :: atom_info
117 : TYPE(cp_logger_type), POINTER :: logger
118 : TYPE(cp_parser_type) :: parser
119 :
120 774 : NULLIFY (logger)
121 1548 : logger => cp_get_default_logger()
122 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
123 774 : extension=".subsysLog")
124 774 : CALL timeset(routineN, handle)
125 :
126 774 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
127 774 : atom_info => topology%atom_info
128 774 : CALL reallocate(atom_info%id_molname, 1, nblock)
129 774 : CALL reallocate(atom_info%id_resname, 1, nblock)
130 774 : CALL reallocate(atom_info%resid, 1, nblock)
131 774 : CALL reallocate(atom_info%id_atmname, 1, nblock)
132 774 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
133 774 : CALL reallocate(atom_info%atm_mass, 1, nblock)
134 774 : CALL reallocate(atom_info%atm_charge, 1, nblock)
135 774 : CALL reallocate(atom_info%occup, 1, nblock)
136 774 : CALL reallocate(atom_info%beta, 1, nblock)
137 774 : CALL reallocate(atom_info%id_element, 1, nblock)
138 :
139 774 : IF (iw > 0) THEN
140 : WRITE (UNIT=iw, FMT="(T2,A)") &
141 1 : "BEGIN of PDB data read from file "//TRIM(topology%coord_file_name)
142 : END IF
143 :
144 774 : id0 = str2id(s2s(""))
145 774 : topology%molname_generated = .FALSE.
146 :
147 774 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
148 :
149 774 : natom = 0
150 774 : inum_mol = 1
151 774 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
152 : DO
153 425289 : line = ""
154 425289 : CALL parser_get_next_line(parser, 1, at_end=my_end)
155 425289 : IF (my_end) EXIT
156 424779 : line = parser%input_line(1:default_path_length)
157 424779 : record = line(1:6)
158 : record = TRIM(record)
159 :
160 424779 : IF ((record == "ATOM") .OR. (record == "HETATM")) THEN
161 382781 : natom = natom + 1
162 382781 : topology%natoms = natom
163 382781 : IF (natom > SIZE(atom_info%id_atmname)) THEN
164 468 : newsize = INT(pfactor*natom)
165 468 : CALL reallocate(atom_info%id_molname, 1, newsize)
166 468 : CALL reallocate(atom_info%id_resname, 1, newsize)
167 468 : CALL reallocate(atom_info%resid, 1, newsize)
168 468 : CALL reallocate(atom_info%id_atmname, 1, newsize)
169 468 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
170 468 : CALL reallocate(atom_info%atm_mass, 1, newsize)
171 468 : CALL reallocate(atom_info%atm_charge, 1, newsize)
172 468 : CALL reallocate(atom_info%occup, 1, newsize)
173 468 : CALL reallocate(atom_info%beta, 1, newsize)
174 468 : CALL reallocate(atom_info%id_element, 1, newsize)
175 : END IF
176 : END IF
177 :
178 510 : SELECT CASE (record)
179 : CASE ("ATOM", "HETATM")
180 382781 : READ (UNIT=line(13:16), FMT=*) strtmp
181 382781 : atom_info%id_atmname(natom) = str2id(s2s(strtmp))
182 382781 : READ (UNIT=line(18:20), FMT=*, IOSTAT=istat) strtmp
183 382781 : IF (istat == 0) THEN
184 379615 : atom_info%id_resname(natom) = str2id(s2s(strtmp))
185 : ELSE
186 3166 : atom_info%id_resname(natom) = id0
187 : END IF
188 : ! Some information is not always given, so we mark it as used to prevent linters from crying
189 : ! regarding never a never used variable 'istat'
190 382781 : READ (UNIT=line(23:26), FMT=*, IOSTAT=istat) atom_info%resid(natom)
191 : MARK_USED(istat)
192 382781 : READ (UNIT=line(31:38), FMT=*, IOSTAT=istat) atom_info%r(1, natom)
193 : MARK_USED(istat)
194 382781 : READ (UNIT=line(39:46), FMT=*, IOSTAT=istat) atom_info%r(2, natom)
195 : MARK_USED(istat)
196 382781 : READ (UNIT=line(47:54), FMT=*, IOSTAT=istat) atom_info%r(3, natom)
197 : MARK_USED(istat)
198 382781 : READ (UNIT=line(55:60), FMT=*, IOSTAT=istat) atom_info%occup(natom)
199 : MARK_USED(istat)
200 382781 : READ (UNIT=line(61:66), FMT=*, IOSTAT=istat) atom_info%beta(natom)
201 : MARK_USED(istat)
202 382781 : READ (UNIT=line(73:76), FMT=*, IOSTAT=istat) strtmp
203 382781 : IF (istat == 0) THEN
204 224920 : atom_info%id_molname(natom) = str2id(s2s(strtmp))
205 : ELSE
206 157861 : atom_info%id_molname(natom) = str2id(s2s(root_mol_name))
207 157861 : topology%molname_generated = .TRUE.
208 : END IF
209 382781 : READ (UNIT=line(77:78), FMT=*, IOSTAT=istat) strtmp
210 382781 : IF (istat == 0) THEN
211 145181 : atom_info%id_element(natom) = str2id(s2s(strtmp))
212 : ELSE
213 237600 : atom_info%id_element(natom) = id0
214 : END IF
215 382781 : atom_info%atm_mass(natom) = 0.0_dp
216 382781 : atom_info%atm_charge(natom) = -HUGE(0.0_dp)
217 382781 : IF (topology%charge_occup) atom_info%atm_charge(natom) = atom_info%occup(natom)
218 382781 : IF (topology%charge_beta) atom_info%atm_charge(natom) = atom_info%beta(natom)
219 382781 : IF (topology%charge_extended) THEN
220 3188 : IF (LEN_TRIM(line) > 80) THEN
221 3188 : READ (UNIT=line(81:), FMT=*) atom_info%atm_charge(natom)
222 : END IF
223 : END IF
224 :
225 382781 : IF (atom_info%id_element(natom) == id0) THEN
226 : ! Element is assigned on the basis of the atm_name
227 237600 : topology%aa_element = .TRUE.
228 237600 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
229 : END IF
230 :
231 382781 : IF (iw > 0) THEN
232 : WRITE (UNIT=iw, FMT="(A6,I5,T13,A4,T18,A3,T23,I4,T31,3F8.3,T73,A4,T77,A2)") &
233 6 : record, natom, &
234 6 : TRIM(id2str(atom_info%id_atmname(natom))), &
235 6 : TRIM(id2str(atom_info%id_resname(natom))), &
236 6 : atom_info%resid(natom), &
237 6 : atom_info%r(1, natom), &
238 6 : atom_info%r(2, natom), &
239 6 : atom_info%r(3, natom), &
240 6 : ADJUSTL(TRIM(id2str(atom_info%id_molname(natom)))), &
241 12 : ADJUSTR(TRIM(id2str(atom_info%id_element(natom))))
242 : END IF
243 382781 : atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "angstrom")
244 382781 : atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "angstrom")
245 382781 : atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "angstrom")
246 : CASE ("TER")
247 41348 : inum_mol = inum_mol + 1
248 41348 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
249 : CASE ("REMARK")
250 280 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) TRIM(line)
251 : CASE ("END")
252 424779 : EXIT
253 : CASE DEFAULT
254 : END SELECT
255 : END DO
256 774 : CALL parser_release(parser)
257 :
258 774 : CALL reallocate(atom_info%id_molname, 1, natom)
259 774 : CALL reallocate(atom_info%id_resname, 1, natom)
260 774 : CALL reallocate(atom_info%resid, 1, natom)
261 774 : CALL reallocate(atom_info%id_atmname, 1, natom)
262 774 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
263 774 : CALL reallocate(atom_info%atm_mass, 1, natom)
264 774 : CALL reallocate(atom_info%atm_charge, 1, natom)
265 774 : CALL reallocate(atom_info%occup, 1, natom)
266 774 : CALL reallocate(atom_info%beta, 1, natom)
267 774 : CALL reallocate(atom_info%id_element, 1, natom)
268 :
269 774 : IF (topology%conn_type /= do_conn_user) THEN
270 910 : IF (.NOT. topology%para_res) atom_info%resid(:) = 1
271 : END IF
272 :
273 774 : IF (iw > 0) THEN
274 : WRITE (UNIT=iw, FMT="(T2,A)") &
275 1 : "END of PDB data read from file "//TRIM(topology%coord_file_name)
276 : END IF
277 :
278 774 : topology%natoms = natom
279 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
280 774 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
281 774 : CALL timestop(handle)
282 :
283 2322 : END SUBROUTINE read_coordinate_pdb
284 :
285 : ! **************************************************************************************************
286 : !> \brief ...
287 : !> \param file_unit ...
288 : !> \param topology ...
289 : !> \param subsys_section ...
290 : ! **************************************************************************************************
291 330 : SUBROUTINE write_coordinate_pdb(file_unit, topology, subsys_section)
292 :
293 : INTEGER, INTENT(IN) :: file_unit
294 : TYPE(topology_parameters_type) :: topology
295 : TYPE(section_vals_type), POINTER :: subsys_section
296 :
297 : CHARACTER(len=*), PARAMETER :: routineN = 'write_coordinate_pdb'
298 :
299 : CHARACTER(LEN=120) :: line
300 : CHARACTER(LEN=default_path_length) :: record
301 : CHARACTER(LEN=default_string_length) :: my_tag1, my_tag2, my_tag3, my_tag4
302 : CHARACTER(LEN=timestamp_length) :: timestamp
303 : INTEGER :: handle, i, id1, id2, idres, iw, natom
304 : LOGICAL :: charge_beta, charge_extended, &
305 : charge_occup, ldum
306 : REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma
307 : REAL(KIND=dp), DIMENSION(3) :: abc
308 : TYPE(atom_info_type), POINTER :: atom_info
309 : TYPE(cp_logger_type), POINTER :: logger
310 : TYPE(section_vals_type), POINTER :: print_key
311 :
312 55 : NULLIFY (logger)
313 55 : logger => cp_get_default_logger()
314 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
315 55 : extension=".subsysLog")
316 55 : print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PDB")
317 55 : CALL timeset(routineN, handle)
318 :
319 55 : CALL section_vals_val_get(print_key, "CHARGE_OCCUP", l_val=charge_occup)
320 55 : CALL section_vals_val_get(print_key, "CHARGE_BETA", l_val=charge_beta)
321 55 : CALL section_vals_val_get(print_key, "CHARGE_EXTENDED", l_val=charge_extended)
322 220 : i = COUNT([charge_occup, charge_beta, charge_extended])
323 55 : IF (i > 1) THEN
324 0 : CPABORT("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
325 : END IF
326 :
327 55 : atom_info => topology%atom_info
328 : record = cp_print_key_generate_filename(logger, print_key, &
329 : extension=".pdb", &
330 55 : my_local=.FALSE.)
331 :
332 55 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) " Writing out PDB file ", TRIM(record)
333 :
334 : ! Write file header
335 55 : CALL m_timestamp(timestamp)
336 : WRITE (UNIT=file_unit, FMT="(A6,T11,A)") &
337 55 : "TITLE ", "PDB file created by "//TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")", &
338 110 : "AUTHOR", TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
339 : ! Write cell information
340 55 : CALL get_cell(cell=topology%cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc)
341 : WRITE (UNIT=file_unit, FMT="(A6,3F9.3,3F7.2)") &
342 220 : "CRYST1", abc(1:3)*angstrom, angle_alpha, angle_beta, angle_gamma
343 :
344 55 : natom = topology%natoms
345 55 : idres = 0
346 55 : id1 = 0
347 55 : id2 = 0
348 :
349 29019 : DO i = 1, natom
350 :
351 28964 : IF (topology%para_res) THEN
352 28268 : idres = atom_info%resid(i)
353 : ELSE
354 696 : IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i))) THEN
355 52 : idres = idres + 1
356 52 : id1 = atom_info%map_mol_num(i)
357 52 : id2 = atom_info%map_mol_typ(i)
358 : END IF
359 : END IF
360 :
361 28964 : line = ""
362 28964 : my_tag1 = id2str(atom_info%id_atmname(i)); ldum = qmmm_ff_precond_only_qm(my_tag1)
363 28964 : my_tag2 = id2str(atom_info%id_resname(i)); ldum = qmmm_ff_precond_only_qm(my_tag2)
364 28964 : my_tag3 = id2str(atom_info%id_molname(i)); ldum = qmmm_ff_precond_only_qm(my_tag3)
365 28964 : my_tag4 = id2str(atom_info%id_element(i)); ldum = qmmm_ff_precond_only_qm(my_tag4)
366 :
367 28964 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
368 28964 : WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(i, 100000)
369 28964 : WRITE (UNIT=line(13:16), FMT="(A4)") ADJUSTL(my_tag1(1:4))
370 28964 : WRITE (UNIT=line(18:20), FMT="(A3)") TRIM(my_tag2)
371 28964 : WRITE (UNIT=line(23:26), FMT="(I4)") MODULO(idres, 10000)
372 115856 : WRITE (UNIT=line(31:54), FMT="(3F8.3)") atom_info%r(1:3, i)*angstrom
373 28964 : IF (ASSOCIATED(atom_info%occup)) THEN
374 28652 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%occup(i)
375 : ELSE
376 312 : WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp
377 : END IF
378 28964 : IF (ASSOCIATED(atom_info%beta)) THEN
379 28652 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%beta(i)
380 : ELSE
381 312 : WRITE (UNIT=line(61:66), FMT="(F6.2)") 0.0_dp
382 : END IF
383 28964 : IF (ASSOCIATED(atom_info%atm_charge)) THEN
384 28964 : IF (ANY([charge_occup, charge_beta, charge_extended]) .AND. &
385 : (atom_info%atm_charge(i) == -HUGE(0.0_dp))) THEN
386 0 : CPABORT("No atomic charges found yet (after the topology setup)")
387 : END IF
388 28964 : IF (charge_occup) THEN
389 0 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%atm_charge(i)
390 28964 : ELSE IF (charge_beta) THEN
391 0 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%atm_charge(i)
392 28964 : ELSE IF (charge_extended) THEN
393 0 : WRITE (UNIT=line(81:), FMT="(F20.16)") atom_info%atm_charge(i)
394 : ELSE
395 : ! Write no atomic charge
396 : END IF
397 : END IF
398 28964 : WRITE (UNIT=line(73:76), FMT="(A4)") ADJUSTL(my_tag3)
399 28964 : WRITE (UNIT=line(77:78), FMT="(A2)") TRIM(my_tag4)
400 29019 : WRITE (UNIT=file_unit, FMT="(A)") TRIM(line)
401 : END DO
402 55 : WRITE (UNIT=file_unit, FMT="(A3)") "END"
403 :
404 55 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) " Exiting "//routineN
405 :
406 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
407 55 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
408 :
409 55 : CALL timestop(handle)
410 :
411 55 : END SUBROUTINE write_coordinate_pdb
412 :
413 : END MODULE topology_pdb
|
