Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : MODULE topology_cp2k
9 :
10 : USE cell_types, ONLY: cell_transform_input_cartesian,&
11 : cell_type,&
12 : scaled_to_real
13 : USE cp_log_handling, ONLY: cp_get_default_logger,&
14 : cp_logger_type
15 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
16 : cp_print_key_unit_nr
17 : USE cp_parser_methods, ONLY: parser_get_object,&
18 : parser_test_next_token,&
19 : read_integer_object
20 : USE cp_parser_types, ONLY: cp_parser_type,&
21 : parser_create,&
22 : parser_release
23 : USE cp_units, ONLY: cp_unit_to_cp2k
24 : USE input_section_types, ONLY: section_get_ival,&
25 : section_get_rval,&
26 : section_vals_get,&
27 : section_vals_get_subs_vals,&
28 : section_vals_type,&
29 : section_vals_val_get,&
30 : section_vals_val_set
31 : USE kinds, ONLY: default_string_length,&
32 : dp,&
33 : max_line_lengTh
34 : USE memory_utilities, ONLY: reallocate
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE periodic_table, ONLY: nelem,&
37 : ptable
38 : USE string_table, ONLY: id2str,&
39 : s2s,&
40 : str2id
41 : USE topology_types, ONLY: atom_info_type,&
42 : topology_parameters_type
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_cp2k'
48 :
49 : PRIVATE
50 : PUBLIC :: read_coordinate_cp2k
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Read the CP2K &COORD section from an external file, i.e. read
56 : !> atomic coordinates and molecule/residue information in CP2K format.
57 : !> \param topology ...
58 : !> \param para_env ...
59 : !> \param subsys_section ...
60 : !> \date 17.01.2011 (Creation, MK)
61 : !> \author Matthias Krack (MK)
62 : !> \version 1.0
63 : ! **************************************************************************************************
64 36 : SUBROUTINE read_coordinate_cp2k(topology, para_env, subsys_section)
65 :
66 : TYPE(topology_parameters_type) :: topology
67 : TYPE(mp_para_env_type), POINTER :: para_env
68 : TYPE(section_vals_type), POINTER :: subsys_section
69 :
70 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_coordinate_cp2k'
71 :
72 : CHARACTER(LEN=default_string_length) :: string
73 : CHARACTER(LEN=max_line_length) :: error_message
74 : INTEGER :: handle, i, ian, iw, natom, newsize, &
75 : number_of_atoms
76 : LOGICAL :: eof, explicit, scaled_coordinates
77 : REAL(KIND=dp) :: pfactor, unit_conv
78 : REAL(KIND=dp), DIMENSION(3) :: r
79 : TYPE(atom_info_type), POINTER :: atom_info
80 : TYPE(cell_type), POINTER :: cell
81 : TYPE(cp_logger_type), POINTER :: logger
82 : TYPE(cp_parser_type), POINTER :: parser
83 : TYPE(section_vals_type), POINTER :: coord_section
84 :
85 12 : CALL timeset(routineN, handle)
86 :
87 12 : NULLIFY (coord_section)
88 12 : NULLIFY (logger)
89 12 : NULLIFY (parser)
90 12 : logger => cp_get_default_logger()
91 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XYZ_INFO", &
92 12 : extension=".subsysLog")
93 :
94 : ! Check if there is a &COORD section
95 12 : coord_section => section_vals_get_subs_vals(subsys_section, "COORD")
96 12 : CALL section_vals_get(coord_section, explicit=explicit)
97 12 : IF (explicit) THEN
98 10 : CALL section_vals_val_get(coord_section, "UNIT", c_val=string)
99 10 : CALL section_vals_val_get(coord_section, "SCALED", l_val=scaled_coordinates)
100 : ELSE
101 : ! The default is Cartesian coordinates in Angstrom
102 2 : scaled_coordinates = .FALSE.
103 2 : string = "angstrom"
104 : END IF
105 12 : unit_conv = cp_unit_to_cp2k(1.0_dp, TRIM(string))
106 :
107 12 : atom_info => topology%atom_info
108 12 : cell => topology%cell_muc
109 :
110 12 : IF (iw > 0) THEN
111 : WRITE (UNIT=iw, FMT="(T2,A)") &
112 6 : "BEGIN of COORD section data read from file "//TRIM(topology%coord_file_name)
113 : END IF
114 :
115 12 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
116 12 : number_of_atoms = section_get_ival(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS")
117 12 : IF (number_of_atoms < 1) THEN
118 12 : newsize = 1000
119 : ELSE
120 0 : newsize = number_of_atoms
121 : END IF
122 :
123 12 : CALL reallocate(atom_info%id_molname, 1, newsize)
124 12 : CALL reallocate(atom_info%id_resname, 1, newsize)
125 12 : CALL reallocate(atom_info%resid, 1, newsize)
126 12 : CALL reallocate(atom_info%id_atmname, 1, newsize)
127 12 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
128 12 : CALL reallocate(atom_info%atm_mass, 1, newsize)
129 12 : CALL reallocate(atom_info%atm_charge, 1, newsize)
130 12 : CALL reallocate(atom_info%occup, 1, newsize)
131 12 : CALL reallocate(atom_info%beta, 1, newsize)
132 12 : CALL reallocate(atom_info%id_element, 1, newsize)
133 :
134 12 : topology%molname_generated = .FALSE.
135 : ! Element is assigned on the basis of the atm_name
136 12 : topology%aa_element = .TRUE.
137 :
138 36 : ALLOCATE (parser)
139 12 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
140 :
141 12 : natom = 0
142 : DO
143 1164 : CALL parser_get_object(parser, object=string, newline=.TRUE., at_end=eof)
144 1164 : IF (eof) EXIT
145 1152 : natom = natom + 1
146 1152 : IF (natom > SIZE(atom_info%id_atmname)) THEN
147 0 : newsize = INT(pfactor*natom)
148 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
149 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
150 0 : CALL reallocate(atom_info%resid, 1, newsize)
151 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
152 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
153 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
154 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
155 0 : CALL reallocate(atom_info%occup, 1, newsize)
156 0 : CALL reallocate(atom_info%beta, 1, newsize)
157 0 : CALL reallocate(atom_info%id_element, 1, newsize)
158 : END IF
159 1152 : error_message = ""
160 1152 : CALL read_integer_object(string, ian, error_message)
161 1152 : IF (LEN_TRIM(error_message) == 0) THEN
162 : ! Integer value found: assume atomic number, check it, and load
163 : ! the corresponding element symbol if valid
164 0 : IF ((ian < 0) .OR. (ian > nelem)) THEN
165 : error_message = "Invalid atomic number <"//TRIM(string)// &
166 0 : "> found in the xyz file <"//TRIM(topology%coord_file_name)//">!"
167 0 : CPABORT(TRIM(error_message))
168 : ELSE
169 0 : atom_info%id_atmname(natom) = str2id(s2s(ptable(ian)%symbol))
170 : END IF
171 : ELSE
172 1152 : atom_info%id_atmname(natom) = str2id(s2s(string))
173 : END IF
174 : ! Read x, y, and z coordinate of the current atom
175 4608 : DO i = 1, 3
176 4608 : CALL parser_get_object(parser, object=r(i))
177 : END DO
178 1152 : IF (scaled_coordinates) THEN
179 576 : CALL scaled_to_real(atom_info%r(1:3, natom), r, cell)
180 : ELSE
181 2304 : atom_info%r(1:3, natom) = r(1:3)*unit_conv
182 576 : CALL cell_transform_input_cartesian(cell, atom_info%r(1:3, natom))
183 : END IF
184 1152 : IF (parser_test_next_token(parser) /= "EOL") THEN
185 768 : CALL parser_get_object(parser, object=string)
186 768 : atom_info%id_molname(natom) = str2id(s2s(string))
187 768 : IF (parser_test_next_token(parser) /= "EOL") THEN
188 384 : CALL parser_get_object(parser, object=string)
189 384 : atom_info%id_resname(natom) = str2id(s2s(string))
190 : ELSE
191 1152 : atom_info%id_resname(natom) = atom_info%id_molname(natom)
192 : END IF
193 : ELSE
194 384 : string = ""
195 384 : WRITE (UNIT=string, FMT="(I0)") natom
196 384 : atom_info%id_molname(natom) = str2id(s2s(TRIM(id2str(atom_info%id_atmname(natom)))//TRIM(string)))
197 384 : atom_info%id_resname(natom) = atom_info%id_molname(natom)
198 1536 : topology%molname_generated = .TRUE.
199 : END IF
200 1152 : atom_info%resid(natom) = 1
201 1152 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
202 1152 : atom_info%atm_mass(natom) = HUGE(0.0_dp)
203 1152 : atom_info%atm_charge(natom) = -HUGE(0.0_dp)
204 1152 : IF (iw > 0) THEN
205 : WRITE (UNIT=iw, FMT="(T2,A,3F20.8,2(2X,A))") &
206 2304 : TRIM(id2str(atom_info%id_atmname(natom))), atom_info%r(1:3, natom), &
207 576 : ADJUSTL(TRIM(id2str(atom_info%id_molname(natom)))), &
208 1152 : ADJUSTL(TRIM(id2str(atom_info%id_resname(natom))))
209 : END IF
210 1164 : IF (natom == number_of_atoms) EXIT
211 : END DO
212 :
213 12 : CALL parser_release(parser)
214 12 : DEALLOCATE (parser)
215 :
216 12 : topology%natoms = natom
217 12 : CALL reallocate(atom_info%id_molname, 1, natom)
218 12 : CALL reallocate(atom_info%id_resname, 1, natom)
219 12 : CALL reallocate(atom_info%resid, 1, natom)
220 12 : CALL reallocate(atom_info%id_atmname, 1, natom)
221 12 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
222 12 : CALL reallocate(atom_info%atm_mass, 1, natom)
223 12 : CALL reallocate(atom_info%atm_charge, 1, natom)
224 12 : CALL reallocate(atom_info%occup, 1, natom)
225 12 : CALL reallocate(atom_info%beta, 1, natom)
226 12 : CALL reallocate(atom_info%id_element, 1, natom)
227 :
228 12 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", i_val=natom)
229 :
230 12 : IF (iw > 0) THEN
231 : WRITE (UNIT=iw, FMT="(T2,A)") &
232 6 : "END of COORD section data read from file "//TRIM(topology%coord_file_name)
233 : END IF
234 :
235 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
236 12 : "PRINT%TOPOLOGY_INFO/XYZ_INFO")
237 :
238 12 : CALL timestop(handle)
239 :
240 12 : END SUBROUTINE read_coordinate_cp2k
241 :
242 : END MODULE topology_cp2k
|
