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