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_xyz
19   use xtb_mctc_accuracy, only : wp
20   use xtb_mctc_convert
21   use xtb_mctc_symbols, only : toNumber, symbolLength
22   use xtb_type_molecule
23   use xtb_pbc_tools
24   use xtb_readin, only : getline => strip_line
25   implicit none
26   private
27
28   public :: readMoleculeXYZ
29
30
31   logical, parameter :: debug = .false.
32
33
34contains
35
36
37subroutine readMoleculeXYZ(mol, unit, status, iomsg)
38   class(TMolecule), intent(out) :: mol
39   integer, intent(in) :: unit
40   logical, intent(out) :: status
41   character(len=:), allocatable, intent(out) :: iomsg
42   integer  :: ii, n, iat
43   character(len=symbolLength), allocatable :: sym(:)
44   real(wp), allocatable :: xyz(:, :)
45   real(wp) :: x, y, z
46   real(wp) :: conv
47
48   character(len=:),allocatable :: message
49   character(len=:),allocatable :: line
50   character(len=symbolLength) :: chdum
51   integer  :: err
52
53   status = .false.
54
55   conv = aatoau
56
57   read(unit,*,iostat=err) n
58   if (err.ne.0) then
59      iomsg = "Could not read number of atoms, check format!"
60      return
61   endif
62
63   if (n.lt.1) then
64      iomsg = "Found no atoms, cannot work without atoms!"
65      return
66   endif
67
68   allocate(sym(n))
69   allocate(xyz(3, n))
70
71   ! drop next record
72   read(unit,'(a)')
73
74   ii = 0
75   do while (ii < n)
76      call getline(unit,line,err)
77      if (is_iostat_end(err)) exit
78      if (err.ne.0) then
79         iomsg = "Could not read geometry from Xmol file"
80         return
81      endif
82      if (debug) print'(">",a)',line
83      read(line,*,iostat=err) chdum, x, y, z
84      if (err.ne.0) then
85         iomsg = "Could not parse coordinates from Xmol file"
86         return
87      endif
88      if (debug) print'("->",a)',chdum
89      if (debug) print'("->",3g0)',x, y, z
90
91      iat = toNumber(chdum)
92      if (debug) print'("->",g0)',iat
93      if (iat > 0) then
94         ii = ii+1
95         sym(ii) = trim(chdum)
96         xyz(:,ii) = [x, y, z]*conv
97      else
98         iomsg = "Unknown element symbol: '"//trim(chdum)//"'"
99         return
100      endif
101   enddo
102
103   if (ii.ne.n) then
104      iomsg = "Atom number missmatch in Xmol file"
105      return
106   endif
107
108   call init(mol, sym, xyz)
109
110   status = .true.
111
112end subroutine readMoleculeXYZ
113
114
115end module xtb_io_reader_xyz
116