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 2304 : 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 768 : NULLIFY (logger)
121 1536 : logger => cp_get_default_logger()
122 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
123 768 : extension=".subsysLog")
124 768 : CALL timeset(routineN, handle)
125 :
126 768 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
127 768 : atom_info => topology%atom_info
128 768 : CALL reallocate(atom_info%id_molname, 1, nblock)
129 768 : CALL reallocate(atom_info%id_resname, 1, nblock)
130 768 : CALL reallocate(atom_info%resid, 1, nblock)
131 768 : CALL reallocate(atom_info%id_atmname, 1, nblock)
132 768 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
133 768 : CALL reallocate(atom_info%atm_mass, 1, nblock)
134 768 : CALL reallocate(atom_info%atm_charge, 1, nblock)
135 768 : CALL reallocate(atom_info%occup, 1, nblock)
136 768 : CALL reallocate(atom_info%beta, 1, nblock)
137 768 : CALL reallocate(atom_info%id_element, 1, nblock)
138 :
139 768 : 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 768 : id0 = str2id(s2s(""))
145 768 : topology%molname_generated = .FALSE.
146 :
147 768 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
148 :
149 768 : natom = 0
150 768 : inum_mol = 1
151 768 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
152 : DO
153 425259 : line = ""
154 425259 : CALL parser_get_next_line(parser, 1, at_end=my_end)
155 425259 : IF (my_end) EXIT
156 424755 : line = parser%input_line(1:default_path_length)
157 424755 : record = line(1:6)
158 : record = TRIM(record)
159 :
160 424755 : IF ((record == "ATOM") .OR. (record == "HETATM")) THEN
161 382757 : natom = natom + 1
162 382757 : topology%natoms = natom
163 382757 : 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 504 : SELECT CASE (record)
179 : CASE ("ATOM", "HETATM")
180 382757 : READ (UNIT=line(13:16), FMT=*) strtmp
181 382757 : atom_info%id_atmname(natom) = str2id(s2s(strtmp))
182 382757 : READ (UNIT=line(18:20), FMT=*, IOSTAT=istat) strtmp
183 382757 : IF (istat == 0) THEN
184 379591 : 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 382757 : READ (UNIT=line(23:26), FMT=*, IOSTAT=istat) atom_info%resid(natom)
191 : MARK_USED(istat)
192 382757 : READ (UNIT=line(31:38), FMT=*, IOSTAT=istat) atom_info%r(1, natom)
193 : MARK_USED(istat)
194 382757 : READ (UNIT=line(39:46), FMT=*, IOSTAT=istat) atom_info%r(2, natom)
195 : MARK_USED(istat)
196 382757 : READ (UNIT=line(47:54), FMT=*, IOSTAT=istat) atom_info%r(3, natom)
197 : MARK_USED(istat)
198 382757 : READ (UNIT=line(55:60), FMT=*, IOSTAT=istat) atom_info%occup(natom)
199 : MARK_USED(istat)
200 382757 : READ (UNIT=line(61:66), FMT=*, IOSTAT=istat) atom_info%beta(natom)
201 : MARK_USED(istat)
202 382757 : READ (UNIT=line(73:76), FMT=*, IOSTAT=istat) strtmp
203 382757 : IF (istat == 0) THEN
204 224896 : 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 382757 : READ (UNIT=line(77:78), FMT=*, IOSTAT=istat) strtmp
210 382757 : IF (istat == 0) THEN
211 145157 : atom_info%id_element(natom) = str2id(s2s(strtmp))
212 : ELSE
213 237600 : atom_info%id_element(natom) = id0
214 : END IF
215 382757 : atom_info%atm_mass(natom) = 0.0_dp
216 382757 : atom_info%atm_charge(natom) = -HUGE(0.0_dp)
217 382757 : IF (topology%charge_occup) atom_info%atm_charge(natom) = atom_info%occup(natom)
218 382757 : IF (topology%charge_beta) atom_info%atm_charge(natom) = atom_info%beta(natom)
219 382757 : IF (topology%charge_extended) THEN
220 3188 : READ (UNIT=line(81:), FMT=*) atom_info%atm_charge(natom)
221 : END IF
222 :
223 382757 : IF (atom_info%id_element(natom) == id0) THEN
224 : ! Element is assigned on the basis of the atm_name
225 237600 : topology%aa_element = .TRUE.
226 237600 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
227 : END IF
228 :
229 382757 : IF (iw > 0) THEN
230 : WRITE (UNIT=iw, FMT="(A6,I5,T13,A4,T18,A3,T23,I4,T31,3F8.3,T73,A4,T77,A2)") &
231 6 : record, natom, &
232 6 : TRIM(id2str(atom_info%id_atmname(natom))), &
233 6 : TRIM(id2str(atom_info%id_resname(natom))), &
234 6 : atom_info%resid(natom), &
235 6 : atom_info%r(1, natom), &
236 6 : atom_info%r(2, natom), &
237 6 : atom_info%r(3, natom), &
238 6 : ADJUSTL(TRIM(id2str(atom_info%id_molname(natom)))), &
239 12 : ADJUSTR(TRIM(id2str(atom_info%id_element(natom))))
240 : END IF
241 382757 : atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "angstrom")
242 382757 : atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "angstrom")
243 382757 : atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "angstrom")
244 : CASE ("TER")
245 41348 : inum_mol = inum_mol + 1
246 41348 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
247 : CASE ("REMARK")
248 280 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) TRIM(line)
249 : CASE ("END")
250 424755 : EXIT
251 : CASE DEFAULT
252 : END SELECT
253 : END DO
254 768 : CALL parser_release(parser)
255 :
256 768 : CALL reallocate(atom_info%id_molname, 1, natom)
257 768 : CALL reallocate(atom_info%id_resname, 1, natom)
258 768 : CALL reallocate(atom_info%resid, 1, natom)
259 768 : CALL reallocate(atom_info%id_atmname, 1, natom)
260 768 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
261 768 : CALL reallocate(atom_info%atm_mass, 1, natom)
262 768 : CALL reallocate(atom_info%atm_charge, 1, natom)
263 768 : CALL reallocate(atom_info%occup, 1, natom)
264 768 : CALL reallocate(atom_info%beta, 1, natom)
265 768 : CALL reallocate(atom_info%id_element, 1, natom)
266 :
267 768 : IF (topology%conn_type /= do_conn_user) THEN
268 904 : IF (.NOT. topology%para_res) atom_info%resid(:) = 1
269 : END IF
270 :
271 768 : IF (iw > 0) THEN
272 : WRITE (UNIT=iw, FMT="(T2,A)") &
273 1 : "END of PDB data read from file "//TRIM(topology%coord_file_name)
274 : END IF
275 :
276 768 : topology%natoms = natom
277 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
278 768 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
279 768 : CALL timestop(handle)
280 :
281 2304 : END SUBROUTINE read_coordinate_pdb
282 :
283 : ! **************************************************************************************************
284 : !> \brief ...
285 : !> \param file_unit ...
286 : !> \param topology ...
287 : !> \param subsys_section ...
288 : ! **************************************************************************************************
289 324 : SUBROUTINE write_coordinate_pdb(file_unit, topology, subsys_section)
290 :
291 : INTEGER, INTENT(IN) :: file_unit
292 : TYPE(topology_parameters_type) :: topology
293 : TYPE(section_vals_type), POINTER :: subsys_section
294 :
295 : CHARACTER(len=*), PARAMETER :: routineN = 'write_coordinate_pdb'
296 :
297 : CHARACTER(LEN=120) :: line
298 : CHARACTER(LEN=default_path_length) :: record
299 : CHARACTER(LEN=default_string_length) :: my_tag1, my_tag2, my_tag3, my_tag4
300 : CHARACTER(LEN=timestamp_length) :: timestamp
301 : INTEGER :: handle, i, id1, id2, idres, iw, natom
302 : LOGICAL :: charge_beta, charge_extended, &
303 : charge_occup, ldum
304 : REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma
305 : REAL(KIND=dp), DIMENSION(3) :: abc
306 : TYPE(atom_info_type), POINTER :: atom_info
307 : TYPE(cp_logger_type), POINTER :: logger
308 : TYPE(section_vals_type), POINTER :: print_key
309 :
310 54 : NULLIFY (logger)
311 54 : logger => cp_get_default_logger()
312 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
313 54 : extension=".subsysLog")
314 54 : print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PDB")
315 54 : CALL timeset(routineN, handle)
316 :
317 54 : CALL section_vals_val_get(print_key, "CHARGE_OCCUP", l_val=charge_occup)
318 54 : CALL section_vals_val_get(print_key, "CHARGE_BETA", l_val=charge_beta)
319 54 : CALL section_vals_val_get(print_key, "CHARGE_EXTENDED", l_val=charge_extended)
320 216 : i = COUNT([charge_occup, charge_beta, charge_extended])
321 54 : IF (i > 1) THEN
322 0 : CPABORT("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
323 : END IF
324 :
325 54 : atom_info => topology%atom_info
326 : record = cp_print_key_generate_filename(logger, print_key, &
327 : extension=".pdb", &
328 54 : my_local=.FALSE.)
329 :
330 54 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) " Writing out PDB file ", TRIM(record)
331 :
332 : ! Write file header
333 54 : CALL m_timestamp(timestamp)
334 : WRITE (UNIT=file_unit, FMT="(A6,T11,A)") &
335 54 : "TITLE ", "PDB file created by "//TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")", &
336 108 : "AUTHOR", TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//timestamp(:19)
337 : ! Write cell information
338 54 : CALL get_cell(cell=topology%cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc)
339 : WRITE (UNIT=file_unit, FMT="(A6,3F9.3,3F7.2)") &
340 216 : "CRYST1", abc(1:3)*angstrom, angle_alpha, angle_beta, angle_gamma
341 :
342 54 : natom = topology%natoms
343 54 : idres = 0
344 54 : id1 = 0
345 54 : id2 = 0
346 :
347 28988 : DO i = 1, natom
348 :
349 28934 : IF (topology%para_res) THEN
350 28238 : idres = atom_info%resid(i)
351 : ELSE
352 696 : IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i))) THEN
353 52 : idres = idres + 1
354 52 : id1 = atom_info%map_mol_num(i)
355 52 : id2 = atom_info%map_mol_typ(i)
356 : END IF
357 : END IF
358 :
359 28934 : line = ""
360 28934 : my_tag1 = id2str(atom_info%id_atmname(i)); ldum = qmmm_ff_precond_only_qm(my_tag1)
361 28934 : my_tag2 = id2str(atom_info%id_resname(i)); ldum = qmmm_ff_precond_only_qm(my_tag2)
362 28934 : my_tag3 = id2str(atom_info%id_molname(i)); ldum = qmmm_ff_precond_only_qm(my_tag3)
363 28934 : my_tag4 = id2str(atom_info%id_element(i)); ldum = qmmm_ff_precond_only_qm(my_tag4)
364 :
365 28934 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
366 28934 : WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(i, 100000)
367 28934 : WRITE (UNIT=line(13:16), FMT="(A4)") ADJUSTL(my_tag1(1:4))
368 28934 : WRITE (UNIT=line(18:20), FMT="(A3)") TRIM(my_tag2)
369 28934 : WRITE (UNIT=line(23:26), FMT="(I4)") MODULO(idres, 10000)
370 115736 : WRITE (UNIT=line(31:54), FMT="(3F8.3)") atom_info%r(1:3, i)*angstrom
371 28934 : IF (ASSOCIATED(atom_info%occup)) THEN
372 28652 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%occup(i)
373 : ELSE
374 282 : WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp
375 : END IF
376 28934 : IF (ASSOCIATED(atom_info%beta)) THEN
377 28652 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%beta(i)
378 : ELSE
379 282 : WRITE (UNIT=line(61:66), FMT="(F6.2)") 0.0_dp
380 : END IF
381 28934 : IF (ASSOCIATED(atom_info%atm_charge)) THEN
382 28934 : IF (ANY([charge_occup, charge_beta, charge_extended]) .AND. &
383 : (atom_info%atm_charge(i) == -HUGE(0.0_dp))) THEN
384 0 : CPABORT("No atomic charges found yet (after the topology setup)")
385 : END IF
386 28934 : IF (charge_occup) THEN
387 0 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%atm_charge(i)
388 28934 : ELSE IF (charge_beta) THEN
389 0 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%atm_charge(i)
390 28934 : ELSE IF (charge_extended) THEN
391 0 : WRITE (UNIT=line(81:), FMT="(F20.16)") atom_info%atm_charge(i)
392 : ELSE
393 : ! Write no atomic charge
394 : END IF
395 : END IF
396 28934 : WRITE (UNIT=line(73:76), FMT="(A4)") ADJUSTL(my_tag3)
397 28934 : WRITE (UNIT=line(77:78), FMT="(A2)") TRIM(my_tag4)
398 28988 : WRITE (UNIT=file_unit, FMT="(A)") TRIM(line)
399 : END DO
400 54 : WRITE (UNIT=file_unit, FMT="(A3)") "END"
401 :
402 54 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) " Exiting "//routineN
403 :
404 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
405 54 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
406 :
407 54 : CALL timestop(handle)
408 :
409 54 : END SUBROUTINE write_coordinate_pdb
410 :
411 : END MODULE topology_pdb
|
