Next: Self Interaction Term Up: Appendix A: Code Previous: Main Program   Contents

```c **********************************************************************
c **********************************************************************
c **********************************************************************

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')

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
print *, "Box Length = ",RBoxL

print *, "cutoff radius = ",rcut

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

print *, "alpha = ",alpha

c     Read p, the order of the spline
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
&        (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)

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)
&          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
```

Thomas G Dimiduk 2004-04-15