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