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_genformat 19 use xtb_mctc_accuracy, only : wp 20 use xtb_mctc_convert 21 use xtb_mctc_strings 22 use xtb_mctc_symbols, only : toNumber, symbolLength 23 use xtb_mctc_systools 24 use xtb_pbc_tools 25 use xtb_type_molecule 26 use xtb_type_reader 27 use xtb_type_vendordata, only : vasp_info 28 implicit none 29 private 30 31 public :: readMoleculeGenFormat 32 public :: readHessianDFTBPlus 33 34 35contains 36 37 38subroutine readMoleculeGenFormat(mol, unit, status, iomsg) 39 logical, parameter :: debug = .false. 40 type(TMolecule),intent(out) :: mol 41 integer,intent(in) :: unit 42 logical, intent(out) :: status 43 character(len=:), allocatable, intent(out) :: iomsg 44 45 character(len=:), allocatable :: line 46 integer :: natoms, nspecies, iatom, dummy, isp, ilat, error 47 logical :: cartesian, periodic 48 real(wp) :: coord(3), lattice(3, 3) 49 character(len=1) :: variant 50 character(len=symbolLength), allocatable :: species(:), sym(:) 51 real(wp), allocatable :: xyz(:, :), abc(:, :) 52 53 status = .false. 54 55 call next_line(unit, line, error) 56 read(line, *, iostat=error) natoms, variant 57 if (error /= 0 .or. natoms < 1) then 58 iomsg = 'could not read number of atoms' 59 return 60 endif 61 62 allocate(species(natoms)) 63 allocate(sym(natoms)) 64 allocate(xyz(3, natoms)) 65 allocate(abc(3, natoms)) 66 67 select case(variant) 68 case('c', 'C') 69 cartesian = .true. 70 periodic = .false. 71 case('s', 'S') 72 cartesian = .true. 73 periodic = .true. 74 case('f', 'F') 75 cartesian = .false. 76 periodic = .true. 77 case default 78 iomsg = 'invalid input version' 79 return 80 endselect 81 82 call next_line(unit, line, error) 83 call parse(line, ' ', species, nspecies) 84 if (any(toNumber(species(:nspecies)) == 0)) then 85 iomsg = 'unknown atom type present' 86 return 87 endif 88 89 do iatom = 1, natoms 90 call next_line(unit, line, error) 91 read(line, *, iostat=error) dummy, isp, coord 92 if (error /= 0) then 93 iomsg = 'could not read coordinates from file' 94 return 95 endif 96 sym(iatom) = species(isp) 97 if (cartesian) then 98 xyz(:, iatom) = coord * aatoau 99 else 100 abc(:, iatom) = coord 101 endif 102 enddo 103 104 if (periodic) then 105 call next_line(unit, line, error) 106 if (error /= 0) then 107 iomsg = 'missing lattice information' 108 return 109 endif 110 do ilat = 1, 3 111 call next_line(unit, line, error) 112 read(line, *, iostat=error) coord 113 if (error /= 0) then 114 iomsg = 'could not read lattice from file' 115 return 116 endif 117 lattice(:, ilat) = coord * aatoau 118 enddo 119 if (.not.cartesian) then 120 call abc_to_xyz(natoms, lattice, abc, xyz) 121 endif 122 call init(mol, sym, xyz, lattice=lattice) 123 else 124 call init(mol, sym, xyz) 125 endif 126 127 mol%vasp = vasp_info(cartesian=cartesian) 128 129 status = .true. 130 131contains 132 133subroutine next_line(unit, line, error) 134 integer,intent(in) :: unit 135 character(len=:), allocatable, intent(out) :: line 136 integer, intent(out) :: error 137 integer :: ihash 138 139 error = 0 140 do while(error == 0) 141 call getline(unit, line, error) 142 ihash = index(line, '#') 143 if (ihash > 0) line = line(:ihash-1) 144 if (len_trim(line) > 0) exit 145 enddo 146 line = trim(adjustl(line)) 147end subroutine next_line 148 149end subroutine readMoleculeGenFormat 150 151 152subroutine readHessianDFTBPlus(hessian, reader, mol, status, iomsg) 153 real(wp), intent(out) :: hessian(:, :) 154 type(TMolecule), intent(in) :: mol 155 type(TReader), intent(inout) :: reader 156 logical, intent(out) :: status 157 character(len=:), allocatable, intent(out) :: iomsg 158 character(len=:), allocatable :: line 159 character(len=32) :: buffer 160 integer :: error, ndim, ii, jj, jbatch, iline 161 162 status = .false. 163 iline = 0 164 165 if (error /= 0) then 166 iomsg = "Could not find $hessian data group" 167 return 168 end if 169 170 ndim = 3*len(mol) 171 172 read(reader%unit, *, iostat=error) ((hessian(jj, ii), jj=1, ndim), ii=1, ndim) 173 174 if (error /= 0) then 175 if (is_iostat_end(error)) then 176 iomsg = "Unexpected end of file while reading hessian" 177 else 178 write(buffer, '(i0)') iline 179 iomsg = "Failed to read hessian in line "//trim(buffer) 180 end if 181 return 182 end if 183 184 status = .true. 185 186end subroutine readHessianDFTBPlus 187 188 189end module xtb_io_reader_genformat 190