(git:3add494)
topology_pdb.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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,&
51  cp2k_version,&
52  r_datx,&
53  r_host_name,&
56  cp_logger_type
61  USE cp_parser_types, ONLY: cp_parser_type,&
64  USE cp_units, ONLY: cp_unit_to_cp2k
65  USE input_constants, ONLY: do_conn_user
68  section_vals_type,&
70  USE kinds, ONLY: default_path_length,&
72  dp
73  USE memory_utilities, ONLY: reallocate
74  USE message_passing, ONLY: mp_para_env_type
75  USE physcon, ONLY: angstrom
77  USE string_table, ONLY: id2str,&
78  s2s,&
79  str2id
80  USE topology_types, ONLY: atom_info_type,&
82 #include "./base/base_uses.f90"
83 
84  IMPLICIT NONE
85 
86  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_pdb'
87 
88  PRIVATE
90 
91 CONTAINS
92 
93 ! **************************************************************************************************
94 !> \brief ...
95 !> \param topology ...
96 !> \param para_env ...
97 !> \param subsys_section ...
98 !> \par History
99 !> TLAINO 05.2004 - Added the TER option to use different non-bonded molecules
100 ! **************************************************************************************************
101  SUBROUTINE read_coordinate_pdb(topology, para_env, subsys_section)
103  TYPE(mp_para_env_type), POINTER :: para_env
104  TYPE(section_vals_type), POINTER :: subsys_section
105 
106  CHARACTER(len=*), PARAMETER :: routinen = 'read_coordinate_pdb'
107  INTEGER, PARAMETER :: nblock = 1000
108 
109  CHARACTER(LEN=default_path_length) :: line
110  CHARACTER(LEN=default_string_length) :: record, root_mol_name, strtmp
111  INTEGER :: handle, id0, inum_mol, istat, iw, natom, &
112  newsize
113  LOGICAL :: my_end
114  REAL(kind=dp) :: pfactor
115  TYPE(atom_info_type), POINTER :: atom_info
116  TYPE(cp_logger_type), POINTER :: logger
117  TYPE(cp_parser_type) :: parser
118 
119  NULLIFY (logger)
120  logger => cp_get_default_logger()
121  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
122  extension=".subsysLog")
123  CALL timeset(routinen, handle)
124 
125  pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
126  atom_info => topology%atom_info
127  CALL reallocate(atom_info%id_molname, 1, nblock)
128  CALL reallocate(atom_info%id_resname, 1, nblock)
129  CALL reallocate(atom_info%resid, 1, nblock)
130  CALL reallocate(atom_info%id_atmname, 1, nblock)
131  CALL reallocate(atom_info%r, 1, 3, 1, nblock)
132  CALL reallocate(atom_info%atm_mass, 1, nblock)
133  CALL reallocate(atom_info%atm_charge, 1, nblock)
134  CALL reallocate(atom_info%occup, 1, nblock)
135  CALL reallocate(atom_info%beta, 1, nblock)
136  CALL reallocate(atom_info%id_element, 1, nblock)
137 
138  IF (iw > 0) THEN
139  WRITE (unit=iw, fmt="(T2,A)") &
140  "BEGIN of PDB data read from file "//trim(topology%coord_file_name)
141  END IF
142 
143  id0 = str2id(s2s(""))
144  topology%molname_generated = .false.
145 
146  CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
147 
148  natom = 0
149  inum_mol = 1
150  WRITE (unit=root_mol_name, fmt='(A3,I0)') "MOL", inum_mol
151  DO
152  line = ""
153  CALL parser_get_next_line(parser, 1, at_end=my_end)
154  IF (my_end) EXIT
155  line = parser%input_line(1:default_path_length)
156  record = line(1:6)
157  record = trim(record)
158 
159  IF ((record == "ATOM") .OR. (record == "HETATM")) THEN
160  natom = natom + 1
161  topology%natoms = natom
162  IF (natom > SIZE(atom_info%id_atmname)) THEN
163  newsize = int(pfactor*natom)
164  CALL reallocate(atom_info%id_molname, 1, newsize)
165  CALL reallocate(atom_info%id_resname, 1, newsize)
166  CALL reallocate(atom_info%resid, 1, newsize)
167  CALL reallocate(atom_info%id_atmname, 1, newsize)
168  CALL reallocate(atom_info%r, 1, 3, 1, newsize)
169  CALL reallocate(atom_info%atm_mass, 1, newsize)
170  CALL reallocate(atom_info%atm_charge, 1, newsize)
171  CALL reallocate(atom_info%occup, 1, newsize)
172  CALL reallocate(atom_info%beta, 1, newsize)
173  CALL reallocate(atom_info%id_element, 1, newsize)
174  END IF
175  END IF
176 
177  SELECT CASE (record)
178  CASE ("ATOM", "HETATM")
179  READ (unit=line(13:16), fmt=*) strtmp
180  atom_info%id_atmname(natom) = str2id(s2s(strtmp))
181  READ (unit=line(18:20), fmt=*, iostat=istat) strtmp
182  IF (istat == 0) THEN
183  atom_info%id_resname(natom) = str2id(s2s(strtmp))
184  ELSE
185  atom_info%id_resname(natom) = id0
186  END IF
187  READ (unit=line(23:26), fmt=*, iostat=istat) atom_info%resid(natom)
188  READ (unit=line(31:38), fmt=*, iostat=istat) atom_info%r(1, natom)
189  READ (unit=line(39:46), fmt=*, iostat=istat) atom_info%r(2, natom)
190  READ (unit=line(47:54), fmt=*, iostat=istat) atom_info%r(3, natom)
191  READ (unit=line(55:60), fmt=*, iostat=istat) atom_info%occup(natom)
192  READ (unit=line(61:66), fmt=*, iostat=istat) atom_info%beta(natom)
193  READ (unit=line(73:76), fmt=*, iostat=istat) strtmp
194  IF (istat == 0) THEN
195  atom_info%id_molname(natom) = str2id(s2s(strtmp))
196  ELSE
197  atom_info%id_molname(natom) = str2id(s2s(root_mol_name))
198  topology%molname_generated = .true.
199  END IF
200  READ (unit=line(77:78), fmt=*, iostat=istat) strtmp
201  IF (istat == 0) THEN
202  atom_info%id_element(natom) = str2id(s2s(strtmp))
203  ELSE
204  atom_info%id_element(natom) = id0
205  END IF
206  atom_info%atm_mass(natom) = 0.0_dp
207  atom_info%atm_charge(natom) = -huge(0.0_dp)
208  IF (topology%charge_occup) atom_info%atm_charge(natom) = atom_info%occup(natom)
209  IF (topology%charge_beta) atom_info%atm_charge(natom) = atom_info%beta(natom)
210  IF (topology%charge_extended) THEN
211  READ (unit=line(81:), fmt=*, iostat=istat) atom_info%atm_charge(natom)
212  END IF
213 
214  IF (atom_info%id_element(natom) == id0) THEN
215  ! Element is assigned on the basis of the atm_name
216  topology%aa_element = .true.
217  atom_info%id_element(natom) = atom_info%id_atmname(natom)
218  END IF
219 
220  IF (iw > 0) THEN
221  WRITE (unit=iw, fmt="(A6,I5,T13,A4,T18,A3,T23,I4,T31,3F8.3,T73,A4,T77,A2)") &
222  record, natom, &
223  trim(id2str(atom_info%id_atmname(natom))), &
224  trim(id2str(atom_info%id_resname(natom))), &
225  atom_info%resid(natom), &
226  atom_info%r(1, natom), &
227  atom_info%r(2, natom), &
228  atom_info%r(3, natom), &
229  adjustl(trim(id2str(atom_info%id_molname(natom)))), &
230  adjustr(trim(id2str(atom_info%id_element(natom))))
231  END IF
232  atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "angstrom")
233  atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "angstrom")
234  atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "angstrom")
235  CASE ("TER")
236  inum_mol = inum_mol + 1
237  WRITE (unit=root_mol_name, fmt='(A3,I0)') "MOL", inum_mol
238  CASE ("REMARK")
239  IF (iw > 0) WRITE (unit=iw, fmt=*) trim(line)
240  CASE ("END")
241  EXIT
242  CASE DEFAULT
243  END SELECT
244  END DO
245  CALL parser_release(parser)
246 
247  CALL reallocate(atom_info%id_molname, 1, natom)
248  CALL reallocate(atom_info%id_resname, 1, natom)
249  CALL reallocate(atom_info%resid, 1, natom)
250  CALL reallocate(atom_info%id_atmname, 1, natom)
251  CALL reallocate(atom_info%r, 1, 3, 1, natom)
252  CALL reallocate(atom_info%atm_mass, 1, natom)
253  CALL reallocate(atom_info%atm_charge, 1, natom)
254  CALL reallocate(atom_info%occup, 1, natom)
255  CALL reallocate(atom_info%beta, 1, natom)
256  CALL reallocate(atom_info%id_element, 1, natom)
257 
258  IF (topology%conn_type /= do_conn_user) THEN
259  IF (.NOT. topology%para_res) atom_info%resid(:) = 1
260  END IF
261 
262  IF (iw > 0) THEN
263  WRITE (unit=iw, fmt="(T2,A)") &
264  "END of PDB data read from file "//trim(topology%coord_file_name)
265  END IF
266 
267  topology%natoms = natom
268  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
269  "PRINT%TOPOLOGY_INFO/PDB_INFO")
270  CALL timestop(handle)
271 
272  END SUBROUTINE read_coordinate_pdb
273 
274 ! **************************************************************************************************
275 !> \brief ...
276 !> \param file_unit ...
277 !> \param topology ...
278 !> \param subsys_section ...
279 ! **************************************************************************************************
280  SUBROUTINE write_coordinate_pdb(file_unit, topology, subsys_section)
281 
282  INTEGER, INTENT(IN) :: file_unit
284  TYPE(section_vals_type), POINTER :: subsys_section
285 
286  CHARACTER(len=*), PARAMETER :: routinen = 'write_coordinate_pdb'
287 
288  CHARACTER(LEN=120) :: line
289  CHARACTER(LEN=default_string_length) :: my_tag1, my_tag2, my_tag3, my_tag4, &
290  record
291  INTEGER :: handle, i, id1, id2, idres, iw, natom
292  LOGICAL :: charge_beta, charge_extended, &
293  charge_occup, ldum
294  REAL(kind=dp) :: angle_alpha, angle_beta, angle_gamma
295  REAL(kind=dp), DIMENSION(3) :: abc
296  TYPE(atom_info_type), POINTER :: atom_info
297  TYPE(cp_logger_type), POINTER :: logger
298  TYPE(section_vals_type), POINTER :: print_key
299 
300  NULLIFY (logger)
301  logger => cp_get_default_logger()
302  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
303  extension=".subsysLog")
304  print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PDB")
305  CALL timeset(routinen, handle)
306 
307  CALL section_vals_val_get(print_key, "CHARGE_OCCUP", l_val=charge_occup)
308  CALL section_vals_val_get(print_key, "CHARGE_BETA", l_val=charge_beta)
309  CALL section_vals_val_get(print_key, "CHARGE_EXTENDED", l_val=charge_extended)
310  i = count((/charge_occup, charge_beta, charge_extended/))
311  IF (i > 1) &
312  cpabort("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
313 
314  atom_info => topology%atom_info
315  record = cp_print_key_generate_filename(logger, print_key, &
316  extension=".pdb", &
317  my_local=.false.)
318 
319  IF (iw > 0) WRITE (unit=iw, fmt=*) " Writing out PDB file ", trim(record)
320 
321  ! Write file header
322  WRITE (unit=file_unit, fmt="(A6,T11,A)") &
323  "TITLE ", "PDB file created by "//trim(cp2k_version)//" (revision "//trim(compile_revision)//")", &
324  "AUTHOR", trim(r_user_name)//"@"//trim(r_host_name)//" "//r_datx(1:19)
325  ! Write cell information
326  CALL get_cell(cell=topology%cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc)
327  WRITE (unit=file_unit, fmt="(A6,3F9.3,3F7.2)") &
328  "CRYST1", abc(1:3)*angstrom, angle_alpha, angle_beta, angle_gamma
329 
330  natom = topology%natoms
331  idres = 0
332  id1 = 0
333  id2 = 0
334 
335  DO i = 1, natom
336 
337  IF (topology%para_res) THEN
338  idres = atom_info%resid(i)
339  ELSE
340  IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i))) THEN
341  idres = idres + 1
342  id1 = atom_info%map_mol_num(i)
343  id2 = atom_info%map_mol_typ(i)
344  END IF
345  END IF
346 
347  line = ""
348  my_tag1 = id2str(atom_info%id_atmname(i)); ldum = qmmm_ff_precond_only_qm(my_tag1)
349  my_tag2 = id2str(atom_info%id_resname(i)); ldum = qmmm_ff_precond_only_qm(my_tag2)
350  my_tag3 = id2str(atom_info%id_molname(i)); ldum = qmmm_ff_precond_only_qm(my_tag3)
351  my_tag4 = id2str(atom_info%id_element(i)); ldum = qmmm_ff_precond_only_qm(my_tag4)
352 
353  WRITE (unit=line(1:6), fmt="(A6)") "ATOM "
354  WRITE (unit=line(7:11), fmt="(I5)") modulo(i, 100000)
355  WRITE (unit=line(13:16), fmt="(A4)") adjustl(my_tag1(1:4))
356  WRITE (unit=line(18:20), fmt="(A3)") trim(my_tag2)
357  WRITE (unit=line(23:26), fmt="(I4)") modulo(idres, 10000)
358  WRITE (unit=line(31:54), fmt="(3F8.3)") atom_info%r(1:3, i)*angstrom
359  IF (ASSOCIATED(atom_info%occup)) THEN
360  WRITE (unit=line(55:60), fmt="(F6.2)") atom_info%occup(i)
361  ELSE
362  WRITE (unit=line(55:60), fmt="(F6.2)") 0.0_dp
363  END IF
364  IF (ASSOCIATED(atom_info%beta)) THEN
365  WRITE (unit=line(61:66), fmt="(F6.2)") atom_info%beta(i)
366  ELSE
367  WRITE (unit=line(61:66), fmt="(F6.2)") 0.0_dp
368  END IF
369  IF (ASSOCIATED(atom_info%atm_charge)) THEN
370  IF (any((/charge_occup, charge_beta, charge_extended/)) .AND. &
371  (atom_info%atm_charge(i) == -huge(0.0_dp))) &
372  cpabort("No atomic charges found yet (after the topology setup)")
373  IF (charge_occup) THEN
374  WRITE (unit=line(55:60), fmt="(F6.2)") atom_info%atm_charge(i)
375  ELSE IF (charge_beta) THEN
376  WRITE (unit=line(61:66), fmt="(F6.2)") atom_info%atm_charge(i)
377  ELSE IF (charge_extended) THEN
378  WRITE (unit=line(81:), fmt="(F20.16)") atom_info%atm_charge(i)
379  ELSE
380  ! Write no atomic charge
381  END IF
382  END IF
383  WRITE (unit=line(73:76), fmt="(A4)") adjustl(my_tag3)
384  WRITE (unit=line(77:78), fmt="(A2)") trim(my_tag4)
385  WRITE (unit=file_unit, fmt="(A)") trim(line)
386  END DO
387  WRITE (unit=file_unit, fmt="(A3)") "END"
388 
389  IF (iw > 0) WRITE (unit=iw, fmt=*) " Exiting "//routinen
390 
391  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
392  "PRINT%TOPOLOGY_INFO/PDB_INFO")
393 
394  CALL timestop(handle)
395 
396  END SUBROUTINE write_coordinate_pdb
397 
398 END 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....
Definition: grid_common.h:117
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:195
some minimal info about CP2K, including its version and license
Definition: cp2k_info.F:16
character(len=default_string_length), public r_host_name
Definition: cp2k_info.F:66
character(len= *), parameter, public compile_revision
Definition: cp2k_info.F:37
character(len= *), parameter, public cp2k_version
Definition: cp2k_info.F:40
character(len=default_string_length), public r_user_name
Definition: cp2k_info.F:66
character(len=26), public r_datx
Definition: cp2k_info.F:64
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:1150
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
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 ...
Definition: qmmm_ff_fist.F:39
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition: string_table.F:72
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition: string_table.F:115
Handles PDB files.
Definition: topology_pdb.F:48
subroutine, public write_coordinate_pdb(file_unit, topology, subsys_section)
...
Definition: topology_pdb.F:281
subroutine, public read_coordinate_pdb(topology, para_env, subsys_section)
...
Definition: topology_pdb.F:102
Control for reading in different topologies and coordinates.
Definition: topology.F:13