1! This file is part of xtb. 2! 3! Copyright (C) 2017-2020 Stefan Grimme 4! 5! xtb is free software: you can redistribute it and/or modify it under 6! the terms of the GNU Lesser General Public License as published by 7! the Free Software Foundation, either version 3 of the License, or 8! (at your option) any later version. 9! 10! xtb is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU Lesser General Public License for more details. 14! 15! You should have received a copy of the GNU Lesser General Public License 16! along with xtb. If not, see <https://www.gnu.org/licenses/>. 17 18module xtb_io_reader_turbomole 19 use xtb_mctc_accuracy, only : wp 20 use xtb_mctc_constants 21 use xtb_mctc_convert 22 use xtb_mctc_resize 23 use xtb_mctc_symbols, only : toNumber, symbolLength 24 use xtb_pbc_tools 25 use xtb_readin, getline => strip_line 26 use xtb_type_molecule 27 use xtb_type_reader 28 use xtb_type_vendordata, only : turbo_info 29 implicit none 30 private 31 32 public :: readMoleculeCoord 33 public :: readHessianTurbomole 34 35 36 logical, parameter :: debug = .false. 37 38 39contains 40 41 42subroutine readMoleculeCoord(mol, unit, status, iomsg) 43 type(TMolecule),intent(inout) :: mol 44 integer,intent(in) :: unit !< file handle 45 logical, intent(out) :: status 46 character(len=:), allocatable, intent(out) :: iomsg 47 character(len=:), allocatable :: line 48 character(len=:), allocatable :: cell_string, lattice_string 49 character(len=1), parameter :: flag = '$' 50 character(len=symbolLength), allocatable :: sym(:) 51 real(wp), allocatable :: coord(:, :), xyz(:, :) 52 real(wp) :: latvec(9), conv 53 integer :: error 54 integer :: iatom, i, j, natoms 55 integer :: periodic, cell_vectors 56 integer, parameter :: p_initial_size = 100 57 integer, parameter :: p_nlv(3) = [1, 4, 9] 58 integer, parameter :: p_ncp(3) = [1, 3, 6] 59 real(wp) :: cellpar(6), lattice(3, 3) 60 logical :: has_coord, has_periodic, has_lattice, has_cell 61 logical :: cartesian, coord_in_bohr, lattice_in_bohr 62 logical :: pbc(3) 63 64 status = .false. 65 66 allocate(sym(p_initial_size), source=' ') 67 allocate(coord(3, p_initial_size), source=0.0_wp) 68 69 iatom = 0 70 periodic = 0 71 has_coord = .false. 72 has_periodic = .false. 73 has_lattice = .false. 74 has_cell = .false. 75 cartesian = .true. 76 coord_in_bohr = .true. 77 lattice_in_bohr = .true. 78 lattice = 0.0_wp 79 pbc = .false. 80 81 error = 0 82 do while(error == 0) 83 call getline(unit, line, error) 84 if (index(line, flag) == 1) then 85 if (index(line, 'end') == 2) then 86 exit 87 88 elseif (.not.has_coord .and. index(line, 'coord') == 2) then 89 has_coord = .true. 90 ! $coord angs / $coord bohr / $coord frac 91 call select_unit(line, coord_in_bohr, cartesian) 92 coord_group: do while(error == 0) 93 call getline(unit, line, error) 94 if (index(line, flag) == 1) then 95 backspace(unit) 96 exit coord_group 97 endif 98 if (iatom >= size(coord, 2)) call resize(coord) 99 if (iatom >= size(sym)) call resize(sym) 100 iatom = iatom + 1 101 read(line, *, iostat=error) coord(:, iatom), sym(iatom) 102 enddo coord_group 103 104 elseif (.not.has_periodic .and. index(line, 'periodic') == 2) then 105 has_periodic = .true. 106 ! $periodic 0/1/2/3 107 read(line(10:), *, iostat=error) periodic 108 109 elseif (.not.has_lattice .and. index(line, 'lattice') == 2) then 110 has_lattice = .true. 111 ! $lattice bohr / $lattice angs 112 call select_unit(line, lattice_in_bohr) 113 cell_vectors = 0 114 lattice_string = '' 115 lattice_group: do while(error == 0) 116 call getline(unit, line, error) 117 if (index(line, flag) == 1) then 118 backspace(unit) 119 exit lattice_group 120 endif 121 cell_vectors = cell_vectors + 1 122 lattice_string = lattice_string // ' ' // line 123 enddo lattice_group 124 125 elseif (.not.has_cell .and. index(line, 'cell') == 2) then 126 has_cell = .true. 127 ! $cell bohr / $cell angs 128 call select_unit(line, lattice_in_bohr) 129 call getline(unit, cell_string, error) 130 if (debug) print*, cell_string 131 132 endif 133 endif 134 enddo 135 136 if (.not.has_coord .or. iatom == 0) then 137 iomsg = "coordinates not present, cannot work without coordinates" 138 return 139 endif 140 141 if (has_cell .and. has_lattice) then 142 iomsg = "both lattice and cell group are present" 143 return 144 endif 145 146 if (.not.has_periodic .and. (has_cell .or. has_lattice)) then 147 iomsg = "cell and lattice definition is present, but periodicity is not given" 148 return 149 endif 150 151 if (periodic > 0 .and. .not.(has_cell .or. has_lattice)) then 152 iomsg = "system is periodic but definition of lattice/cell is missing" 153 return 154 endif 155 156 if (.not.cartesian .and. periodic == 0) then 157 iomsg = "fractional coordinates do not work for molecular systems" 158 return 159 endif 160 161 natoms = iatom 162 allocate(xyz(3, natoms)) 163 164 if (any(toNumber(sym(:natoms)) == 0)) then 165 iomsg = "unknown element present" 166 return 167 endif 168 169 if (periodic > 0) pbc(:periodic) = .true. 170 171 if (has_cell) then 172 read(cell_string, *, iostat=error) latvec(:p_ncp(periodic)) 173 if (debug) print*, latvec(:p_ncp(periodic)) 174 if (lattice_in_bohr) then 175 conv = 1.0_wp 176 else 177 conv = aatoau 178 endif 179 select case(periodic) 180 case(1) 181 cellpar = [latvec(1)*conv, 1.0_wp, 1.0_wp, & 182 & pi/2, pi/2, pi/2] 183 case(2) 184 cellpar = [latvec(1)*conv, latvec(2)*conv, 1.0_wp, & 185 & pi/2, pi/2, latvec(3)*pi/180.0_wp] 186 case(3) 187 cellpar = [latvec(1:3)*conv, latvec(4:6)*pi/180.0_wp] 188 end select 189 call cell_to_dlat(cellpar, lattice) 190 endif 191 192 if (has_lattice) then 193 if (cell_vectors /= periodic) then 194 iomsg = "number of cell vectors does not match periodicity" 195 return 196 endif 197 read(lattice_string, *, iostat=error) latvec(:p_nlv(periodic)) 198 if (lattice_in_bohr) then 199 conv = 1.0_wp 200 else 201 conv = aatoau 202 endif 203 j = 0 204 do i = 1, periodic 205 lattice(:periodic,i) = latvec(j+1:j+periodic) * conv 206 j = j + periodic 207 enddo 208 endif 209 210 if (cartesian) then 211 if (coord_in_bohr) then 212 conv = 1.0_wp 213 else 214 conv = aatoau 215 endif 216 xyz(:, :) = coord(:, :natoms) * conv 217 else 218 call abc_to_xyz(natoms, lattice, coord, xyz) 219 endif 220 221 call init(mol, sym(:natoms), xyz, lattice=lattice, pbc=pbc) 222 ! save data on input format 223 mol%turbo = turbo_info(cartesian=cartesian, lattice=has_lattice, & 224 & angs_lattice=.not.lattice_in_bohr, & 225 & angs_coord=.not.coord_in_bohr) 226 227 status = .true. 228 229contains 230 231 subroutine select_unit(line, in_bohr, cartesian) 232 character(len=*), intent(in) :: line 233 logical, intent(out) :: in_bohr 234 logical, intent(out), optional :: cartesian 235 in_bohr = index(line, ' angs') == 0 236 if (present(cartesian)) cartesian = index(line, ' frac') == 0 237 end subroutine select_unit 238 239end subroutine readMoleculeCoord 240 241 242subroutine readHessianTurbomole(hessian, reader, mol, status, iomsg) 243 real(wp), intent(out) :: hessian(:, :) 244 type(TMolecule), intent(in) :: mol 245 type(TReader), intent(inout) :: reader 246 logical, intent(out) :: status 247 character(len=:), allocatable, intent(out) :: iomsg 248 character(len=:), allocatable :: line 249 character(len=32) :: buffer 250 integer :: error, ndim, ii, jj, jbatch, iline 251 252 status = .false. 253 iline = 1 254 255 call reader%read(line, error) 256 do while(error == 0) 257 if (index(line, '$hessian') == 1) exit 258 call reader%read(line, error) 259 iline = iline + 1 260 end do 261 262 if (error /= 0) then 263 iomsg = "Could not find $hessian data group" 264 return 265 end if 266 267 ndim = 3*len(mol) 268 269 rdlp: do ii = 1, ndim 270 do jj = 1, ndim, 5 271 call reader%read(line, error) 272 iline = iline + 1 273 if (error /= 0) exit rdlp 274 jbatch = min(jj+4, ndim) 275 read(line, '(5x, 5f15.10)', iostat=error) hessian(jj:jbatch, ii) 276 if (error /= 0) exit rdlp 277 end do 278 end do rdlp 279 280 if (error /= 0) then 281 if (is_iostat_end(error)) then 282 iomsg = "Unexpected end of file while reading hessian" 283 else 284 write(buffer, '(i0)') iline 285 iomsg = "Failed to read hessian in line "//trim(buffer) 286 end if 287 return 288 end if 289 290 status = .true. 291 292end subroutine readHessianTurbomole 293 294 295end module xtb_io_reader_turbomole 296