82 #include "./base/base_uses.f90"
86 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'topology_pdb'
103 TYPE(mp_para_env_type),
POINTER :: para_env
104 TYPE(section_vals_type),
POINTER :: subsys_section
106 CHARACTER(len=*),
PARAMETER :: routinen =
'read_coordinate_pdb'
107 INTEGER,
PARAMETER :: nblock = 1000
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, &
114 REAL(kind=
dp) :: pfactor
116 TYPE(cp_logger_type),
POINTER :: logger
117 TYPE(cp_parser_type) :: parser
122 extension=
".subsysLog")
123 CALL timeset(routinen, handle)
125 pfactor =
section_get_rval(subsys_section,
"TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
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)
139 WRITE (unit=iw, fmt=
"(T2,A)") &
140 "BEGIN of PDB data read from file "//trim(
topology%coord_file_name)
144 topology%molname_generated = .false.
150 WRITE (unit=root_mol_name, fmt=
'(A3,I0)')
"MOL", inum_mol
157 record = trim(record)
159 IF ((record ==
"ATOM") .OR. (record ==
"HETATM"))
THEN
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)
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
183 atom_info%id_resname(natom) =
str2id(
s2s(strtmp))
185 atom_info%id_resname(natom) = id0
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
195 atom_info%id_molname(natom) =
str2id(
s2s(strtmp))
197 atom_info%id_molname(natom) =
str2id(
s2s(root_mol_name))
200 READ (unit=line(77:78), fmt=*, iostat=istat) strtmp
202 atom_info%id_element(natom) =
str2id(
s2s(strtmp))
204 atom_info%id_element(natom) = id0
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)
211 READ (unit=line(81:), fmt=*, iostat=istat) atom_info%atm_charge(natom)
214 IF (atom_info%id_element(natom) == id0)
THEN
217 atom_info%id_element(natom) = atom_info%id_atmname(natom)
221 WRITE (unit=iw, fmt=
"(A6,I5,T13,A4,T18,A3,T23,I4,T31,3F8.3,T73,A4,T77,A2)") &
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))))
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")
236 inum_mol = inum_mol + 1
237 WRITE (unit=root_mol_name, fmt=
'(A3,I0)')
"MOL", inum_mol
239 IF (iw > 0)
WRITE (unit=iw, fmt=*) trim(line)
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)
259 IF (.NOT.
topology%para_res) atom_info%resid(:) = 1
263 WRITE (unit=iw, fmt=
"(T2,A)") &
264 "END of PDB data read from file "//trim(
topology%coord_file_name)
269 "PRINT%TOPOLOGY_INFO/PDB_INFO")
270 CALL timestop(handle)
282 INTEGER,
INTENT(IN) :: file_unit
284 TYPE(section_vals_type),
POINTER :: subsys_section
286 CHARACTER(len=*),
PARAMETER :: routinen =
'write_coordinate_pdb'
288 CHARACTER(LEN=120) :: line
289 CHARACTER(LEN=default_string_length) :: my_tag1, my_tag2, my_tag3, my_tag4, &
291 INTEGER :: handle, i, id1, id2, idres, iw, natom
292 LOGICAL :: charge_beta, charge_extended, &
294 REAL(kind=
dp) :: angle_alpha, angle_beta, angle_gamma
295 REAL(kind=
dp),
DIMENSION(3) :: abc
297 TYPE(cp_logger_type),
POINTER :: logger
298 TYPE(section_vals_type),
POINTER :: print_key
303 extension=
".subsysLog")
305 CALL timeset(routinen, handle)
310 i = count((/charge_occup, charge_beta, charge_extended/))
312 cpabort(
"Either only CHARGE_OCCUP, CHARGE_BETA, or CHARGE_EXTENDED can be selected")
319 IF (iw > 0)
WRITE (unit=iw, fmt=*)
" Writing out PDB file ", trim(record)
322 WRITE (unit=file_unit, fmt=
"(A6,T11,A)") &
327 WRITE (unit=file_unit, fmt=
"(A6,3F9.3,3F7.2)") &
328 "CRYST1", abc(1:3)*
angstrom, angle_alpha, angle_beta, angle_gamma
338 idres = atom_info%resid(i)
340 IF ((id1 /= atom_info%map_mol_num(i)) .OR. (id2 /= atom_info%map_mol_typ(i)))
THEN
342 id1 = atom_info%map_mol_num(i)
343 id2 = atom_info%map_mol_typ(i)
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)
362 WRITE (unit=line(55:60), fmt=
"(F6.2)") 0.0_dp
364 IF (
ASSOCIATED(atom_info%beta))
THEN
365 WRITE (unit=line(61:66), fmt=
"(F6.2)") atom_info%beta(i)
367 WRITE (unit=line(61:66), fmt=
"(F6.2)") 0.0_dp
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)
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)
387 WRITE (unit=file_unit, fmt=
"(A3)")
"END"
389 IF (iw > 0)
WRITE (unit=iw, fmt=*)
" Exiting "//routinen
392 "PRINT%TOPOLOGY_INFO/PDB_INFO")
394 CALL timestop(handle)
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.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
some minimal info about CP2K, including its version and license
character(len=default_string_length), public r_host_name
character(len= *), parameter, public compile_revision
character(len= *), parameter, public cp2k_version
character(len=default_string_length), public r_user_name
character(len=26), public r_datx
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.
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Utility routines for the memory handling.
Interface to the message passing library MPI.
Definition of physical constants:
real(kind=dp), parameter, public angstrom
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
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.