(git:07a6c39)
Loading...
Searching...
No Matches
topology_pdb.F
Go to the documentation of this file.
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! **************************************************************************************************
49 USE cell_types, ONLY: get_cell
50 USE cp2k_info, ONLY: compile_revision,&
63 USE cp_units, ONLY: cp_unit_to_cp2k
69 USE kinds, ONLY: default_path_length,&
71 dp
72 USE machine, ONLY: m_timestamp,&
76 USE physcon, ONLY: angstrom
78 USE string_table, ONLY: id2str,&
79 s2s,&
80 str2id
83#include "./base/base_uses.f90"
84
85 IMPLICIT NONE
86
87 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_pdb'
88
89 PRIVATE
91
92CONTAINS
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 SUBROUTINE read_coordinate_pdb(topology, para_env, subsys_section)
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 NULLIFY (logger)
121 logger => cp_get_default_logger()
122 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
123 extension=".subsysLog")
124 CALL timeset(routinen, handle)
125
126 pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
127 atom_info => topology%atom_info
128 CALL reallocate(atom_info%id_molname, 1, nblock)
129 CALL reallocate(atom_info%id_resname, 1, nblock)
130 CALL reallocate(atom_info%resid, 1, nblock)
131 CALL reallocate(atom_info%id_atmname, 1, nblock)
132 CALL reallocate(atom_info%r, 1, 3, 1, nblock)
133 CALL reallocate(atom_info%atm_mass, 1, nblock)
134 CALL reallocate(atom_info%atm_charge, 1, nblock)
135 CALL reallocate(atom_info%occup, 1, nblock)
136 CALL reallocate(atom_info%beta, 1, nblock)
137 CALL reallocate(atom_info%id_element, 1, nblock)
138
139 IF (iw > 0) THEN
140 WRITE (unit=iw, fmt="(T2,A)") &
141 "BEGIN of PDB data read from file "//trim(topology%coord_file_name)
142 END IF
143
144 id0 = str2id(s2s(""))
145 topology%molname_generated = .false.
146
147 CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
148
149 natom = 0
150 inum_mol = 1
151 WRITE (unit=root_mol_name, fmt='(A3,I0)') "MOL", inum_mol
152 DO
153 line = ""
154 CALL parser_get_next_line(parser, 1, at_end=my_end)
155 IF (my_end) EXIT
156 line = parser%input_line(1:default_path_length)
157 record = line(1:6)
158 record = trim(record)
159
160 IF ((record == "ATOM") .OR. (record == "HETATM")) THEN
161 natom = natom + 1
162 topology%natoms = natom
163 IF (natom > SIZE(atom_info%id_atmname)) THEN
164 newsize = int(pfactor*natom)
165 CALL reallocate(atom_info%id_molname, 1, newsize)
166 CALL reallocate(atom_info%id_resname, 1, newsize)
167 CALL reallocate(atom_info%resid, 1, newsize)
168 CALL reallocate(atom_info%id_atmname, 1, newsize)
169 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
170 CALL reallocate(atom_info%atm_mass, 1, newsize)
171 CALL reallocate(atom_info%atm_charge, 1, newsize)
172 CALL reallocate(atom_info%occup, 1, newsize)
173 CALL reallocate(atom_info%beta, 1, newsize)
174 CALL reallocate(atom_info%id_element, 1, newsize)
175 END IF
176 END IF
177
178 SELECT CASE (record)
179 CASE ("ATOM", "HETATM")
180 READ (unit=line(13:16), fmt=*) strtmp
181 atom_info%id_atmname(natom) = str2id(s2s(strtmp))
182 READ (unit=line(18:20), fmt=*, iostat=istat) strtmp
183 IF (istat == 0) THEN
184 atom_info%id_resname(natom) = str2id(s2s(strtmp))
185 ELSE
186 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 READ (unit=line(23:26), fmt=*, iostat=istat) atom_info%resid(natom)
191 mark_used(istat)
192 READ (unit=line(31:38), fmt=*, iostat=istat) atom_info%r(1, natom)
193 mark_used(istat)
194 READ (unit=line(39:46), fmt=*, iostat=istat) atom_info%r(2, natom)
195 mark_used(istat)
196 READ (unit=line(47:54), fmt=*, iostat=istat) atom_info%r(3, natom)
197 mark_used(istat)
198 READ (unit=line(55:60), fmt=*, iostat=istat) atom_info%occup(natom)
199 mark_used(istat)
200 READ (unit=line(61:66), fmt=*, iostat=istat) atom_info%beta(natom)
201 mark_used(istat)
202 READ (unit=line(73:76), fmt=*, iostat=istat) strtmp
203 IF (istat == 0) THEN
204 atom_info%id_molname(natom) = str2id(s2s(strtmp))
205 ELSE
206 atom_info%id_molname(natom) = str2id(s2s(root_mol_name))
207 topology%molname_generated = .true.
208 END IF
209 READ (unit=line(77:78), fmt=*, iostat=istat) strtmp
210 IF (istat == 0) THEN
211 atom_info%id_element(natom) = str2id(s2s(strtmp))
212 ELSE
213 atom_info%id_element(natom) = id0
214 END IF
215 atom_info%atm_mass(natom) = 0.0_dp
216 atom_info%atm_charge(natom) = -huge(0.0_dp)
217 IF (topology%charge_occup) atom_info%atm_charge(natom) = atom_info%occup(natom)
218 IF (topology%charge_beta) atom_info%atm_charge(natom) = atom_info%beta(natom)
219 IF (topology%charge_extended) THEN
220 IF (len_trim(line) > 80) THEN
221 READ (unit=line(81:), fmt=*) atom_info%atm_charge(natom)
222 END IF
223 END IF
224
225 IF (atom_info%id_element(natom) == id0) THEN
226 ! Element is assigned on the basis of the atm_name
227 topology%aa_element = .true.
228 atom_info%id_element(natom) = atom_info%id_atmname(natom)
229 END IF
230
231 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 record, natom, &
234 trim(id2str(atom_info%id_atmname(natom))), &
235 trim(id2str(atom_info%id_resname(natom))), &
236 atom_info%resid(natom), &
237 atom_info%r(1, natom), &
238 atom_info%r(2, natom), &
239 atom_info%r(3, natom), &
240 adjustl(trim(id2str(atom_info%id_molname(natom)))), &
241 adjustr(trim(id2str(atom_info%id_element(natom))))
242 END IF
243 atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "angstrom")
244 atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "angstrom")
245 atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "angstrom")
246 CASE ("TER")
247 inum_mol = inum_mol + 1
248 WRITE (unit=root_mol_name, fmt='(A3,I0)') "MOL", inum_mol
249 CASE ("REMARK")
250 IF (iw > 0) WRITE (unit=iw, fmt=*) trim(line)
251 CASE ("END")
252 EXIT
253 CASE DEFAULT
254 END SELECT
255 END DO
256 CALL parser_release(parser)
257
258 CALL reallocate(atom_info%id_molname, 1, natom)
259 CALL reallocate(atom_info%id_resname, 1, natom)
260 CALL reallocate(atom_info%resid, 1, natom)
261 CALL reallocate(atom_info%id_atmname, 1, natom)
262 CALL reallocate(atom_info%r, 1, 3, 1, natom)
263 CALL reallocate(atom_info%atm_mass, 1, natom)
264 CALL reallocate(atom_info%atm_charge, 1, natom)
265 CALL reallocate(atom_info%occup, 1, natom)
266 CALL reallocate(atom_info%beta, 1, natom)
267 CALL reallocate(atom_info%id_element, 1, natom)
268
269 IF (topology%conn_type /= do_conn_user) THEN
270 IF (.NOT. topology%para_res) atom_info%resid(:) = 1
271 END IF
272
273 IF (iw > 0) THEN
274 WRITE (unit=iw, fmt="(T2,A)") &
275 "END of PDB data read from file "//trim(topology%coord_file_name)
276 END IF
277
278 topology%natoms = natom
279 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
280 "PRINT%TOPOLOGY_INFO/PDB_INFO")
281 CALL timestop(handle)
282
283 END SUBROUTINE read_coordinate_pdb
284
285! **************************************************************************************************
286!> \brief ...
287!> \param file_unit ...
288!> \param topology ...
289!> \param subsys_section ...
290! **************************************************************************************************
291 SUBROUTINE write_coordinate_pdb(file_unit, topology, subsys_section)
292
293 INTEGER, INTENT(IN) :: file_unit
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 NULLIFY (logger)
313 logger => cp_get_default_logger()
314 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
315 extension=".subsysLog")
316 print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PDB")
317 CALL timeset(routinen, handle)
318
319 CALL section_vals_val_get(print_key, "CHARGE_OCCUP", l_val=charge_occup)
320 CALL section_vals_val_get(print_key, "CHARGE_BETA", l_val=charge_beta)
321 CALL section_vals_val_get(print_key, "CHARGE_EXTENDED", l_val=charge_extended)
322 i = count([charge_occup, charge_beta, charge_extended])
323 IF (i > 1) THEN
324 cpabort("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
325 END IF
326
327 atom_info => topology%atom_info
328 record = cp_print_key_generate_filename(logger, print_key, &
329 extension=".pdb", &
330 my_local=.false.)
331
332 IF (iw > 0) WRITE (unit=iw, fmt=*) " Writing out PDB file ", trim(record)
333
334 ! Write file header
335 CALL m_timestamp(timestamp)
336 WRITE (unit=file_unit, fmt="(A6,T11,A)") &
337 "TITLE ", "PDB file created by "//trim(cp2k_version)//" (revision "//trim(compile_revision)//")", &
338 "AUTHOR", trim(r_user_name)//"@"//trim(r_host_name)//" "//timestamp(:19)
339 ! Write cell information
340 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 "CRYST1", abc(1:3)*angstrom, angle_alpha, angle_beta, angle_gamma
343
344 natom = topology%natoms
345 idres = 0
346 id1 = 0
347 id2 = 0
348
349 DO i = 1, natom
350
351 IF (topology%para_res) THEN
352 idres = atom_info%resid(i)
353 ELSE
354 IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i))) THEN
355 idres = idres + 1
356 id1 = atom_info%map_mol_num(i)
357 id2 = atom_info%map_mol_typ(i)
358 END IF
359 END IF
360
361 line = ""
362 my_tag1 = id2str(atom_info%id_atmname(i)); ldum = qmmm_ff_precond_only_qm(my_tag1)
363 my_tag2 = id2str(atom_info%id_resname(i)); ldum = qmmm_ff_precond_only_qm(my_tag2)
364 my_tag3 = id2str(atom_info%id_molname(i)); ldum = qmmm_ff_precond_only_qm(my_tag3)
365 my_tag4 = id2str(atom_info%id_element(i)); ldum = qmmm_ff_precond_only_qm(my_tag4)
366
367 WRITE (unit=line(1:6), fmt="(A6)") "ATOM "
368 WRITE (unit=line(7:11), fmt="(I5)") modulo(i, 100000)
369 WRITE (unit=line(13:16), fmt="(A4)") adjustl(my_tag1(1:4))
370 WRITE (unit=line(18:20), fmt="(A3)") trim(my_tag2)
371 WRITE (unit=line(23:26), fmt="(I4)") modulo(idres, 10000)
372 WRITE (unit=line(31:54), fmt="(3F8.3)") atom_info%r(1:3, i)*angstrom
373 IF (ASSOCIATED(atom_info%occup)) THEN
374 WRITE (unit=line(55:60), fmt="(F6.2)") atom_info%occup(i)
375 ELSE
376 WRITE (unit=line(55:60), fmt="(F6.2)") 0.0_dp
377 END IF
378 IF (ASSOCIATED(atom_info%beta)) THEN
379 WRITE (unit=line(61:66), fmt="(F6.2)") atom_info%beta(i)
380 ELSE
381 WRITE (unit=line(61:66), fmt="(F6.2)") 0.0_dp
382 END IF
383 IF (ASSOCIATED(atom_info%atm_charge)) THEN
384 IF (any([charge_occup, charge_beta, charge_extended]) .AND. &
385 (atom_info%atm_charge(i) == -huge(0.0_dp))) THEN
386 cpabort("No atomic charges found yet (after the topology setup)")
387 END IF
388 IF (charge_occup) THEN
389 WRITE (unit=line(55:60), fmt="(F6.2)") atom_info%atm_charge(i)
390 ELSE IF (charge_beta) THEN
391 WRITE (unit=line(61:66), fmt="(F6.2)") atom_info%atm_charge(i)
392 ELSE IF (charge_extended) THEN
393 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 WRITE (unit=line(73:76), fmt="(A4)") adjustl(my_tag3)
399 WRITE (unit=line(77:78), fmt="(A2)") trim(my_tag4)
400 WRITE (unit=file_unit, fmt="(A)") trim(line)
401 END DO
402 WRITE (unit=file_unit, fmt="(A3)") "END"
403
404 IF (iw > 0) WRITE (unit=iw, fmt=*) " Exiting "//routinen
405
406 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
407 "PRINT%TOPOLOGY_INFO/PDB_INFO")
408
409 CALL timestop(handle)
410
411 END SUBROUTINE write_coordinate_pdb
412
413END MODULE topology_pdb
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:22
character(len=default_string_length), public r_host_name
Definition cp2k_info.F:74
character(len= *), parameter, public compile_revision
Definition cp2k_info.F:45
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:49
character(len=default_string_length), public r_user_name
Definition cp2k_info.F:74
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1222
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_conn_user
objects that represent the structure of input sections and the data contained in an input section
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
integer, parameter, public timestamp_length
Definition machine.F:46
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:381
Utility routines for the memory handling.
Interface to the message passing library MPI.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
generates a unique id number for a string (str2id) that can be used two compare two strings....
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Handles PDB files.
subroutine, public write_coordinate_pdb(file_unit, topology, subsys_section)
...
subroutine, public read_coordinate_pdb(topology, para_env, subsys_section)
...
Control for reading in different topologies and coordinates.
Definition topology.F:13
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment