107 CHARACTER(len=*),
PARAMETER :: routinen =
'read_coordinate_pdb'
108 INTEGER,
PARAMETER :: nblock = 1000
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, &
115 REAL(kind=
dp) :: pfactor
123 extension=
".subsysLog")
124 CALL timeset(routinen, handle)
126 pfactor =
section_get_rval(subsys_section,
"TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
128 CALL reallocate(atom_info%id_molname, 1, nblock)
129 CALL reallocate(atom_info%id_resname, 1, nblock)
131 CALL reallocate(atom_info%id_atmname, 1, nblock)
133 CALL reallocate(atom_info%atm_mass, 1, nblock)
134 CALL reallocate(atom_info%atm_charge, 1, nblock)
137 CALL reallocate(atom_info%id_element, 1, nblock)
140 WRITE (unit=iw, fmt=
"(T2,A)") &
141 "BEGIN of PDB data read from file "//trim(
topology%coord_file_name)
145 topology%molname_generated = .false.
151 WRITE (unit=root_mol_name, fmt=
'(A3,I0)')
"MOL", inum_mol
158 record = trim(record)
160 IF ((record ==
"ATOM") .OR. (record ==
"HETATM"))
THEN
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)
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)
174 CALL reallocate(atom_info%id_element, 1, newsize)
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
184 atom_info%id_resname(natom) =
str2id(
s2s(strtmp))
186 atom_info%id_resname(natom) = id0
190 READ (unit=line(23:26), fmt=*, iostat=istat) atom_info%resid(natom)
192 READ (unit=line(31:38), fmt=*, iostat=istat) atom_info%r(1, natom)
194 READ (unit=line(39:46), fmt=*, iostat=istat) atom_info%r(2, natom)
196 READ (unit=line(47:54), fmt=*, iostat=istat) atom_info%r(3, natom)
198 READ (unit=line(55:60), fmt=*, iostat=istat) atom_info%occup(natom)
200 READ (unit=line(61:66), fmt=*, iostat=istat) atom_info%beta(natom)
202 READ (unit=line(73:76), fmt=*, iostat=istat) strtmp
204 atom_info%id_molname(natom) =
str2id(
s2s(strtmp))
206 atom_info%id_molname(natom) =
str2id(
s2s(root_mol_name))
209 READ (unit=line(77:78), fmt=*, iostat=istat) strtmp
211 atom_info%id_element(natom) =
str2id(
s2s(strtmp))
213 atom_info%id_element(natom) = id0
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)
220 IF (len_trim(line) > 80)
THEN
221 READ (unit=line(81:), fmt=*) atom_info%atm_charge(natom)
225 IF (atom_info%id_element(natom) == id0)
THEN
228 atom_info%id_element(natom) = atom_info%id_atmname(natom)
232 WRITE (unit=iw, fmt=
"(A6,I5,T13,A4,T18,A3,T23,I4,T31,3F8.3,T73,A4,T77,A2)") &
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))))
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")
247 inum_mol = inum_mol + 1
248 WRITE (unit=root_mol_name, fmt=
'(A3,I0)')
"MOL", inum_mol
250 IF (iw > 0)
WRITE (unit=iw, fmt=*) trim(line)
258 CALL reallocate(atom_info%id_molname, 1, natom)
259 CALL reallocate(atom_info%id_resname, 1, natom)
261 CALL reallocate(atom_info%id_atmname, 1, natom)
264 CALL reallocate(atom_info%atm_charge, 1, natom)
267 CALL reallocate(atom_info%id_element, 1, natom)
270 IF (.NOT.
topology%para_res) atom_info%resid(:) = 1
274 WRITE (unit=iw, fmt=
"(T2,A)") &
275 "END of PDB data read from file "//trim(
topology%coord_file_name)
280 "PRINT%TOPOLOGY_INFO/PDB_INFO")
281 CALL timestop(handle)
293 INTEGER,
INTENT(IN) :: file_unit
297 CHARACTER(len=*),
PARAMETER :: routinen =
'write_coordinate_pdb'
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, &
306 REAL(kind=
dp) :: angle_alpha, angle_beta, angle_gamma
307 REAL(kind=
dp),
DIMENSION(3) :: abc
315 extension=
".subsysLog")
317 CALL timeset(routinen, handle)
322 i = count([charge_occup, charge_beta, charge_extended])
324 cpabort(
"Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
332 IF (iw > 0)
WRITE (unit=iw, fmt=*)
" Writing out PDB file ", trim(record)
336 WRITE (unit=file_unit, fmt=
"(A6,T11,A)") &
341 WRITE (unit=file_unit, fmt=
"(A6,3F9.3,3F7.2)") &
342 "CRYST1", abc(1:3)*
angstrom, angle_alpha, angle_beta, angle_gamma
352 idres = atom_info%resid(i)
354 IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i)))
THEN
356 id1 = atom_info%map_mol_num(i)
357 id2 = atom_info%map_mol_typ(i)
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)
376 WRITE (unit=line(55:60), fmt=
"(F6.2)") 0.0_dp
378 IF (
ASSOCIATED(atom_info%beta))
THEN
379 WRITE (unit=line(61:66), fmt=
"(F6.2)") atom_info%beta(i)
381 WRITE (unit=line(61:66), fmt=
"(F6.2)") 0.0_dp
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)")
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)
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)
402 WRITE (unit=file_unit, fmt=
"(A3)")
"END"
404 IF (iw > 0)
WRITE (unit=iw, fmt=*)
" Exiting "//routinen
407 "PRINT%TOPOLOGY_INFO/PDB_INFO")
409 CALL timestop(handle)