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