78 CHARACTER(len=*),
PARAMETER :: routinen =
'read_coordinate_xtl'
79 INTEGER,
PARAMETER :: nblock = 1000
80 REAL(kind=
dp),
PARAMETER :: threshold = 1.0e-6_dp
82 CHARACTER(LEN=default_string_length) :: strtmp
83 INTEGER :: dimensions, handle, icol, ii, isym, iw, &
84 jj, natom, natom_orig, newsize
85 INTEGER,
DIMENSION(3) :: periodic
86 LOGICAL :: check, found, my_end
87 REAL(kind=
dp) :: pfactor, threshold2
88 REAL(kind=
dp),
DIMENSION(3) :: cell_angles, cell_lengths, r, r1, r2, s, &
90 REAL(kind=
dp),
DIMENSION(3, 3) :: rot_mat
99 extension=
".subsysLog")
100 CALL timeset(routinen, handle)
102 pfactor =
section_get_rval(subsys_section,
"TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
107 CALL reallocate(atom_info%id_molname, 1, nblock)
108 CALL reallocate(atom_info%id_resname, 1, nblock)
110 CALL reallocate(atom_info%id_atmname, 1, nblock)
112 CALL reallocate(atom_info%atm_mass, 1, nblock)
113 CALL reallocate(atom_info%atm_charge, 1, nblock)
116 CALL reallocate(atom_info%id_element, 1, nblock)
118 IF (iw > 0)
WRITE (iw, *)
" Reading in XTL file ", trim(
topology%coord_file_name)
123 begin_line=.false., search_from_begin_of_file=.true.)
125 IF (iw > 0)
WRITE (iw,
'(/,A)')
" XTL_INFO| TITLE :: "//trim(parser%input_line(parser%icol:))
130 begin_line=.false., search_from_begin_of_file=.true.)
132 IF (iw > 0)
WRITE (iw,
'(A)')
" XTL_INFO| DIMENSION :: "//trim(parser%input_line(parser%icol:))
134 IF (dimensions /= 3)
THEN
135 cpabort(
"XTL file with working DIMENSION different from 3 cannot be parsed!")
146 begin_line=.false., search_from_begin_of_file=.true.)
148 cpabort(
"The field CELL was not found in XTL file! ")
173 CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
180 begin_line=.false., search_from_begin_of_file=.true.)
182 cpabort(
"The field ATOMS was not found in XTL file! ")
185 found = (index(parser%input_line,
"NAME X Y Z") /= 0)
187 cpabort(
"The field ATOMS in XTL file, is not followed by name and coordinates tags! ")
191 DO WHILE (index(parser%input_line,
"EOF") == 0)
194 IF (natom >
SIZE(atom_info%id_molname))
THEN
195 newsize = int(pfactor*natom)
196 CALL reallocate(atom_info%id_molname, 1, newsize)
197 CALL reallocate(atom_info%id_resname, 1, newsize)
199 CALL reallocate(atom_info%id_atmname, 1, newsize)
200 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
201 CALL reallocate(atom_info%atm_mass, 1, newsize)
202 CALL reallocate(atom_info%atm_charge, 1, newsize)
205 CALL reallocate(atom_info%id_element, 1, newsize)
209 atom_info%id_atmname(natom) =
str2id(strtmp)
211 atom_info%id_resname(natom) = atom_info%id_molname(natom)
212 atom_info%resid(natom) = 1
213 atom_info%id_element(natom) = atom_info%id_atmname(natom)
220 s = atom_info%r(1:3, natom)
226 threshold2 = threshold*threshold
230 r1 = atom_info%r(1:3, ii)
231 DO jj = ii + 1, natom
232 r2 = atom_info%r(1:3, jj)
233 r =
pbc(r1 - r2, cell)
235 check = (dot_product(r, r) >= threshold2)
242 begin_line=.false., search_from_begin_of_file=.true.)
244 IF (iw > 0)
WRITE (iw,
'(A)')
" XTL_INFO| Symmetry Infos :: "//trim(parser%input_line(parser%icol:))
249 begin_line=.false., search_from_begin_of_file=.true.)
251 cpwarn(
"The field SYM MAT was not found in XTL file! ")
252 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" XTL_INFO| Number of atoms before applying symmetry operations :: ", natom
253 IF (iw > 0)
WRITE (iw,
'(A10,1X,3F12.6)') (trim(
id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
260 icol = index(parser%input_line,
"SYM MAT") + 8
261 READ (parser%input_line(icol:), *) ((rot_mat(ii, jj), jj=1, 3), ii=1, 3), transl_vec(1:3)
262 loop_over_unique_atoms:
DO ii = 1, natom_orig
264 r1 = matmul(rot_mat, atom_info%r(1:3, ii)) + transl_vec
268 r2 = atom_info%r(1:3, jj)
269 r =
pbc(r1 - r2, cell)
271 IF (dot_product(r, r) <= threshold2)
THEN
280 IF (natom >
SIZE(atom_info%id_molname))
THEN
281 newsize = int(pfactor*natom)
282 CALL reallocate(atom_info%id_molname, 1, newsize)
283 CALL reallocate(atom_info%id_resname, 1, newsize)
285 CALL reallocate(atom_info%id_atmname, 1, newsize)
286 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
287 CALL reallocate(atom_info%atm_mass, 1, newsize)
288 CALL reallocate(atom_info%atm_charge, 1, newsize)
291 CALL reallocate(atom_info%id_element, 1, newsize)
293 atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
294 atom_info%id_molname(natom) = atom_info%id_molname(ii)
295 atom_info%id_resname(natom) = atom_info%id_resname(ii)
296 atom_info%resid(natom) = atom_info%resid(ii)
297 atom_info%id_element(natom) = atom_info%id_element(ii)
298 atom_info%r(1:3, natom) = r1
300 END DO loop_over_unique_atoms
302 begin_line=.false., search_from_begin_of_file=.false.)
305 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" XTL_INFO| Number of symmetry operations :: ", isym
306 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" XTL_INFO| Number of total atoms :: ", natom
307 IF (iw > 0)
WRITE (iw,
'(A10,1X,3F12.6)') (trim(
id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
314 CALL reallocate(atom_info%id_molname, 1, natom)
315 CALL reallocate(atom_info%id_resname, 1, natom)
317 CALL reallocate(atom_info%id_atmname, 1, natom)
320 CALL reallocate(atom_info%atm_charge, 1, natom)
323 CALL reallocate(atom_info%id_element, 1, natom)
328 "PRINT%TOPOLOGY_INFO/XTL_INFO")
329 CALL timestop(handle)