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