Line data Source code
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 : ! **************************************************************************************************
48 : MODULE topology_pdb
49 : USE cell_types, ONLY: get_cell
50 : USE cp2k_info, ONLY: compile_revision,&
51 : cp2k_version,&
52 : r_datx,&
53 : r_host_name,&
54 : r_user_name
55 : USE cp_log_handling, ONLY: cp_get_default_logger,&
56 : cp_logger_type
57 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
58 : cp_print_key_generate_filename,&
59 : cp_print_key_unit_nr
60 : USE cp_parser_methods, ONLY: parser_get_next_line
61 : USE cp_parser_types, ONLY: cp_parser_type,&
62 : parser_create,&
63 : parser_release
64 : USE cp_units, ONLY: cp_unit_to_cp2k
65 : USE input_constants, ONLY: do_conn_user
66 : USE input_section_types, ONLY: section_get_rval,&
67 : section_vals_get_subs_vals,&
68 : section_vals_type,&
69 : section_vals_val_get
70 : USE kinds, ONLY: default_path_length,&
71 : default_string_length,&
72 : dp
73 : USE memory_utilities, ONLY: reallocate
74 : USE message_passing, ONLY: mp_para_env_type
75 : USE physcon, ONLY: angstrom
76 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
77 : USE string_table, ONLY: id2str,&
78 : s2s,&
79 : str2id
80 : USE topology_types, ONLY: atom_info_type,&
81 : topology_parameters_type
82 : #include "./base/base_uses.f90"
83 :
84 : IMPLICIT NONE
85 :
86 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_pdb'
87 :
88 : PRIVATE
89 : PUBLIC :: read_coordinate_pdb, write_coordinate_pdb
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 1604 : SUBROUTINE read_coordinate_pdb(topology, para_env, subsys_section)
102 : TYPE(topology_parameters_type) :: topology
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 802 : NULLIFY (logger)
120 1604 : logger => cp_get_default_logger()
121 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
122 802 : extension=".subsysLog")
123 802 : CALL timeset(routineN, handle)
124 :
125 802 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
126 802 : atom_info => topology%atom_info
127 802 : CALL reallocate(atom_info%id_molname, 1, nblock)
128 802 : CALL reallocate(atom_info%id_resname, 1, nblock)
129 802 : CALL reallocate(atom_info%resid, 1, nblock)
130 802 : CALL reallocate(atom_info%id_atmname, 1, nblock)
131 802 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
132 802 : CALL reallocate(atom_info%atm_mass, 1, nblock)
133 802 : CALL reallocate(atom_info%atm_charge, 1, nblock)
134 802 : CALL reallocate(atom_info%occup, 1, nblock)
135 802 : CALL reallocate(atom_info%beta, 1, nblock)
136 802 : CALL reallocate(atom_info%id_element, 1, nblock)
137 :
138 802 : IF (iw > 0) THEN
139 : WRITE (UNIT=iw, FMT="(T2,A)") &
140 1 : "BEGIN of PDB data read from file "//TRIM(topology%coord_file_name)
141 : END IF
142 :
143 802 : id0 = str2id(s2s(""))
144 802 : topology%molname_generated = .FALSE.
145 :
146 802 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
147 :
148 802 : natom = 0
149 802 : inum_mol = 1
150 802 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
151 : DO
152 440683 : line = ""
153 440683 : CALL parser_get_next_line(parser, 1, at_end=my_end)
154 440683 : IF (my_end) EXIT
155 440183 : line = parser%input_line(1:default_path_length)
156 440183 : record = line(1:6)
157 : record = TRIM(record)
158 :
159 440183 : IF ((record == "ATOM") .OR. (record == "HETATM")) THEN
160 397989 : natom = natom + 1
161 397989 : topology%natoms = natom
162 397989 : IF (natom > SIZE(atom_info%id_atmname)) THEN
163 500 : newsize = INT(pfactor*natom)
164 500 : CALL reallocate(atom_info%id_molname, 1, newsize)
165 500 : CALL reallocate(atom_info%id_resname, 1, newsize)
166 500 : CALL reallocate(atom_info%resid, 1, newsize)
167 500 : CALL reallocate(atom_info%id_atmname, 1, newsize)
168 500 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
169 500 : CALL reallocate(atom_info%atm_mass, 1, newsize)
170 500 : CALL reallocate(atom_info%atm_charge, 1, newsize)
171 500 : CALL reallocate(atom_info%occup, 1, newsize)
172 500 : CALL reallocate(atom_info%beta, 1, newsize)
173 500 : CALL reallocate(atom_info%id_element, 1, newsize)
174 : END IF
175 : END IF
176 :
177 500 : SELECT CASE (record)
178 : CASE ("ATOM", "HETATM")
179 397989 : READ (UNIT=line(13:16), FMT=*) strtmp
180 397989 : atom_info%id_atmname(natom) = str2id(s2s(strtmp))
181 397989 : READ (UNIT=line(18:20), FMT=*, IOSTAT=istat) strtmp
182 397989 : IF (istat == 0) THEN
183 394823 : atom_info%id_resname(natom) = str2id(s2s(strtmp))
184 : ELSE
185 3166 : atom_info%id_resname(natom) = id0
186 : END IF
187 397989 : READ (UNIT=line(23:26), FMT=*, IOSTAT=istat) atom_info%resid(natom)
188 397989 : READ (UNIT=line(31:38), FMT=*, IOSTAT=istat) atom_info%r(1, natom)
189 397989 : READ (UNIT=line(39:46), FMT=*, IOSTAT=istat) atom_info%r(2, natom)
190 397989 : READ (UNIT=line(47:54), FMT=*, IOSTAT=istat) atom_info%r(3, natom)
191 397989 : READ (UNIT=line(55:60), FMT=*, IOSTAT=istat) atom_info%occup(natom)
192 397989 : READ (UNIT=line(61:66), FMT=*, IOSTAT=istat) atom_info%beta(natom)
193 397989 : READ (UNIT=line(73:76), FMT=*, IOSTAT=istat) strtmp
194 397989 : IF (istat == 0) THEN
195 239544 : atom_info%id_molname(natom) = str2id(s2s(strtmp))
196 : ELSE
197 158445 : atom_info%id_molname(natom) = str2id(s2s(root_mol_name))
198 158445 : topology%molname_generated = .TRUE.
199 : END IF
200 397989 : READ (UNIT=line(77:78), FMT=*, IOSTAT=istat) strtmp
201 397989 : IF (istat == 0) THEN
202 160389 : atom_info%id_element(natom) = str2id(s2s(strtmp))
203 : ELSE
204 237600 : atom_info%id_element(natom) = id0
205 : END IF
206 397989 : atom_info%atm_mass(natom) = 0.0_dp
207 397989 : atom_info%atm_charge(natom) = -HUGE(0.0_dp)
208 397989 : IF (topology%charge_occup) atom_info%atm_charge(natom) = atom_info%occup(natom)
209 397989 : IF (topology%charge_beta) atom_info%atm_charge(natom) = atom_info%beta(natom)
210 397989 : IF (topology%charge_extended) THEN
211 3188 : READ (UNIT=line(81:), FMT=*, IOSTAT=istat) atom_info%atm_charge(natom)
212 : END IF
213 :
214 397989 : IF (atom_info%id_element(natom) == id0) THEN
215 : ! Element is assigned on the basis of the atm_name
216 237600 : topology%aa_element = .TRUE.
217 237600 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
218 : END IF
219 :
220 397989 : 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 6 : record, natom, &
223 6 : TRIM(id2str(atom_info%id_atmname(natom))), &
224 6 : TRIM(id2str(atom_info%id_resname(natom))), &
225 6 : atom_info%resid(natom), &
226 6 : atom_info%r(1, natom), &
227 6 : atom_info%r(2, natom), &
228 6 : atom_info%r(3, natom), &
229 6 : ADJUSTL(TRIM(id2str(atom_info%id_molname(natom)))), &
230 12 : ADJUSTR(TRIM(id2str(atom_info%id_element(natom))))
231 : END IF
232 397989 : atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "angstrom")
233 397989 : atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "angstrom")
234 397989 : atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "angstrom")
235 : CASE ("TER")
236 41506 : inum_mol = inum_mol + 1
237 41506 : WRITE (UNIT=root_mol_name, FMT='(A3,I0)') "MOL", inum_mol
238 : CASE ("REMARK")
239 278 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) TRIM(line)
240 : CASE ("END")
241 440183 : EXIT
242 : CASE DEFAULT
243 : END SELECT
244 : END DO
245 802 : CALL parser_release(parser)
246 :
247 802 : CALL reallocate(atom_info%id_molname, 1, natom)
248 802 : CALL reallocate(atom_info%id_resname, 1, natom)
249 802 : CALL reallocate(atom_info%resid, 1, natom)
250 802 : CALL reallocate(atom_info%id_atmname, 1, natom)
251 802 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
252 802 : CALL reallocate(atom_info%atm_mass, 1, natom)
253 802 : CALL reallocate(atom_info%atm_charge, 1, natom)
254 802 : CALL reallocate(atom_info%occup, 1, natom)
255 802 : CALL reallocate(atom_info%beta, 1, natom)
256 802 : CALL reallocate(atom_info%id_element, 1, natom)
257 :
258 802 : IF (topology%conn_type /= do_conn_user) THEN
259 938 : IF (.NOT. topology%para_res) atom_info%resid(:) = 1
260 : END IF
261 :
262 802 : IF (iw > 0) THEN
263 : WRITE (UNIT=iw, FMT="(T2,A)") &
264 1 : "END of PDB data read from file "//TRIM(topology%coord_file_name)
265 : END IF
266 :
267 802 : topology%natoms = natom
268 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
269 802 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
270 802 : CALL timestop(handle)
271 :
272 2406 : END SUBROUTINE read_coordinate_pdb
273 :
274 : ! **************************************************************************************************
275 : !> \brief ...
276 : !> \param file_unit ...
277 : !> \param topology ...
278 : !> \param subsys_section ...
279 : ! **************************************************************************************************
280 318 : SUBROUTINE write_coordinate_pdb(file_unit, topology, subsys_section)
281 :
282 : INTEGER, INTENT(IN) :: file_unit
283 : TYPE(topology_parameters_type) :: topology
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 53 : NULLIFY (logger)
301 53 : logger => cp_get_default_logger()
302 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PDB_INFO", &
303 53 : extension=".subsysLog")
304 53 : print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PDB")
305 53 : CALL timeset(routineN, handle)
306 :
307 53 : CALL section_vals_val_get(print_key, "CHARGE_OCCUP", l_val=charge_occup)
308 53 : CALL section_vals_val_get(print_key, "CHARGE_BETA", l_val=charge_beta)
309 53 : CALL section_vals_val_get(print_key, "CHARGE_EXTENDED", l_val=charge_extended)
310 212 : i = COUNT((/charge_occup, charge_beta, charge_extended/))
311 53 : IF (i > 1) &
312 0 : CPABORT("Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
313 :
314 53 : atom_info => topology%atom_info
315 : record = cp_print_key_generate_filename(logger, print_key, &
316 : extension=".pdb", &
317 53 : my_local=.FALSE.)
318 :
319 53 : 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 53 : "TITLE ", "PDB file created by "//TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")", &
324 106 : "AUTHOR", TRIM(r_user_name)//"@"//TRIM(r_host_name)//" "//r_datx(1:19)
325 : ! Write cell information
326 53 : 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 212 : "CRYST1", abc(1:3)*angstrom, angle_alpha, angle_beta, angle_gamma
329 :
330 53 : natom = topology%natoms
331 53 : idres = 0
332 53 : id1 = 0
333 53 : id2 = 0
334 :
335 28984 : DO i = 1, natom
336 :
337 28931 : IF (topology%para_res) THEN
338 28235 : idres = atom_info%resid(i)
339 : ELSE
340 696 : IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i))) THEN
341 49 : idres = idres + 1
342 49 : id1 = atom_info%map_mol_num(i)
343 49 : id2 = atom_info%map_mol_typ(i)
344 : END IF
345 : END IF
346 :
347 28931 : line = ""
348 28931 : my_tag1 = id2str(atom_info%id_atmname(i)); ldum = qmmm_ff_precond_only_qm(my_tag1)
349 28931 : my_tag2 = id2str(atom_info%id_resname(i)); ldum = qmmm_ff_precond_only_qm(my_tag2)
350 28931 : my_tag3 = id2str(atom_info%id_molname(i)); ldum = qmmm_ff_precond_only_qm(my_tag3)
351 28931 : my_tag4 = id2str(atom_info%id_element(i)); ldum = qmmm_ff_precond_only_qm(my_tag4)
352 :
353 28931 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
354 28931 : WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(i, 100000)
355 28931 : WRITE (UNIT=line(13:16), FMT="(A4)") ADJUSTL(my_tag1(1:4))
356 28931 : WRITE (UNIT=line(18:20), FMT="(A3)") TRIM(my_tag2)
357 28931 : WRITE (UNIT=line(23:26), FMT="(I4)") MODULO(idres, 10000)
358 115724 : WRITE (UNIT=line(31:54), FMT="(3F8.3)") atom_info%r(1:3, i)*angstrom
359 28931 : IF (ASSOCIATED(atom_info%occup)) THEN
360 28652 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%occup(i)
361 : ELSE
362 279 : WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp
363 : END IF
364 28931 : IF (ASSOCIATED(atom_info%beta)) THEN
365 28652 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%beta(i)
366 : ELSE
367 279 : WRITE (UNIT=line(61:66), FMT="(F6.2)") 0.0_dp
368 : END IF
369 28931 : IF (ASSOCIATED(atom_info%atm_charge)) THEN
370 28931 : IF (ANY((/charge_occup, charge_beta, charge_extended/)) .AND. &
371 : (atom_info%atm_charge(i) == -HUGE(0.0_dp))) &
372 0 : CPABORT("No atomic charges found yet (after the topology setup)")
373 28931 : IF (charge_occup) THEN
374 0 : WRITE (UNIT=line(55:60), FMT="(F6.2)") atom_info%atm_charge(i)
375 28931 : ELSE IF (charge_beta) THEN
376 0 : WRITE (UNIT=line(61:66), FMT="(F6.2)") atom_info%atm_charge(i)
377 28931 : ELSE IF (charge_extended) THEN
378 0 : 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 28931 : WRITE (UNIT=line(73:76), FMT="(A4)") ADJUSTL(my_tag3)
384 28931 : WRITE (UNIT=line(77:78), FMT="(A2)") TRIM(my_tag4)
385 28984 : WRITE (UNIT=file_unit, FMT="(A)") TRIM(line)
386 : END DO
387 53 : WRITE (UNIT=file_unit, FMT="(A3)") "END"
388 :
389 53 : IF (iw > 0) WRITE (UNIT=iw, FMT=*) " Exiting "//routineN
390 :
391 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
392 53 : "PRINT%TOPOLOGY_INFO/PDB_INFO")
393 :
394 53 : CALL timestop(handle)
395 :
396 53 : END SUBROUTINE write_coordinate_pdb
397 :
398 : END MODULE topology_pdb
|
