Read Input

c **********************************************************************
c     ReadInput.F
c **********************************************************************
c     ... adaptation of ReadLattice from Vernet
c **********************************************************************

      subroutine ReadInput

      use global_parameters
      use scalars
      use serial_arrays
      use pmeVar
      implicit none

      integer, dimension(48)   :: itext
      integer meshExp,sum
      integer itemp, isum, ir, ia, itrap
      character*11 fn
      real scaleFactor
      itrap = 0

c     hardwire filename
      
      fn   = 'lattice.dat'

c     ... initialize arrays

      a1     = zero
      a2     = zero
      a3     = zero
      b1     = zero
      b2     = zero
      b3     = zero
      iatnum = 0
      qatom  = zero
      atom   = zero
      numa   = 0
      itype  = 0

      if (whoami .eq. 0) then
         open(5,file=fn,status='old')
         read(5,'(48a1)') itext
         
c     get real box length (RboxL), box length is assumed to be 1 in the code
c     but this will allow numbers to be converted to actual length by
c     by multiplying by this number
         read(5,'(f16.10)') RboxL
         print *, "Box Length = ",RBoxL

c     Read rcut
         read(5,'(f16.10)') rcut
         print *, "cutoff radius = ",rcut

c     Read mesh density, convert to a power of two
          read(5,'(i2)') meshExp
          print *, 'mesh Exp =', meshExp
          meshSize = 2**meshExp
          print *, 'mesh size = ', meshSize

c     Read alpha
         read(5,'f15.10') alpha
         print *, "alpha = ",alpha

c     Read p, the order of the spline
         read(5,'i2') p
         print *, 'p ', p

c     ... read direct (a1, a2, and a3) and reciprocal (b1, b2, and b3) lattice
c         basis vectors, and total number of atoms (numa(0)).

c     get vectors defining suit direct and recpriprocal cells
         read(5,'(3f17.10)')
     &        (a1(ir),ir=1,3), (a2(ir),ir=1,3), (a3(ir),ir=1,3),
     &        (b1(ir),ir=1,3), (b2(ir),ir=1,3), (b3(ir),ir=1,3)
         scaleFactor = one/a1(1)

         read(5,'(i4)') numa(0)
         print *, "number of atoms = ",numa(0)
         if (numa(0) .gt. maxAtoms) then
           write(6,'(/,''==> Error:: numa(0) exceeds maxAtoms'')')
           write(6,'(4x,''numa(0) = '',i3)') numa(0)
           write(6,'(4x,''maxAtoms    = '',i3)') maxAtoms
           itrap = 1
         endif
      endif

c     ... each atomic position line in lattice.dat is of the following form
c         (example line from Al Sigma= 3 grain boundary file)
c
c  numt  at #   at. chg.      x                 y                 z
c   1    13      3.00     1.6662812455     8.6737885974     7.9316165296
c
c         numt is a label shared by atoms of the same type.  The value of this
c         label is stored in itype(atom #) and the total number of atoms of a
c         given type is stored in numa(<type label>).  The total number
c         of atoms overall (of any type) is stored in numa0).

      if (whoami .eq. 0) then
         qtot  = zero
         itemp = 1
         isum  = 0

         do ia = 1,numa(0)
           read(5,'(i4,i6,f10.2,3f17.10)')
     &          numt, iatnum(ia), qatom(ia), (atom(ir,ia),ir=1,3) 
           do ir =1,3
              atom(ir,ia) = atom(ir,ia)*scaleFactor
           end do
! ir (1,2,3) is the (x, y, z) coordinates of an atom, ia is the atom number

c     ... the next 'if' block assumes that all atoms of a given type
c         are grouped together in the input file.   First atoms of type
c         '1' are processed, then type '2', etc.  isum is re-initialized
c         when the next type is encountered to re-start number-of-atoms-of-
c         this-type counter.

           if (numt .gt. itemp) then
             itemp = numt
             isum = 0
           end if
           isum = isum + 1
           numa(itemp) = isum
           qtot = qtot + qatom(ia)
           itype(ia) = numt
         end do

         close(5)

         do ia = 1,numa(0)
           if (itype(ia) .gt. maxt) then
             write(6,'(/,''==> Error:: numt exceeds maxt'')')
             write(6,'(4x,''numt = '',i3)') itype(ia)
             write(6,'(4x,''maxt = '',i3)') maxt
             itrap = 1
           endif
         enddo
         sum = 0
         do ia = 1,numa(0)
            sum = sum + qatom(ia)
         end do
         if (sum .ne. 0) then
            print *, "error lattice has net charge"
            stop
         end if
      endif

c     ... given vectors a1, a2, a3, volume= a1 dotted into (a2 x a3)
      
      unitcv = abs( a1(1)*(a2(2)*a3(3) - a2(3)*a3(2))
     &            + a1(2)*(a2(3)*a3(1) - a2(1)*a3(3))
     &            + a1(3)*(a2(1)*a3(2) - a2(2)*a3(1)) )
      
      if(whoami.eq.0) then
        write(6,'(/,'' Lattice input:  '',48a1)') itext
        write(6,'(/,2x,'' primitive lattice vectors:'')')
        write(6,'(/,4x,'' a1:'',3f17.10)') (a1(ir),ir=1,3)
        write(6,'(4x,'' a2:'',3f17.10)') (a2(ir),ir=1,3)
        write(6,'(4x,'' a3:'',3f17.10)') (a3(ir),ir=1,3)
        write(6,'(/,4x,'' b1:'',3f17.10)') (b1(ir),ir=1,3)
        write(6,'(4x,'' b2:'',3f17.10)') (b2(ir),ir=1,3)
        write(6,'(4x,'' b3:'',3f17.10)') (b3(ir),ir=1,3)
        write(6,'(/,2x,'' unit cell volume:'',f17.10)') unitcv
        write(6,'(/,2x,'' atoms:'')')
        write(6,'(/,5x,''number'',5x,''atomic number'',23x,''position'')
     &')
        do ia = 1,numa(0), 1 
          write(6,'(/,6x,i3,11x,i3,5x,3f17.10)')
     &          ia, iatnum(ia), (atom(ir,ia),ir=1,3)
        end do
      endif

      if(itrap .eq. 1) then
         print *, 'error'
         stop
      end if
c     put in error handler to terminate program

      return
      end