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