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 : MODULE topology_xyz
9 : USE cp_log_handling, ONLY: cp_get_default_logger,&
10 : cp_logger_type
11 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
12 : cp_print_key_unit_nr
13 : USE cp_parser_methods, ONLY: parser_get_next_line,&
14 : parser_get_object,&
15 : read_integer_object
16 : USE cp_parser_types, ONLY: cp_parser_type,&
17 : parser_create,&
18 : parser_release
19 : USE cp_units, ONLY: cp_unit_to_cp2k
20 : USE input_section_types, ONLY: section_vals_type
21 : USE kinds, ONLY: default_string_length,&
22 : dp,&
23 : max_line_length
24 : USE memory_utilities, ONLY: reallocate
25 : USE message_passing, ONLY: mp_para_env_type
26 : USE periodic_table, ONLY: nelem,&
27 : ptable
28 : USE string_table, ONLY: id2str,&
29 : s2s,&
30 : str2id
31 : USE topology_types, ONLY: atom_info_type,&
32 : topology_parameters_type
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_xyz'
38 :
39 : PRIVATE
40 : PUBLIC :: read_coordinate_xyz
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief ...
46 : !> \param topology ...
47 : !> \param para_env ...
48 : !> \param subsys_section ...
49 : !> \author Teodoro Laino
50 : ! **************************************************************************************************
51 2022 : SUBROUTINE read_coordinate_xyz(topology, para_env, subsys_section)
52 : TYPE(topology_parameters_type) :: topology
53 : TYPE(mp_para_env_type), POINTER :: para_env
54 : TYPE(section_vals_type), POINTER :: subsys_section
55 :
56 : CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_xyz'
57 :
58 : CHARACTER(LEN=default_string_length) :: my_default_index, strtmp
59 : CHARACTER(LEN=max_line_length) :: error_message
60 : INTEGER :: frame, handle, ian, iw, j, natom
61 : LOGICAL :: my_end
62 : TYPE(atom_info_type), POINTER :: atom_info
63 : TYPE(cp_logger_type), POINTER :: logger
64 : TYPE(cp_parser_type) :: parser
65 :
66 1011 : CALL timeset(routineN, handle)
67 :
68 1011 : NULLIFY (logger)
69 1011 : logger => cp_get_default_logger()
70 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XYZ_INFO", &
71 1011 : extension=".subsysLog")
72 :
73 1011 : atom_info => topology%atom_info
74 :
75 1011 : IF (iw > 0) THEN
76 : WRITE (UNIT=iw, FMT="(T2,A)") &
77 3 : "BEGIN of XYZ data read from file "//TRIM(topology%coord_file_name)
78 : END IF
79 :
80 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env, &
81 1011 : parse_white_lines=.TRUE.)
82 :
83 : ! Element is assigned on the basis of the atm_name
84 1011 : topology%aa_element = .TRUE.
85 :
86 : natom = 0
87 1011 : frame = 0
88 1011 : CALL parser_get_next_line(parser, 1)
89 1047 : Frames: DO
90 : ! Atom numbers
91 1047 : CALL parser_get_object(parser, natom)
92 1047 : frame = frame + 1
93 1047 : IF (frame == 1) THEN
94 1011 : CALL reallocate(atom_info%id_molname, 1, natom)
95 1011 : CALL reallocate(atom_info%id_resname, 1, natom)
96 1011 : CALL reallocate(atom_info%resid, 1, natom)
97 1011 : CALL reallocate(atom_info%id_atmname, 1, natom)
98 1011 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
99 1011 : CALL reallocate(atom_info%atm_mass, 1, natom)
100 1011 : CALL reallocate(atom_info%atm_charge, 1, natom)
101 1011 : CALL reallocate(atom_info%occup, 1, natom)
102 1011 : CALL reallocate(atom_info%beta, 1, natom)
103 1011 : CALL reallocate(atom_info%id_element, 1, natom)
104 36 : ELSE IF (natom > SIZE(atom_info%id_atmname)) THEN
105 0 : CPABORT("Atom number differs in different frames!")
106 : END IF
107 : ! Dummy line
108 1047 : CALL parser_get_next_line(parser, 2)
109 34256 : DO j = 1, natom
110 : ! Atom coordinates
111 34220 : READ (parser%input_line, *) strtmp, &
112 34220 : atom_info%r(1, j), &
113 34220 : atom_info%r(2, j), &
114 68440 : atom_info%r(3, j)
115 34220 : error_message = ""
116 34220 : CALL read_integer_object(strtmp, ian, error_message)
117 34220 : IF (LEN_TRIM(error_message) == 0) THEN
118 : ! Integer value found: assume atomic number, check it, and load
119 : ! the corresponding element symbol if valid
120 32 : IF ((ian < 0) .OR. (ian > nelem)) THEN
121 : error_message = "Invalid atomic number <"//TRIM(strtmp)// &
122 0 : "> found in the xyz file <"//TRIM(topology%coord_file_name)//">!"
123 0 : CPABORT(TRIM(error_message))
124 : ELSE
125 32 : atom_info%id_atmname(j) = str2id(s2s(ptable(ian)%symbol))
126 : END IF
127 : ELSE
128 34188 : atom_info%id_atmname(j) = str2id(s2s(strtmp))
129 : END IF
130 : ! For default, set atom name to residue name to molecule name
131 34220 : WRITE (my_default_index, '(I0)') j
132 34220 : atom_info%id_molname(j) = str2id(s2s(TRIM(id2str(atom_info%id_atmname(j)))//TRIM(my_default_index)))
133 34220 : atom_info%id_resname(j) = atom_info%id_molname(j)
134 34220 : atom_info%resid(j) = 1
135 34220 : atom_info%id_element(j) = atom_info%id_atmname(j)
136 34220 : atom_info%atm_mass(j) = HUGE(0.0_dp)
137 34220 : atom_info%atm_charge(j) = -HUGE(0.0_dp)
138 34220 : IF (iw > 0) THEN
139 : WRITE (UNIT=iw, FMT="(T2,A4,3F8.3,2X,A)") &
140 12 : TRIM(id2str(atom_info%id_atmname(j))), &
141 12 : atom_info%r(1, j), &
142 12 : atom_info%r(2, j), &
143 12 : atom_info%r(3, j), &
144 24 : ADJUSTL(TRIM(id2str(atom_info%id_molname(j))))
145 : END IF
146 34220 : atom_info%r(1, j) = cp_unit_to_cp2k(atom_info%r(1, j), "angstrom")
147 34220 : atom_info%r(2, j) = cp_unit_to_cp2k(atom_info%r(2, j), "angstrom")
148 34220 : atom_info%r(3, j) = cp_unit_to_cp2k(atom_info%r(3, j), "angstrom")
149 : ! If there's a white line or end of file exit.. otherwise read other available
150 : ! snapshots
151 34220 : CALL parser_get_next_line(parser, 1, at_end=my_end)
152 34220 : my_end = my_end .OR. (LEN_TRIM(parser%input_line) == 0)
153 34256 : IF (my_end) THEN
154 1011 : IF (j /= natom) &
155 : CALL cp_abort(__LOCATION__, &
156 : "Number of lines in XYZ format not equal to the number of atoms."// &
157 : " Error in XYZ format. Very probably the line with title is missing or is empty."// &
158 0 : " Please check the XYZ file and rerun your job!")
159 : EXIT Frames
160 : END IF
161 : END DO
162 : END DO Frames
163 1011 : CALL parser_release(parser)
164 :
165 1011 : IF (iw > 0) THEN
166 : WRITE (UNIT=iw, FMT="(T2,A)") &
167 3 : "END of XYZ frame data read from file "//TRIM(topology%coord_file_name)
168 : END IF
169 :
170 1011 : topology%natoms = natom
171 1011 : topology%molname_generated = .TRUE.
172 :
173 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
174 1011 : "PRINT%TOPOLOGY_INFO/XYZ_INFO")
175 :
176 1011 : CALL timestop(handle)
177 :
178 3033 : END SUBROUTINE read_coordinate_xyz
179 :
180 : END MODULE topology_xyz
|
