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_vasp
19   use xtb_mctc_accuracy, only : wp
20   use xtb_mctc_convert
21   use xtb_mctc_strings
22   use xtb_mctc_systools
23   use xtb_type_molecule
24   use xtb_pbc_tools
25   use xtb_mctc_symbols, only : toNumber, symbolLength
26   use xtb_type_molecule
27   use xtb_type_vendordata, only : vasp_info
28   implicit none
29   private
30
31   public :: readMoleculeVasp
32
33
34   logical, parameter :: debug = .false.
35
36
37contains
38
39
40subroutine readMoleculeVasp(mol, unit, status, iomsg)
41   type(TMolecule),intent(out) :: mol
42   integer,intent(in) :: unit
43   logical, intent(out) :: status
44   character(len=:), allocatable, intent(out) :: iomsg
45   real(wp) :: lattice(3,3)
46   integer, allocatable :: ncount(:)
47   logical              :: selective=.false. ! Selective dynamics
48   logical              :: cartesian=.true.  ! Cartesian or direct
49
50   integer :: natoms
51   real(wp) :: ddum,latvec(3)
52   real(wp) xx(10),scalar
53   real(wp) :: coord(3)
54   character(len=:),allocatable :: line
55   character(len=2*symbolLength) :: args(256),args2(256)
56   real(wp), allocatable :: xyz(:, :)
57   character(len=symbolLength), allocatable :: sym(:)
58
59   integer i,j,k,nn,ntype,ntype2,atnum,i_dummy1,i_dummy2,ncheck
60
61   integer :: iat, idum, err
62
63   status = .false.
64
65   lattice=0
66
67   ncheck=0
68   ntype=0
69   ! first line contains the symbols of different atom types
70   call getline(unit,line,err)
71   if (err.ne.0) then
72      iomsg = "Could not read POSCAR"
73      return
74   endif
75   if (debug) print'(">",a)',line
76   call parse(line,' ',args,ntype)
77
78   ! this line contains the global scaling factor,
79   call getline(unit,line,err)
80   if (err.ne.0) then
81      iomsg = "Could not read POSCAR"
82      return
83   endif
84   if (debug) print'(">",a)',line
85   read(line,*,iostat=err) ddum
86   if (err.ne.0) then
87      iomsg = "Could not read POSCAR"
88      return
89   endif
90   ! the Ang->au conversion is included in the scaling factor
91   if (debug) print'("->",g0)',ddum
92   scalar = ddum*aatoau
93
94   ! reading the lattice constants
95   do i=1,3
96      call getline(unit,line,err)
97      if (err.ne.0) then
98         iomsg = "Could not read lattice from POSCAR"
99         return
100      endif
101      if (debug) print'("->",a)',line
102      read(line,*,iostat=err) latvec
103      if (err.ne.0) then
104         iomsg = "Could not read lattice from POSCAR"
105         return
106      endif
107      lattice(:,i) = latvec * scalar
108   enddo
109   ! Either here are the numbers of each element,
110   ! or (>vasp.5.1) here are the element symbols
111   call getline(unit,line,err)
112   if (err.ne.0) then
113      iomsg = "Could not read POSCAR"
114      return
115   endif
116   if (debug) print'(">",a)',line
117   ! try to read first element
118   read(line,*,iostat=err) idum
119   ! CONTCAR files have additional Element line here since vasp.5.1
120   if (err.ne.0) then
121      call parse(line,' ',args,ntype)
122      call getline(unit,line,err)
123      if (debug) print'("->",a)',line
124      if (err.ne.0) then
125         iomsg = "Could not read POSCAR"
126         return
127      endif
128   endif
129   call parse(line,' ',args2,nn)
130   if (nn.ne.ntype) then
131      iomsg = 'Error reading number of atomtypes'
132      return
133   endif
134
135   allocate(ncount(nn), source = 0)
136   ncheck = 0
137   do i = 1, nn
138      read(args2(i), *, iostat=err) ncount(i)
139      iat = toNumber(args(i))
140      if (iat < 1 .or. ncount(i) < 1) then
141         iomsg = 'unknown element.'
142         return
143      endif
144   enddo
145
146   natoms = sum(ncount)
147   allocate(sym(natoms))
148   allocate(xyz(3, natoms))
149
150   k = 0
151   do i = 1, nn
152      do j = 1, ncount(i)
153         k = k+1
154         sym(k) = trim(args(i))
155      enddo
156   enddo
157
158   call getline(unit, line, err)
159   if (err.ne.0) then
160      iomsg = "Could not read POSCAR"
161      return
162   endif
163   if (debug) print'(">",a)', line
164   line=adjustl(line)
165   if (line(:1).eq.'s' .or. line(:1).eq.'S') then
166      selective = .true.
167      call getline(unit,line,err)
168      if (debug) print'("->",a)', line
169      if (err.ne.0) then
170         iomsg = "Could not read POSCAR"
171         return
172      endif
173      line = adjustl(line)
174   endif
175
176   cartesian=(line(:1).eq.'c' .or. line(:1).eq.'C' .or. &
177      &       line(:1).eq.'k' .or. line(:1).eq.'K')
178   do i = 1, natoms
179      call getline(unit, line, err)
180      if (err.ne.0) then
181         iomsg = "Could not read geometry from POSCAR"
182         return
183      endif
184      if (debug) print'("-->",a)',line
185      read(line, *, iostat=err) coord
186      if (err.ne.0) then
187         iomsg = "Could not read geometry from POSCAR"
188         return
189      endif
190
191      if (cartesian) then
192         xyz(:,i) = coord*scalar
193      else
194         xyz(:,i) = matmul(lattice, coord)
195      endif
196
197   enddo
198
199   call init(mol, sym, xyz, lattice=lattice)
200   ! save information about this POSCAR for later
201   mol%vasp = vasp_info(scale=ddum, selective=selective, cartesian=cartesian)
202
203   status = .true.
204
205end subroutine readMoleculeVasp
206
207
208end module xtb_io_reader_vasp
209