program atoms
      implicit real(a-h,o-z)
      implicit integer(i-n)
c       include 'version.h'
c-*-fortran-*-
      character*9 vrsion
      parameter (vrsion='2.50 ')
c----------------------------------------------------------------------
c        copyright 1998-2001  Bruce Ravel
c        copyright 1993-1997  University of Washington
c             written by      Bruce Ravel
c                  e-mail     ravel@phys.washington.edu
c                     WWW     http://feff.phys.washington.edu/~ravel
c           please use email for communication with the author
c
c    This program is free software; you can redistribute it and/or modify
c    it under the terms of the GNU General Public License as published by
c    the Free Software Foundation; either version 2, or (at your option)
c    any later version.

c    This program is distributed in the hope that it will be useful,
c    but WITHOUT ANY WARRANTY; without even the implied warranty of
c    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c    GNU General Public License for more details.

c    You should have received a copy of the GNU General Public License
c    along with GNU Emacs; see the file COPYING.  If not, write to
c    the Free Software Foundation, 675 Massachusettes Ave,
c    Cambridge, MA 02139, USA.

c    Everyone is granted permission to copy, modify and redistribute this
c    and related files provided:
c      1. All copies contain this copyright notice.
c      2. All modified copies shall carry a prominant notice stating who
c         made modifications and the date of such modifications.
c      3. The name of the modified file be changed.
c      4. No charge is made for this software or works derived from it.
c         This clause shall not be construed as constraining other software
c         distributed on the same medium as this software, nor is a
c         distribution fee considered a charge.
c----------------------------------------------------------------------
c      brief description of the code:
c
c  atoms writes a list of atomic coordinates for any crystal given its
c  crystallographic information. the list will be sorted by radial
c  distance from an atom chosen as the central atom. atoms also estimates
c  the bulk absorption and density of the material and various corrections
c  to xafs data due to experimental effects.
c----------------------------------------------------------------------
c  comments blocks that follow:
c       glossary of variables
c       descriptions of runtime error codes
c       sample input file
c       version history
c
c  the code begins after the first occurance of this string:  %%%%
c  the first executable statement begins after its second occurance
c----------------------------------------------------------------------
c >>> glossary of variables
c       include 'glossary'
c=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  this is a list of globaly or commonly used parameters and
c  variables with brief descriptions of their purposes
c  () = dimension, [] = set of values, {} = alternate name
c
c=-=-=-=-=-=-=-=-=-=-=-=-=-= parameters =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  eps:    a small number for use in floating point logicals
c  etok:   conversion between ev and invang
c  iat:    maximum number of unique atoms
c  natx:   maximum number of atoms in cluster
c  ndopx:  maximum number of dopants at a site
c  nfit:   maximum number of points for linear regression
c  ngeomx: maximum number of one bounce flags
c  ntitx:  maximum number of title lines
c  nvals:  maximum number of parameters for linear regression
c  nwdx:   maximum number of words for bwords
c  nlogx:  maximum number of logic flags in array logic
c  nexafs: maximum number of items calculated from mcmaster data
c  vrsion: version number
c
c=-=-=-=-=-=-=-=-=-=-=-=-=-= integers =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  iatom:  number of atoms in atom list
c  ibasis: number of atoms in basis list
c  ipt:    (iat), number of positions of a unique atom in unit cell
c  iptful: (iat), number of positions of a unique atom in overfull cell
c  imult:  (iat), multiplicity at each site
c  ispa:   [1..230], number of space group in itxc
c  isyst:  [1..8], identifies crystal system
c  itot:   number of atoms found in cluster
c  job:    counts complete runs through code with multiple inputs
c  ngeom:  (ngeomx) one-bounce flags for geom.dat
c  ns:     {also ng}, simple multiplicity at a point
c  nsites: equal iatom or, if using basis, equals ibasis
c  ntit:   (9), number of title lines
c
c=-=-=-=-=-=-=-=-=-=-=-=-=-= logicals =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  vaxflg: true for compilation on vms, for open status of files
c  stdout: true for stdin/stdout operation, false for disk operation
c  expnd:  true to use expanded atoms list, false for normal
c  logic:  (nlogx) array of various run control and diagnostic flags,
c          see arrays.map
c
c=-=-=-=-=-=-=-=-=-=-=-=-=  characters  =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  core:   atomic symbol of central atom
c  cnum:   used to write out numbers with messag
c  dopant: (iat,ndopx) atomic symbols of atoms at site (main+dopants)
c  edge:   k or l3 edge of core atom
c  messg:  string sent to messag
c  outfil: output file name
c  spcgrp: space group of crystal
C  title:  (ntitx), user input title lines
c
c=-=-=-=-=-=-=-=-=-=-=-=-=-=-= reals =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c  ampslf: self-absorption amplitude correction
c  atlis:  (natx,8), list of atoms in cluster, fractional, and cartesian
c  cell:   (6), a,b,c,alpha,beta, and gamma of cell
c  dmax:   radius of cluster
c  fulcel: (iat,192,3), coordinates of atoms in overfull cell
c  fs:     (3,3,24),
c  gasses: (3) percent by pressure in i0 chamber of argon, krypton, nitrogen
c  percnt: (iat,ndopx) percent substitution of dopants for replaced atoms
c  exafs:  (10) array of values from mcamster calculation
c  st:     (iat,192,3), fractioanl coords of points after bravais transl.
c  trmtx:  (3,3), transformation matrix between cartesian and cell-axis
c  ts:     (3,24), fractional positions of points before bravais transl.
c  x:      (iat), fractional coordinate of unique atoms
c  y:      (iat), fractional coordinate of unique atoms
c  z:      (iat), fractional coordinate of unique atoms
c
c=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c >>> runtime error codes
c       include 'runtime'
c----------------------------------------------------------------------
c			     Error Codes
c
c  All error messages are of 1 or more lines.  They begin with three
c  asterixes (***) followed by a token and end with a period
c  followed by a dash (.-).  If an input file keyword appears in the
c  message, it should be surrounded by double quotes.
c
c  When an error is caught, it should set an error code.  The codes are
c    0   no error (and no message)
c    1   informational (tokens: Caution  Position  Ending)
c    2   warning (token: Warning)
c    3   error (token: Error)
c
c  Informational messages are those that do not indicate a problem,
c    but do indicate something the user should know about.  Currently
c    used informational messages are for debugging messages,
c    positional messages, impending early termination of the program,
c    and non-serious cautions.
c  Warnings are for things that should not impede the progress of the
c    program, but which may produce output that is not what the user
c    intended.  Warnings include unknown keywords and exceeded parameters.
c  Errors are things the preclude the continuation of the program.
c----------------------------------------------------------------------
c=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c >>> sample input file:
c       include 'sample'
c
c           title pbtio3 n&k,10k,a=3.885,c=4.139
c           space  p 4 m m
c           a=3.885   c=4.139
c           rmax=10 core=ti
c           atom
c           ! at.type  x        y       z      tag
c              pb     0.0      0.0     0.0
c              ti     0.5      0.5     0.5377
c              o      0.5      0.5     0.1118   o1
c              o      0.0      0.5     0.6174   o2
c           -------------------------------------------------
c=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c >>> history
c       include 'history'
c    1.0    got it working
c    2.0    satisfied copyright criteria in cluster expansion routines
c    2.1    added indexing (8/92)
c    2.2    added absorption length (9/92)
c    2.3    improved handling of spc grp notation, fixed bcc bug,
c           fixed floating point logic (1/93)
c    2.4    handle bases, atom list option, shift (3/93), mcmaster
c           correction (4/93)
c  approaching 3 too quickly, switch to hundredths
c    2.41   mcmaster correction corrected, i0 and self absorption
c           corrections, new cluster expansion routine, overfull
c           unit cell, simple dopant, more error checking on input,
c           multiple title lines, untab, shift warning and
c           version number in feff.inp (6/93)
c    2.42   modularized in anticipation of dafs applications (6/93)
c           general housekeeping, tags, improved dopants, core keyword
c           required, case insensitivity in code, geom.dat + subshell
c           sorting, self-absorption bug fixed (8/93 and 10/93)
c    2.42.a lower case and true robustness, revamped input parsing for
c           better error checking, handle overlapping atoms better, hex
c           C, V symmetry groups, corrected tri/hex bug in atinpt, no
c           mcmaster for Z<13, diagnostic messages from input file (1/94)
c           corrected bug in dopants,
c           more housekeeping (2/94) calculate delta mu from core rather
c           than total mu and reinstate mcmaster for Z<13 (3/94)
c           silly bug in atinpt that trashed nonorthoganal groups fixed,
c           also improved error messages (3/94)
c    2.42.b fixed mishandling of rhombohedral groups (6/94)
c    2.42.c pass parameters as arguments for easier hacking, had to pass
c           some memory space to certain subroutines, fixed a bug
c           in subroutine ref that was affecting geout, 72 character
c           names for output files, pass io number for feff.inp,
c           l1&l2, move origin into readin (8/94)
c    2.42.d contains full fuctionality of 2.43.c except for dafs apps (1/95)
c    2.42.e contains full fuctionality of 2.43.d except for dafs apps (2/95)
c    2.42.f contains full fuctionality of 2.43.f except for dafs apps
c     & up   parallel development to 2.43 without dafs apps
c
c    2.43   tabulate chantler's data for interpolation, wrote fcal --
c           an engine for using chantler's data, wrote a0.f for reading
c           chantler and cromer-mann and making a0, write a0 as a five
c           column file suitable for feffit, no energy corrections (8/94)
c    2.43.b new a0 filename, allow neg. q's, fix null atoms handling,
c           fixed mishandling of trigonal (not rhomb.) groups (9/94)
c           added warning message for low symmetry cells, added my name
c           to feff.inp and run time messages(11/94), fix some error
c           messages (1/95)
c    2.43.c full sasaki tables, fixed espilon bug in multip, added p1
c           option  (1/95)
c    2.43.d after UWXAFS 3.0 release -- fixed bug with dmax < min distance
c           (2/95)
c    skip subversion e
c    2.43.f fixed bug in self-absorption calculation, added krypton to I0,
c           geom.dat and feff.dat write same number of atoms (4/95)
c    2.43.g no file name conflicts, much improved internal documentation
c           (7/95)
c
c    2.44. & 2.45.  handles different settings of low symmetry crystals,
c                   improve handling of schoenflies notation, better error
c                   messages for mcmaster calculations (8/95), compile time
c                   switch for stdin/stdout (9/95), typo in syschk, in
c                   multip minor bug re. imult for atom with coord=1
c                   edge energy in feff.inp, allow L1&l2 (10/95)
c    2.44 & up parallels 2.45 without dafs apps
c    2.45.a fixed bug in f.p. comparisons in indexing, five digits in
c           feff.inp, geom.dat no longer depends on scratch file (11/95)
c    2.45.b central atom gets tag, reflect.dat (12/95)
c
c    2.46 & up parallels 2.47 without dafs apps and extended atoms list
c    2.47   tighten modularity & combine variables into arrays, fix problem
c           with nofx and polyft by explicitly carrying around misc.f file,
c           information for fluorescence normalization,  move atoms list
c           parsing to external subroutine to allow expanded atoms list,
c           rmax card to 5 digits + 0.00001, fixed potential numerical
c           problem in metric (2/96) made groups module (3/96) angles in
c           monoclinic settings,
c    2.4?.b rework makefile using uninclude perl script and unzipped
c           modules, added xanes card
c    2.4?.c added keywords "feff8" and "correction", added optional
c           feff8 output (BR Jan 16 1998)
c
c    to do: new crystal, fix monoclinic angles
c
c    2.50   change structure of source tree, use libraries and include
c           files in a more sensible manner, (BR Mar  2 1998) (goals
c           for 2.50pre1: modularize crystl+groups+clustr, remove logic()
c
c    changes in functionality from versions before 2.50
c    1. no longer support more than one job per file
c    2. consistent message syntax for post-processing purposes
c
c    goals for 2.50pre2: clean up remaining code, allow spcgrp=ispa,
c                        end keyord wade through code cleaning up
c                        run-time messages
c    goals for 2.50pre3: begin work on monoclinic
c    other goals for 2.50: revise document, write crystl user guide,
c       more output types, see wishlist
c----------------------------------------------------------------------
c%%%%
c       include 'atparm.h'
c-*-fortran-*-
c  These parameters are the variable size declarations for the program
      parameter (iat=50, natx=800, ntitx=9, ndopx=4, ngeomx=natx)
      parameter (neptx=2**11, maxln=natx)
      parameter (nlogx=28, nexafs=13, ndbgx=10)
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:
c
c  iat:    maximum number of unique atom positions
c  natx:   maximum size of atomic cluster
c  ntitx:  maximum number of title lines
c  ndopx:  maximum number of dopants at any site
c  ngeomx: maximum number of lines written to geom.dat
c  neptx:  maximum number of energy points in dafs output files
c  maxln:  maximum number of lines written to feff.inp
c  nlogx:  number of logical parameters in logic array
c  nexafs: number of mcmaster paramters in exafs array
c  ndbgx:  maximum size of debugging code numbers = 2**ndbgx
c------------------------------------------------------------------------

c  stdout=.true. for reading from standard input and writing to standard
c  output
      logical stdout
      parameter (stdout=.false.)

c  expnd=.true. for reading expanded atoms list, false for normal atoms list
      logical expnd
      parameter (expnd=.false.)

c  vaxflg=.true. for compilation on a VMS machine, used for opening
c  files in such a way that VMS version numbers are used
      logical vaxflg
      parameter (vaxflg=.false.)

c  crystl.h contains all information relevant to a unit cell
c       include 'crystl.h'
c-*-fortran-*-

c  various parameters used by module crystl
      common /cryint/ iabs, iatom, ibasis, isystm, ispa, iperm, nsites,
     $            ipt(iat), idop(iat), imult(iat)
      save /cryint/

      parameter(nsysm=4, nshwrn=4)
      character*2  dopant(iat,ndopx)
      character*10 spcgrp, tag(iat)
      character*74 shwarn(nshwrn)
      character*77 sysmes(nsysm)
      common /crystr/ shwarn, sysmes, dopant, tag, spcgrp
      save /crystr/

      logical syserr, shift
      common /crylog/ syserr, shift
      save /crylog/

      dimension trmtx(3,3), st(iat,192,3)
      dimension cell(6), x(iat), y(iat), z(iat)
      dimension percnt(iat,ndopx)
      common /cryflt/ trmtx, st, cell, x, y, z, percnt
      save /cryflt/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:  (* = user input, % = error handling, ! = output needed to
c              construct cluster, the rest are used internally)
c
c * iabs:   index of absorber in unique coordinate list
c * iatom:  >=1 if atoms list is used, else =1
c * ibasis: =1 if basis list is used, else =0
c   isystm: index of crystal system (1..7)=(mono,orth,<not used>,tetr,
c           cubic,hex,triclinic)
c   ispa:   space group index, 1-230 from IXTC
c             0:       not recognized, error in input symbol
c             1-2:     triclinic
c             3-15:    monoclinic
c             16-74:   orthorhombic
c             75-142:  tetragonal
c             143-167: trigonal
c             168-194: hexagonal
c             195-230: cubic
c   iperm:  permutation index for non-standard settings, used in crystl
c             1:     default value -- no permutation necessary
c             1-6:   6 orthorhombic settings (abc, cab, bca, a-cb, ba-c, -cba)
c             11-12: 2 monoclinic settings, z-axis unique, y-axis unique
c             21-22: 2 tetragonal settings, standard and rotated
c ! ipt:    (iat) number of positions of unique atom in unit cell (1..192)
c   imult:  (iat) workspace for calculating multiplicities
c
c % sysmes: 4 line message if space group and axes/angles don't match
c % syserr: true if space group and axes/angles don't match
c % shwarn: 3 line message if space group may require a shift vector
c % shift:  true if space group may require a shift vector
c
c * dopant*2:  (iat,ndopx) matrix with all host and dopant atomic symbols
c * percnt:    (iat,ndopx) matrix with occupancies of hosts and dopants
c * tag*10:    (iat) character tag for each unique site in cell
c * spcgrp*10: space group symbol.  On output it is the short
c              Hermann-Maguin symbol in standard setting.  On input
c              spcgrp can be any short HM, Schoenflies, a number
c              between 1 and 230, or one of a small set of special
c              words (fcc, bcc, etc.).  Other symbol conventions
c              (full HM symbol, Shubnikov, 1935 ITXC, etc.) are not
c              and never will be used.
c
c * x,y,z:  (iat) arrays of fractional coordinates of unique positions
c                 in unit cell
c * cell:   (6) array of a,b,c,alpha,beta,gamma
c
c ! trmtx:  (3,3) transformation matrix between cell-axis and cartesian
c                 bases, see subroutine trans in clustr
c ! st:     (iat,192,3) fractional coordinates of all atoms in unit cell,
c                       first arg refers to unique atom list, second
c                       to position in cell, third is xyz.
c
c           fyi: 192 is the largest possible number of equivalent
c                positions in a cell of any symmetry. see, for example,
c                cubic f m 3 c.
c----------------------------------------------------------------------
c  exafs.h contains information relevant to exafs calculations
c       include 'exafs.h'
c-*-fortran-*-
      common /exaflt/ gasses(3), exafs(nexafs), rmax
      save /exaflt/

      common /exaint/ iedge, iexerr
      save /exaint/

      character*10 core, edge*2
      common /exastr/ core, edge
      save /exastr/

      logical lfluo
      common /exalog/ lfluo
      save /exalog/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:  (* = user input, ! = output, % = error handling)
c
c  * gasses:  (3) percent pressure of argon, krypton, nitrogen in i0 chamber
c         gasses(1)       percentage of argon
c         gasses(2)       percentage of krypton
c         gasses(3)       percentage of nitrogen
c  ! exafs:   (13) amu,delmu,spgrav,sigmm,qrtmm,ampslf,sigslf,qrtslf,
c                  sigi0,qrti0,muf,mub,mue
c         exafs(1):       total mu
c         exafs(2):       delta mu
c         exafs(3):       speciffic gravity
c         exafs(4):       mcmaster sigma^2
c         exafs(5):       mcmaster C4
c         exafs(6):       self absorption amplitude correction
c         exafs(7):       self absorption sigma^2
c         exafs(8):       self absorption C4
c         exafs(9):       i0 corrcetion sigma^2
c         exafs(10):      i0 correction C4
c  * iedge:   edge for calulation, 1=K 2=L1 3=L2 4=L3
c  * iexerr:  exit error code, 0=no prob, 1=info, 2=warning, 3=error
c  * core*10: tag of absorbing atom
c    lfluo:   true if fluorescence corrections were calculated
c----------------------------------------------------------------------
c  unit.h contains information over-full unit cell information
c       include 'unit.h'
c-*-fortran-*-

      common /uninum/ iptful(iat), fullcl(iat,192,3)
      save /uninum/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:
c
c  * iptful: (iat) number of positions of unique atom in overfull unit cell
c  * fullcl: (iat,192,3) fractional coords of all atoms in overfull unit cell
c----------------------------------------------------------------------

c  see glossary for descriptions of variables, see arrays.map for
c  description of logic, gasses, exafs
      parameter(m=8)
      character*2  elemnt(iat), test, vrsn*9
      character*10 inpgrp, noantg(iat)
      character*72 title(ntitx), outfil, afname, refile
      character*78 messg
      character*78 module, string
      logical      logic(nlogx)
      complex      anot(neptx)
      dimension    ngeom(ngeomx)
      dimension    atlis(natx,8), qvect(3)
      dimension    nrefl(3)


c  reserve some array space for use in modules cluster and output
      dimension    cltmp1(natx),  cltmp2(natx),  cltmp3(natx,m)
      logical      cltmp4(natx)
      character*10 optmp1(maxln)
      dimension    optmp2(maxln), optmp3(maxln), optmp4(maxln),
     $             optmp5(maxln), ioptp6(maxln), ioptp7(maxln)

c------------------------------------------------------------------------
c-- dafs stuff
c       dimension    f0(iat,ndopx), usqr(iat)
c       complex      fcore(neptx)
c c  from block data sasaki
c       parameter(nelem=92, ndatx=233)
c       common /fdat/ nfdata(nelem), engrd(nelem,ndatx),
c      $            fp(nelem,ndatx), fpp(nelem,ndatx)
c c  this is the crommer-mann common block
c       common /sk/ aa(1926)
c------------------------------------------------------------------------

c------------------------------------------------------------------------
c  formats needed at top level: version # (4000), line of '=' (4010),
c                 ... aN,(55-N)x ...
 4000 format('ATOMS ',a9,46x,'by Bruce Ravel')
 4010 format(75('='))

c%%%%
c---------------------------------------------------------------------
c  run time messages: lines of '=' and version number
c  set unit number of feff.inp
      if (.not.stdout) then
          ifeff = 2
          write(messg,4010)
          call messag(messg)
          write(messg,4000)vrsion
          call messag(messg)
          write(messg,4010)
          call messag(messg)
      else
c  this clause may trigger a compiler warning regarding there being no
c  possible path to this spot in the code.  this is an unavoidable
c  result of using F77 and no preprocessor and should pose no problem
c  during execution of the code.
          ifeff = 6
      endif

c---------------------------------------------------------------------
c  initialize some top-level variables
      vrsn  = vrsion
      test  = 'ab'
      logic(1) = .false.
      module = 'Atoms'

c 10    continue

c=====================================================================
c  atoms module 1:  initialize, read input file, error checking
c=====================================================================
c      print*,'beginning module 1'
      call readin(iat, natx, ntitx, ndopx, ngeomx, neptx, nlogx,
     $            ifeff, ntit, iatom, ibasis, iabs, isystm,
     $            iperm, ipt, iptful, idop, ngeom, iedge,
     $            nepts, nsites, nrefl, nnoan,
     $            vrsn, spcgrp, inpgrp, title, outfil, afname,
     $            refile, tag, noantg, edge, core, elemnt, dopant,
     $            x, y, z, cell, rmax, st, fullcl, atlis,
     $            percnt, gasses, qvect, egr, anot,
     $            logic, stdout, vaxflg, expnd)

      if (logic(1)) goto 99

c=====================================================================
c  atoms module 2:  decode space group, determine unit cell contents
c=====================================================================
      string = 'the crystl module'
      if (logic(20)) call positn(module, string)
      iunidb = 0
      if (logic(21)) icrydb = icrydb + 2**0
      call crystl(icrydb, iercry)

      if ( (.not.stdout).and.syserr ) then
          call messag(' ')
          call messag(' *** Warning ')
          do 20 i=1,nsysm
            call messag(sysmes(i))
 20       continue
          call messag('The calculation will be finished, but you '//
     $                'might want to edit your')
          call messag('crystallographic input data and try again.-')
          call messag(' ')
      endif

c=====================================================================
c  atoms module 3:  perform various calculations using mcmaster tables
c=====================================================================
      if (logic(28)) then
          string = 'the mcm module'
          if (logic(7)) then
              if (logic(20)) call positn(module, string)
              iexadb = 0
              if (logic(22)) iexadb = iexadb + 2**0
              if (logic(14)) iexadb = iexadb + 2**1
              if (logic(15)) iexadb = iexadb + 2**2
              if (logic(16)) iexadb = iexadb + 2**3
              call mcm(iexadb)
          endif
          if (.not.lfluo) logic(3)=.false.
      endif

c=====================================================================
c  atoms module 4:  construct unit.dat
c=====================================================================
      string = 'the unit module'
      if (logic(20)) call positn(module, string)
      iunidb = 0
      if (logic(23)) iunidb = iunidb + 2**0
      if (logic(8))  iunidb = iunidb + 2**1
      if (logic(9))  iunidb = iunidb + 2**2
      call unit(iunidb, ntit, title, vaxflg, ieruni)

c=====================================================================
c  atoms module 5:  construct structure factor for dafs applications
c=====================================================================
c       if (logic(10).or.logic(11)) then
c           string = 'the ascat module'
c           if (logic(20)) call positn(module, string)
c           call ascat(iat, ntitx, ndopx, neptx, nlogx,
c      $             nsites, ipt, idop, ntit, nepts, nrefl, nnoan,
c      $             edge, core, dopant, title, afname, refile, vrsn,
c      $             tag, noantg,
c      $             st, usqr, cell, percnt, qvect, egr, fcore, f0,
c      $             anot, logic, vaxflg)
c       endif

c=====================================================================
c  atoms module 6:  expand cluster around central atom
c=====================================================================
      if (logic(7)) then
          string = 'the clustr module'
          if (logic(20)) call positn(module, string)
          call clustr(iat,natx,ngeomx,nlogx,
     $            iabs,nsites,iperm,ipt,itot,ngeom,
     $            cell,trmtx,rmax,st,atlis,
     $            cltmp1,cltmp2,cltmp3,cltmp4,
     $            logic)

c=====================================================================
c  atoms module 7:  write feff.inp, geom.dat
c=====================================================================
          string = 'the output module'
          if (logic(20)) call positn(module, string)
          call output(iat, natx, ntitx, ndopx, ngeomx, maxln, nlogx,
     $            nexafs,
     $            ifeff, iabs, itot, ntit, idop, ngeom, imult,
     $            title, tag, edge, core, dopant, outfil, elemnt,
     $            vrsn, percnt, exafs, atlis,
     $            logic, vaxflg,
     $            optmp1,optmp2,optmp3,optmp4,optmp5,ioptp6,ioptp7)

c#####################################################################

c---------------------------------------------------------------------
c  run time message: output file and line of '='
          if ((logic(6)).and.(outfil.ne.'list')) then
              if (logic(13)) then
                  call messag('geom.dat overwritten.')
              else
                  call messag('geom.dat written.')
              endif
          endif
          if (.not.stdout) then
              ii = istrln(outfil)
              if (logic(12)) then
                  call messag(outfil(:ii)//' overwritten.')
              else
                  call messag('Output written to '//outfil(:ii))
              endif
          endif
      else
          call messag('not writing feff.inp.')
      endif
      if (.not.stdout) then
          write(messg,4010)
          call messag(messg)
      endif
c---------------------------------------------------------------------

c      if (.not.logic(1)) goto 10
99    continue
      stop
c  end main program atom
      end
      subroutine readin(iat,natx,ntitx,ndopx,ngeomx,neptx,nlogx,ifeff,
     $      ntit,iatom,ibasis,iabs,isystm,iperm,
     $      ipt,iptful,idop,ngeom,iedge,nepts,nsites,nrefl,nnoan,
     $      vrsion,spcgrp,inpgrp,title,outfil,afname,refile,
     $      tag,noantg,edge,core,elemnt,dopant,
     $      x,y,z,cell,dmax,st,fullcl,atlis,
     $      percnt,gasses,qvect,egr,anot,
     $      logic,stdout,vaxflg,expnd)
c=====================================================================
c  atoms module 1:  initialize, read input file, error checking
c=====================================================================
c  this module consists of the following subroutines and functions:
c     readin atchck atinit atinpt atspec dopfix getatm groups origin
c     rh2hex schfix settng spcchk systm
c=====================================================================
      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

c      parameter (iat=50, natx=800, ntitx=9, ndopx=4, ngeomx=800)
c      parameter (neptx=2**11)

      character*2  elemnt(iat)
      character*2  edge, dopant(iat,ndopx), vrsion*9, test
      character*10 spcgrp, inpgrp, tag(iat), core, geodat, noantg(iat)
      character*72 title(ntitx),outfil,afname,outf,refile
      character*74 messg
      logical      logic(nlogx), stdout, vaxflg, expnd, shift
      complex      anot(neptx)
      dimension    ipt(iat), iptful(iat), idop(iat),
     $             ngeom(ngeomx), nrefl(3)
      dimension    x(iat), y(iat), z(iat)
      dimension    cell(6), qvect(3), gasses(3)
      dimension    st(iat,192,3), fullcl(iat,192,3),  atlis(natx,8)
      dimension    percnt(iat,ndopx)

      parameter(nshwrn=4)
      character*74 shwarn(nshwrn)
      logical      shft

 4000 format(35('*-'),'*')
 4010 format(1x,'* ',a75)
 4020 format(1x,' ')
 4100 format(1x,'* This feff.inp file generated by ATOMS, version ',
     $       a9,/,' * ATOMS written by and copyright (c) Bruce Ravel',
     $       ', 1992-1999',/)

c------------------------------------------------------------
c  initialize everything and read the input file
c  then check the consistency of the input values and
c  determine system from cell constants
c      call messag('  initializing...')
      call atinit(iat,natx,ntitx,ndopx,ngeomx,neptx,nlogx,
     $            title,elemnt,tag,noantg,spcgrp,edge,core,
     $            dopant,outfil,afname,geodat,refile,
     $            iatom,ibasis,dmax,ispa,iperm,idop,ngeom,iedge,iabs,
     $            ipt,iptful,isyst,nepts,nrefl,nnoan,
     $            x,y,z,cell,st,fullcl,atlis,
     $            percnt,gasses,qvect,egr,anot,
     $            logic,stdout)

c      call messag('  reading atom.inp...')
      call atinpt(iat,ntitx,ndopx,nlogx,
     $            title,ntit,iatom,ibasis,iabs,dmax,ispa,iperm,
     $            idop,iedge,nepts,nrefl,isystm,nnoan,
     $            elemnt,tag,noantg,edge,core,spcgrp,inpgrp,
     $            outfil,shwarn,afname,refile,dopant,
     $            x,y,z,cell,percnt,gasses,qvect,egr,
     $            logic, stdout, expnd, shift)

      if (logic(1)) goto 99

c  this is a lame work-around for not yet having headers in readin
      call igtisp(ispa)
c      call messag('  consistancy checks...')
      call atchck(iat,ndopx,core,dopant,edge,
     $            iatom,ibasis,ispa,idop,
     $            x,y,z,cell,dmax,gasses,qvect)

c a few more chores before leaving
      nsites = iatom
      if (logic(2)) nsites=ibasis

      inquire(file=outfil,exist=logic(12))
      if (logic(6)) inquire(file=geodat,exist=logic(13))

      test = 'ab'
      outf=outfil
      call case(test,outf)
      if (stdout) then
          write(ifeff,4100)vrsion
c          call origin(spcgrp, warn, wrning)
      else
          if (outf.eq.'list') then
              inquire(file='atoms.lis',exist=logic(12))
          else
              if (.not.vaxflg) then
                  open(unit=ifeff,file=outfil,status='unknown')
              else
                  open(unit=ifeff,file=outfil,status='new')
              endif
              write(ifeff,4100)vrsion
          endif
      endif

c             --- give warning about space groups that need shift
c      call origin(spcgrp, warn, wrning)
      call gtshft(shft,shwarn)
      if (shft) then
          if (ifeff.ne.6) then
              write(messg,4000)
c              call messag(messg)
              call messag(' ')
              call messag(' *** Warning:')
          endif
          do 20 i=1,nshwrn
            if (ifeff.eq.6) then
                call messag('* '//shwarn(i))
            else
                call messag(shwarn(i))
                write(ifeff,4010)shwarn(i)
            endif
 20       continue
          write(ifeff,4020)
          if (ifeff.ne.6) then
              write(messg,4000)
c              call messag(messg)
              call messag(' ')
          endif
          ierr = 2
      endif

99    continue

      return
c  end of module readin
      end

      subroutine igtisp(isp)
c       include 'atparm.h'
c-*-fortran-*-
c  These parameters are the variable size declarations for the program
      parameter (iat=50, natx=800, ntitx=9, ndopx=4, ngeomx=natx)
      parameter (neptx=2**11, maxln=natx)
      parameter (nlogx=28, nexafs=13, ndbgx=10)
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:
c
c  iat:    maximum number of unique atom positions
c  natx:   maximum size of atomic cluster
c  ntitx:  maximum number of title lines
c  ndopx:  maximum number of dopants at any site
c  ngeomx: maximum number of lines written to geom.dat
c  neptx:  maximum number of energy points in dafs output files
c  maxln:  maximum number of lines written to feff.inp
c  nlogx:  number of logical parameters in logic array
c  nexafs: number of mcmaster paramters in exafs array
c  ndbgx:  maximum size of debugging code numbers = 2**ndbgx
c------------------------------------------------------------------------
c       include 'crystl.h'
c-*-fortran-*-

c  various parameters used by module crystl
      common /cryint/ iabs, iatom, ibasis, isystm, ispa, iperm, nsites,
     $            ipt(iat), idop(iat), imult(iat)
      save /cryint/

      parameter(nsysm=4, nshwrn=4)
      character*2  dopant(iat,ndopx)
      character*10 spcgrp, tag(iat)
      character*74 shwarn(nshwrn)
      character*77 sysmes(nsysm)
      common /crystr/ shwarn, sysmes, dopant, tag, spcgrp
      save /crystr/

      logical syserr, shift
      common /crylog/ syserr, shift
      save /crylog/

      dimension trmtx(3,3), st(iat,192,3)
      dimension cell(6), x(iat), y(iat), z(iat)
      dimension percnt(iat,ndopx)
      common /cryflt/ trmtx, st, cell, x, y, z, percnt
      save /cryflt/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:  (* = user input, % = error handling, ! = output needed to
c              construct cluster, the rest are used internally)
c
c * iabs:   index of absorber in unique coordinate list
c * iatom:  >=1 if atoms list is used, else =1
c * ibasis: =1 if basis list is used, else =0
c   isystm: index of crystal system (1..7)=(mono,orth,<not used>,tetr,
c           cubic,hex,triclinic)
c   ispa:   space group index, 1-230 from IXTC
c             0:       not recognized, error in input symbol
c             1-2:     triclinic
c             3-15:    monoclinic
c             16-74:   orthorhombic
c             75-142:  tetragonal
c             143-167: trigonal
c             168-194: hexagonal
c             195-230: cubic
c   iperm:  permutation index for non-standard settings, used in crystl
c             1:     default value -- no permutation necessary
c             1-6:   6 orthorhombic settings (abc, cab, bca, a-cb, ba-c, -cba)
c             11-12: 2 monoclinic settings, z-axis unique, y-axis unique
c             21-22: 2 tetragonal settings, standard and rotated
c ! ipt:    (iat) number of positions of unique atom in unit cell (1..192)
c   imult:  (iat) workspace for calculating multiplicities
c
c % sysmes: 4 line message if space group and axes/angles don't match
c % syserr: true if space group and axes/angles don't match
c % shwarn: 3 line message if space group may require a shift vector
c % shift:  true if space group may require a shift vector
c
c * dopant*2:  (iat,ndopx) matrix with all host and dopant atomic symbols
c * percnt:    (iat,ndopx) matrix with occupancies of hosts and dopants
c * tag*10:    (iat) character tag for each unique site in cell
c * spcgrp*10: space group symbol.  On output it is the short
c              Hermann-Maguin symbol in standard setting.  On input
c              spcgrp can be any short HM, Schoenflies, a number
c              between 1 and 230, or one of a small set of special
c              words (fcc, bcc, etc.).  Other symbol conventions
c              (full HM symbol, Shubnikov, 1935 ITXC, etc.) are not
c              and never will be used.
c
c * x,y,z:  (iat) arrays of fractional coordinates of unique positions
c                 in unit cell
c * cell:   (6) array of a,b,c,alpha,beta,gamma
c
c ! trmtx:  (3,3) transformation matrix between cell-axis and cartesian
c                 bases, see subroutine trans in clustr
c ! st:     (iat,192,3) fractional coordinates of all atoms in unit cell,
c                       first arg refers to unique atom list, second
c                       to position in cell, third is xyz.
c
c           fyi: 192 is the largest possible number of equivalent
c                positions in a cell of any symmetry. see, for example,
c                cubic f m 3 c.
c----------------------------------------------------------------------
      isp = ispa
      return
      end

      subroutine gtshft(shft,shw)
      logical shft
      character*74 shw(4)
c       include 'atparm.h'
c-*-fortran-*-
c  These parameters are the variable size declarations for the program
      parameter (iat=50, natx=800, ntitx=9, ndopx=4, ngeomx=natx)
      parameter (neptx=2**11, maxln=natx)
      parameter (nlogx=28, nexafs=13, ndbgx=10)
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:
c
c  iat:    maximum number of unique atom positions
c  natx:   maximum size of atomic cluster
c  ntitx:  maximum number of title lines
c  ndopx:  maximum number of dopants at any site
c  ngeomx: maximum number of lines written to geom.dat
c  neptx:  maximum number of energy points in dafs output files
c  maxln:  maximum number of lines written to feff.inp
c  nlogx:  number of logical parameters in logic array
c  nexafs: number of mcmaster paramters in exafs array
c  ndbgx:  maximum size of debugging code numbers = 2**ndbgx
c------------------------------------------------------------------------
c       include 'crystl.h'
c-*-fortran-*-

c  various parameters used by module crystl
      common /cryint/ iabs, iatom, ibasis, isystm, ispa, iperm, nsites,
     $            ipt(iat), idop(iat), imult(iat)
      save /cryint/

      parameter(nsysm=4, nshwrn=4)
      character*2  dopant(iat,ndopx)
      character*10 spcgrp, tag(iat)
      character*74 shwarn(nshwrn)
      character*77 sysmes(nsysm)
      common /crystr/ shwarn, sysmes, dopant, tag, spcgrp
      save /crystr/

      logical syserr, shift
      common /crylog/ syserr, shift
      save /crylog/

      dimension trmtx(3,3), st(iat,192,3)
      dimension cell(6), x(iat), y(iat), z(iat)
      dimension percnt(iat,ndopx)
      common /cryflt/ trmtx, st, cell, x, y, z, percnt
      save /cryflt/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:  (* = user input, % = error handling, ! = output needed to
c              construct cluster, the rest are used internally)
c
c * iabs:   index of absorber in unique coordinate list
c * iatom:  >=1 if atoms list is used, else =1
c * ibasis: =1 if basis list is used, else =0
c   isystm: index of crystal system (1..7)=(mono,orth,<not used>,tetr,
c           cubic,hex,triclinic)
c   ispa:   space group index, 1-230 from IXTC
c             0:       not recognized, error in input symbol
c             1-2:     triclinic
c             3-15:    monoclinic
c             16-74:   orthorhombic
c             75-142:  tetragonal
c             143-167: trigonal
c             168-194: hexagonal
c             195-230: cubic
c   iperm:  permutation index for non-standard settings, used in crystl
c             1:     default value -- no permutation necessary
c             1-6:   6 orthorhombic settings (abc, cab, bca, a-cb, ba-c, -cba)
c             11-12: 2 monoclinic settings, z-axis unique, y-axis unique
c             21-22: 2 tetragonal settings, standard and rotated
c ! ipt:    (iat) number of positions of unique atom in unit cell (1..192)
c   imult:  (iat) workspace for calculating multiplicities
c
c % sysmes: 4 line message if space group and axes/angles don't match
c % syserr: true if space group and axes/angles don't match
c % shwarn: 3 line message if space group may require a shift vector
c % shift:  true if space group may require a shift vector
c
c * dopant*2:  (iat,ndopx) matrix with all host and dopant atomic symbols
c * percnt:    (iat,ndopx) matrix with occupancies of hosts and dopants
c * tag*10:    (iat) character tag for each unique site in cell
c * spcgrp*10: space group symbol.  On output it is the short
c              Hermann-Maguin symbol in standard setting.  On input
c              spcgrp can be any short HM, Schoenflies, a number
c              between 1 and 230, or one of a small set of special
c              words (fcc, bcc, etc.).  Other symbol conventions
c              (full HM symbol, Shubnikov, 1935 ITXC, etc.) are not
c              and never will be used.
c
c * x,y,z:  (iat) arrays of fractional coordinates of unique positions
c                 in unit cell
c * cell:   (6) array of a,b,c,alpha,beta,gamma
c
c ! trmtx:  (3,3) transformation matrix between cell-axis and cartesian
c                 bases, see subroutine trans in clustr
c ! st:     (iat,192,3) fractional coordinates of all atoms in unit cell,
c                       first arg refers to unique atom list, second
c                       to position in cell, third is xyz.
c
c           fyi: 192 is the largest possible number of equivalent
c                positions in a cell of any symmetry. see, for example,
c                cubic f m 3 c.
c----------------------------------------------------------------------
      shft = shift
      shw(1) = shwarn(1)
      shw(2) = shwarn(2)
      shw(3) = shwarn(3)
      shw(4) = shwarn(4)
      return
      end
      subroutine atchck(iat,ndopx,core,dopant,edge,
     $                  iatom,ibasis,ispa,idop,
     $                  x,y,z,cell,dmax,gasses,qvect)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c--------------------------------------------------------------------
c input:
c   dopant: matrix of element symbols
c   x,y,z : arrays of unique atom fractional coordinates
c   cell:   array of cell constants
c   iatom:  number of unique atoms
c   ibasis: number of atoms in basis
c   dmax:   radius of desired cluster
c   ispa:   number between 1 and 230 denoting space group
c   pargon: percent argon in the i0 chamber
c   pnitro: percent nitrogen in the i0 chamber
c   pkrypt: percent krypton in the i0 chamber
c--------------------------------------------------------------------
c  check the consistancy of all the input parameters.
c  if any are funny write a run-time error message and
c  die gracefully.
c--------------------------------------------------------------------
c      parameter(iat=50,ndopx=4)
      parameter(epsi=0.001)
      character*2  el,dopant(iat,ndopx),dp,test,edge
      character*10 core
c      character*77 messg
      dimension    x(iat), y(iat), z(iat), cell(6), qvect(3)
      dimension    idop(iat), gasses(3)
      logical      ldie

      ldie  = .false.
      icent = 0
      inull = 0
      icore = 0
      test  = 'ab'
      iall  = iatom
      if (ibasis.gt.0) iall=ibasis

      do 100 i=1,iall
        el = dopant(i,1)
        call case(test,el)
        if ((el.eq.'nu').and.(i.ne.1)) inull = inull+1
        if (core.eq.el)                icore = icore+1

        if ((x(i).gt.1.0).or.(x(i).lt.-1.0)) then
            call messag(' ')
            call messag('Atom positions must be real numbers '//
     $                  'between -1 and 1.')
            ldie=.true.
        endif
        if ((y(i).gt.1.0).or.(y(i).lt.-1.0)) then
            call messag(' ')
            call messag('Atom positions must be real numbers '//
     $                  'between -1 and 1.')
            ldie=.true.
        endif
        if ((z(i).gt.1.0).or.(z(i).lt.-1.0)) then
            call messag(' ')
            call messag('Atom positions must be real numbers '//
     $                  'between -1 and 1.')
            ldie=.true.
        endif

        if ((is2z(el).eq.0).and.(el.ne.'nu')) then
            call messag(' ')
            call messag(dopant(i,1)//'???')
            call messag('One of your elements is not in the '//
     $                  'periodic table.')
            ldie=.true.
        endif
        if (core.eq.'nu') then
            call messag(' ')
            call messag('The core atom cannot be a null site.')
            ldie=.true.
        endif
 100  continue

c%%%        if (icent.ne.1) then
C%%%              call messag(' ')
c%%%            call messag('Feff requires one and only one absorption '//
c%%%       $                'site.')
c%%%            ldie=.true.
c%%%        endif

      if (inull.ge.1) then
          call messag(' ')
          call messag('Only one empty site is allowed and it '//
     $                'must be the first site listed.')
          ldie=.true.
      endif

c  check if core is a dopant
      do 104 i=1,iat
        do 102 j=2,idop(i)
          dp=dopant(i,j)
          call case(test,dp)
          if (core.eq.dp) icore=icore+1
 102    continue
 104  continue

      if (icore.eq.0) then
          call messag(' ')
          call messag('The central atom specified by the keyword '//
     $                '"core" is not')
          call messag('in the atom or dopant lists.')
          ldie=.true.
      endif

C%%%        if (icore.ge.2) then
C%%%            call messag(' ')
C%%%            call messag('Error reading the core atom.')
C%%%            ii=istrln(core)
C%%%            messg = core(:ii)//' appears more than once in the atom list.')
C%%%            call messag(messg)
C%%%            ldie=.true.
C%%%        endif

      if ((iatom.eq.0).and.(ibasis.eq.0)) then
          call messag(' ')
          call messag('You included no atoms in your input file.')
          ldie=.true.
      endif

      if (ispa.eq.0) then
          call messag(' ')
          call messag('Your space group does not exist.  Check '//
     $                'the international tables.')
          call messag('The ATOMS document explains adapting '//
     $                'notation for the keyboard.')
          ldie=.true.
      endif

      do 110 i=1,3
        if (cell(i).le.0) then
            call messag(' ')
            call messag('Cell constants cannot be negative or zero.')
            ldie=.true.
        endif
        if ((cell(i+3).lt.0).or.(cell(i+3)-180.0.gt.epsi)) then
            call messag(' ')
            call messag('Cell angles must be stated between '//
     $                  '0 and 180 degrees.')
            ldie=.true.
        endif
c%%%          if (qvect(i).lt.0) then
c%%%              call messag(' ')
c%%%              call messag('Values for the components of the q vector '//
c%%%       $              'in dafs must be non-negative.')
c%%%              ldie=.true.
c%%%          endif
110   continue

      if (dmax.le.0) then
          call messag(' ')
          call messag('Rmax is a radial distance.  It must be '//
     $                'positive.')
          ldie=.true.
      endif

      if ( (iatom.gt.1).and.(ibasis.gt.1) ) then
          call messag(' ')
          call messag('You may not specify both an atom list and a '//
     $                'basis list.')
          ldie=.true.
      endif

c       if ( (gasses(1)+gasses(3)+gasses(2))-1.0.gt.epsi ) then
c           call messag(' ')
c           call messag('The sum of the percentages of nitrogen, '//
c      $                'argon and krypton in the ')
c           call messag('I0 chamber cannot exceed 1.0')
c           ldie=.true.
c       endif
c
c       if ( (gasses(1).lt.-epsi).or.(gasses(3).lt.-epsi).or.
c      $            (gasses(3).lt.-epsi) ) then
c           call messag(' ')
c           call messag('The percentage of nitrogen, argon or '//
c      $         'krypton in the i0 chamber')
c           call messag('cannot be less than 0.')
c           ldie=.true.
c       endif

      call case(test,edge)
      if ( .not.((edge.eq.'k').or.(edge.eq.'l3').or.(edge.eq.'l2')
     $            .or.(edge.eq.'l1')) ) then
          call messag(' ')
          call messag('ATOMS only recognizes K and L edges.')
          ldie=.true.
      endif

      if (ldie) goto 666


      return

666   continue
      call messag('ATOMS cannot continue.  Please edit atoms.inp '//
     $            'and try again.')
      call messag(' ')
      stop

c end subroutine atchck
      end
      subroutine atinit(iat,natx,ntitx,ndopx,ngeomx,neptx,nlogx,
     $            title,elemnt,tag,noantg,spcgrp,edge,core,
     $            dopant,outfil,afname,geodat,refile,
     $            iatom,ibasis,dmax,ispa,iperm,idop,ngeom,iedge,iabs,
     $            ipt,iptful,isyst,nepts,nrefl,nnoan,
     $            x,y,z,cell,st,fullcl,atlis,
     $            percnt,gasses,qvect,egr,anot,
     $            logic,stdout)

c      $            lbasis,lindex,lfluor,ldwarf,lmm,lunit,lp1,
c      $            lself,li0,lgeom,ldafs,lf,lmod,lfeff,lfex,lgex,
c      $            lcrys, lmcm, lun, la0, lclus, lout, lrefl)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c-------------------------------------------------------
c  title   : user supplied comment line
c  elemnt  : character array of atom types
c  tag     : character array of site tags
c  spcgrp  : chararcter string with hermann-maguin space group designation
c  edge    : absorption edge of core atom, k or l3
c  x,y,z   : positions of atoms in cell-axis coordinates
c  cell    : lattice constants, a,b,c,alpha,beta & gamma
c  iatom   : number of unique atoms in cell
c  ibasis  : number of atoms in basis
c  dmax    : maximum distance in cluster
c  ispa    : number between 1 and 230, index of s.g. in int'l tables
c  iperm   : index of permutation matrix for non-standard setting
c  outfil  : output file name
c  refile  : reflection amplitudes file name
c  st:     : arrays containing positions in unit cell
c  isystm  : bravais lattice type
c  dopant  : dopant element symbol
c  percnt  : percent of dopant
c  gasses  : percent of argon, krypton, nitrogen in i0 chamber
c  logic   : logic flags, see arrays map
c-------------------------------------------------------
c  initialize everything there is to initialize
c-------------------------------------------------------
      parameter (zero=0.000000000, one=1.000000000)
c      parameter (iat=50,natx=800,neptx=2**11)
c      parameter (ntitx=9,ndopx=4,ngeomx=800)
      character*2  elemnt(iat),edge,dopant(iat,ndopx)
      character*10 spcgrp,tag(iat),core,geodat,noantg(iat)
      character*72 title(ntitx),outfil, afname, refile
      complex      anot(neptx)
      dimension    x(iat), y(iat), z(iat), cell(6), qvect(3)
      dimension    st(iat,192,3),
     $             fullcl(iat,192,3), atlis(natx,8)
      dimension    ipt(iat), iptful(iat), idop(iat),
     $             ngeom(ngeomx), nrefl(3)
      dimension    percnt(iat,ndopx), gasses(3)
      logical      logic(nlogx), stdout

c--------- characters -------------------------------------------------
      do 10 i=1,ntitx
        title(i) = ' '
10    continue
      spcgrp = ' '
      edge   = ' '
      core   = ' '
      outfil = 'feff.inp'
      call lower(outfil)
      afname = ' '
      geodat = 'geom.dat'
      call lower(geodat)
      refile = 'reflect.dat'
      call lower(refile)

c--------- logicals ---------------------------------------------------
      do 15 i=1,nlogx
        logic(i) = .false.
 15   continue
c     7: feff   19: feff8   28: mcm
      logic(7)  = .true.
      logic(19) = .false.
      logic(28) = .true.
c      stdout = .true.

c--------- reals ------------------------------------------------------
      dmax   = 5.
      do 17 i=1,3
        gasses(i) = zero
 17   continue
      egr    = zero

c--------- integers ---------------------------------------------------
      ispa   = 0
      iperm  = 1
      iatom  = 0
      ibasis = 0
      isyst  = 0
      iabs   = 0
      iedge  = 0
      nepts  = 0
      nnoan  = 0
      do 20 i=1,3
        nrefl(i) = 0
 20   continue

c--------- arrays ------------------------------------------------------
      do 100 i=1,iat
        elemnt(i) = ' '
        tag(i)    = ' '
        noantg(i) = ' '
        x(i)      = zero
        y(i)      = zero
        z(i)      = zero
        ipt(i)    = 0
        iptful(i) = 0
        idop(i)   = 1
        do 90 j=1,ndopx
          dopant(i,j) = ' '
          percnt(i,j) = zero
 90     continue
100   continue
      do 110 i=1,3
        cell(i)   = zero
        cell(i+3) = 90.
        qvect(i)  = zero
110   continue
      do 120 i=1,ngeomx
        ngeom(i) = 0
 120  continue

c  unit transformation for trivial position k=1
      do 180 i=1,3
        do 170 j=1,iat
          do 160 k=1,192
            st(j,k,i)     = zero
            fullcl(j,k,i) = zero
160       continue
170     continue
180   continue

      do 200 i=1,natx
        do 190 j=1,8
          atlis(i,j) = zero
190     continue
200   continue

      do 210 i=1,neptx
        anot(i) = cmplx(zero,zero)
 210  continue

      return
c end subroutine atinit
      end
      subroutine atinpt(iat,ntitx,ndopx,nlogx,
     $            title,ntit,iatom,ibasis,iabs,dmax,ispa,iperm,
     $            idop,iedge,nepts,nrefl,isystm,nnoan,
     $            elemnt,tag,noantg,edge,core,spcgrp,inpgrp,
     $            outfil,shwarn,afname,refile,dopant,
     $            x,y,z,cell,percnt,gasses,qvect,egr,
     $            logic, stdout, expnd, shift)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c-------------------------------------------------------
c  iatom   : number of unique atoms in cell
c  ibasis  : number of atoms in basis
c  ispa    : number between 1 and 230, index of s.g. in int'l tables
c  iperm   : index of permutation matrix for non-standard setting
c  title   : user supplied comment line
c  elemnt  : character array of atom types
c  tag     : character array of site tags
c  spcgrp  : chararcter string with hermann-maguin space group designation
c  edge    : absorption edge of core atom, k or l3
c  outfil  : output file name
c  refile  : reflection amplitude file name
c  dopant  : dopant element symbol
c  x,y,z   : positions of atoms in cell-axis coordinates
c  cell    : lattice constants, a,b,c,alpha,beta & gamma
c  dmax    : maximum distance in cluster
c  percent : percent of dopant
c  logic   : various flags, see arrays.map
c  noantg  : tags for which anomalous correction is to be neglected
c---------------------------------------------------------------------
c  this parses the lines of the command file looking for keywords then
c  reads in the value for that keyword from the next word.  the
c  structure is nearly free format except for the atom list which must
c  come at the end of the command file,  the reason for this is that
c  "b" and "c" are keywords and atomic symbols thus enforcing some
c  structure in the command file is the easiest way to distinguish them.
c---------------------------------------------------------------------

      parameter(ndpmax=10, nwdx=20)
c      parameter(iat=50, ntitx=9, ndopx=4 )
      parameter(zero=0)
      character*2  elemnt(iat),edge,dopant(iat,ndopx),
     $             doplis(ndpmax),test
      character*10 spcgrp,inpgrp,tag(iat),replcd(ndpmax),core,co,tagup
      character*10 noantg(iat)
      character*20 words(nwdx)
      character*72 title(ntitx), titln, outfil, fname, afname, refile
      character*80 string,toss,messg*78
      logical      logic(nlogx), stdout, expnd
      logical      inpnul, lcore, lshift, there
      dimension    x(iat), y(iat), z(iat), cell(6), qvect(3)
      dimension    idop(iat), nrefl(3), gasses(3)
      dimension    percnt(iat,ndopx), pclis(ndpmax)
      parameter(nshwrn=4)
      character*74 shwarn(nshwrn)
      logical      shift

 1410 format(bn,f10.0)
 1440 format(a)
 1450 format(a,i2)
 4000 format(i2)
 4010 format(' *** Notice at line ',i3)
 4020 format(' *** Warning at line ',i3)

c---------------------------------------------------------------------
c  initialize some things used only in this routine
      inpnul = .true.
      lcore  = .false.
      lshift = .false.
      icore  = 0
      ndop   = 0
      ntit   = 0
      xshift = zero
      yshift = zero
      zshift = zero
      nnoan  = 0
      test   = 'ab'
      nline  = 0
      ierr   = 0
      iertot = 0
c  the value of the variable test must be in the same case as the
c  keyword names in the long block of elseif's

c---------------------------------------------------------------------
c  open atom.inp, look for upper and lower case file names for both
c  atom.inp and atoms.inp
      if (.not.stdout) then
          fname = 'atoms.inp'
          call lower(fname)
          inquire(file=fname,exist=there)
          if (.not.there) then
              call upper(fname)
              inquire(file=fname,exist=there)
          endif
          if (.not.there) then
              fname = 'atom.inp'
              call lower(fname)
              inquire(file=fname,exist=there)
              if (.not.there) then
                  call upper(fname)
                  inquire(file=fname,exist=there)
              endif
          endif
          if (.not.there) then
              call messag('Input file for ATOMS is not found. '//
     $                    'Hasta luego.')
              stop
          endif
          open(unit=1,file=fname,status='old')
      endif

c---------------------------------------------------------------------
c  begin reading the input file, words must be cleared each time to
c  avoid unintentionally labling more than one atom type as the
c  central atom.
 101  continue
      if (stdout) then
          read (*,1440,end=191,err=191)string
      else
          read (1,1440,end=191,err=191)string
      endif
      nline = nline+1
      call untab(string)
c       call uncomm(string)

 120  continue
      nwds = nwdx
      do 122 iw=1,nwds
        words(iw)=' '
 122  continue
      call triml(string)
c                                           - denotes end of job
      if  (string(1:1).eq.'-') goto 191
c                                           skip a comment line
      if  ((string(1:1).eq.'!').or.(string(1:1).eq.'*')
     $ .or.(string(1:1).eq.'%').or.(string(1:1).eq.'#')) goto 101
c                                           skip a blank line
      if  (string.eq.' ') goto 101
c                                           begin reading line
      i=1
      call bwords(string,nwds,words)
      inpnul = .false.

c  ******************** input file parsing *************************
c  read a word, identify it, assign the value from the following word(s)
c  increment i and come back.  i points to position in string, when i
c  exceeds nwds go read a new line.
130   continue
      call case(test,words(i))
c                                           skip a blank line
      if     (words(i).eq.' ') then
          goto 101
c                                           ignore everything after !,*,%
      elseif ((words(i)(1:1).eq.'!').or.(words(i)(1:1).eq.'*').or.
     $        (words(i)(1:1).eq.'%').or.(words(i)(1:1).eq.'#')) then
          goto 101
      elseif (words(i).eq.'end') then
          goto 191
c                                           title and comment are synonyms
      elseif ((words(i).eq.'comment').or.(words(i).eq.'title')) then
          if (ntit.lt.ntitx) then
              call gettit(words(i), string, titln, ntit, stdout)
              title(ntit) = titln
          endif
          goto 101
c                                           read the next ten characters
c                                           into space, continue reading
c                                           rest of the line.
c                                           this one is perverse and must
c                                           be handled specially
      elseif (words(i).eq.'space') then
          toss = string
          call case(test,toss)
          m = index(toss, 'space')
          toss = string(m+6:)
          call triml(toss)
          string = toss
          if (string(:1).eq.'=') toss = string(2:)
          call triml(toss)
          spcgrp = toss(1:10)
          call triml(spcgrp)
          call case(test,spcgrp)
          string = toss(11:)
          inpgrp = spcgrp
          goto 120
c                                           outfile, default=feff.inp
      elseif (words(i)(1:3).eq.'out') then
          outfil=words(i+1)
          i=i+2
c                                           specified edge, default by z
      elseif ((words(i).eq.'edge').or.(words(i).eq.'hole'))  then
          edge=words(i+1)
          call case(test,edge)
          i=i+2
c                                           specified core
      elseif ((words(i).eq.'core').or.(words(i)(1:4).eq.'cent')) then
          core=words(i+1)
          call case(test,core)
          lcore = .true.
          i=i+2
c                                           dopants
      elseif (words(i)(1:3).eq.'dop') then
          ndop = ndop + 1
          if (ndop.gt.ndopx) then
              write(messg, 4020)nline
              call messag(messg)
              call messag('     You have exceeded the '//
     $                    'maximum number of dopants.')
              call messag('     ATOMS will ignore this and all '//
     $                    'further dopants.-')
              ierr = 2
              goto 137
          endif
          doplis(ndop) = words(i+1)
          call case(test,doplis(ndop))
          replcd(ndop) = words(i+2)
          call case(test,replcd(ndop))
          call getrea(words(i), words(i+3), pclis(ndop),
     $                nline, ierr)
 137      continue
          i=i+4
c                                           argon, fluorescence
      elseif (words(i)(1:3).eq.'arg') then
          call getrea(words(i), words(i+1), gasses(1), nline, ierr)
          logic(3) = .true.
          i=i+2
c                                           krypton, fluorescence
      elseif (words(i)(1:3).eq.'kry') then
          call getrea(words(i), words(i+1), gasses(2), nline, ierr)
          logic(3) = .true.
          i=i+2
c                                           nitrogen, fluorescence
      elseif (words(i)(1:3).eq.'nit') then
          call getrea(words(i), words(i+1), gasses(3), nline, ierr)
          logic(3) = .true.
          i=i+2
c                                           flag for indexing
      elseif (words(i)(1:3).eq.'ind') then
          call getlgc(words(i), words(i+1), logic(4), nline, ierr)
          i=i+2
c                                           geom.dat
      elseif (words(i)(1:3).eq.'geo') then
          call getlgc(words(i), words(i+1), logic(6), nline, ierr)
          i=i+2
c                                           xanes keywords (& not feff8)
c       elseif (words(i)(1:3).eq.'xan') then
c           call getlgc(words(i), words(i+1), logic(18), nline, ierr)
c           logic(19) = .false.
c           i=i+2
c                                           feff8 keywords (& not xanes)
      elseif (words(i)(1:5).eq.'feff8') then
          call getlgc(words(i), words(i+1), logic(19), nline, ierr)
          logic(18) = .false.
          i=i+2
c                                           run clustr & output (note order)
      elseif (words(i)(1:4).eq.'feff') then
          call getlgc(words(i), words(i+1), logic(7), nline, ierr)
          i=i+2
c                                           calc and write McM corrs.
      elseif (words(i)(1:4).eq.'corr') then
          call getlgc(words(i), words(i+1), logic(28), nline, ierr)
          i=i+2
c                                           heh, heh, heh!
      elseif (words(i)(1:3).eq.'dwa') then
          call getlgc(words(i), words(i+1), logic(5), nline, ierr)
          i=i+2
c * * * * * * * * * * diagnostic functions * * * * * * * * * * * *
c                                           mcmast.dat
      elseif (words(i)(1:3).eq.'mcm') then
          call getlgc(words(i), words(i+1), logic(14), nline, ierr)
          i=i+2
c                                           self.dat
      elseif (words(i).eq.'self') then
          call getlgc(words(i), words(i+1), logic(15), nline, ierr)
          i=i+2
c                                           i0.dat
      elseif (words(i).eq.'i0') then
          call getlgc(words(i), words(i+1), logic(16), nline, ierr)
          i=i+2
c                                           unit.dat
      elseif (words(i)(1:3).eq.'uni') then
          call getlgc(words(i), words(i+1), logic(9), nline, ierr)
          i=i+2
c                                           p1.inp
      elseif (words(i).eq.'p1') then
          call getlgc(words(i), words(i+1), logic(8), nline, ierr)
          i=i+2
c                                           f.dat, diagnostic for f'/"
      elseif (words(i)(1:4).eq.'fdat') then
          call getlgc(words(i), words(i+1), logic(17), nline, ierr)
          i=i+2
c                                           print module messages
      elseif (words(i)(1:3).eq.'mod') then
          call getlgc(words(i), words(i+1), logic(20), nline, ierr)
          i=i+2
c                                           print location messages
      elseif (words(i).eq.'message') then
          call getint(words(1),words(i+1), imess, nline, ierr)
          if ((imess.eq.0).or.(imess.eq.2)) logic(21) = .true.
          if ((imess.eq.0).or.(imess.eq.3)) logic(22) = .true.
          if ((imess.eq.0).or.(imess.eq.4)) logic(23) = .true.
          if ((imess.eq.0).or.(imess.eq.5)) logic(24) = .true.
          if ((imess.eq.0).or.(imess.eq.6)) logic(25) = .true.
          if ((imess.eq.0).or.(imess.eq.7)) logic(26) = .true.
          if (imess.eq.0) logic(20) = .true.
          i=i+2
c * * * * * * * * end diagnostic functions * * * * * * * * * * * *
c                                           rmax, default=5.0
      elseif (words(i)(1:3).eq.'rma') then
          call getrea(words(i), words(i+1), dmax, nline, ierr)
          i=i+2
c                                           the lattice constants
      elseif (words(i).eq.'a') then
          call getrea(words(i), words(i+1), cell(1), nline, ierr)
          i=i+2
      elseif (words(i).eq.'b') then
          call getrea(words(i), words(i+1), cell(2), nline, ierr)
          i=i+2
      elseif (words(i).eq.'c') then
          call getrea(words(i), words(i+1), cell(3), nline, ierr)
          i=i+2
c                                           the latice angles
      elseif (words(i)(1:3).eq.'alp') then
          call getrea(words(i), words(i+1), cell(4), nline, ierr)
          i=i+2
      elseif (words(i)(1:3).eq.'bet') then
          call getrea(words(i), words(i+1), cell(5), nline, ierr)
          i=i+2
      elseif (words(i)(1:3).eq.'gam') then
          call getrea(words(i), words(i+1), cell(6), nline, ierr)
          i=i+2
c                                           shift vector
      elseif (words(i).eq.'shift') then
          call getrea(words(i), words(i+1), xshift, nline, ierr)
          call getrea(words(i), words(i+2), yshift, nline, ierr)
          call getrea(words(i), words(i+3), zshift, nline, ierr)
          lshift = .true.
          i=i+4

c     ************ DAFS STUFF ****************

c                                           q vector for dafs
      elseif ((words(i)(1:4).eq.'qvec').or.(words(i).eq.'dafs')) then
          call getrea(words(i), words(i+1), qvect(1), nline, ierr)
          call getrea(words(i), words(i+2), qvect(2), nline, ierr)
          call getrea(words(i), words(i+3), qvect(3), nline, ierr)
          logic(10) = .true.
          i=i+4
      elseif (words(i).eq.'feout') then
          afname = words(i+1)
          i=i+2
c                                           reflection amplitudes
      elseif (words(i)(1:4).eq.'refl') then
          call getint(words(i), words(i+1), nrefl(1), nline, ierr)
          call getint(words(i), words(i+2), nrefl(2), nline, ierr)
          call getint(words(i), words(i+3), nrefl(3), nline, ierr)
          logic(11) = .true.
          i=i+4
c                                           reflection amplitude file name
      elseif (words(i).eq.'refile') then
          refile = words(i+1)
          i=i+2
c                                           # of grid points for dafs
      elseif (words(i)(1:3).eq.'nep') then
          call getint(words(i), words(i+1), nepts, nline, ierr)
          i=i+2
c                                           grid spacing for dafs
      elseif (words(i)(1:3).eq.'egr') then
          call getrea(words(i), words(i+1), egr, nline, ierr)
          i=i+2
c                                           neglect anomalous correction
      elseif (words(i)(1:4).eq.'noan') then
          nnoan = nnoan+1
          noantg(nnoan) =  words(i+1)
          i=i+2

c                                           beneath the word atom is a
c                                           5 col. list of unique atoms
c                                           in this order:
c                                              sym  x  y  z   center?
c                                           atom info is read in until eof
c                                           or - is found
      elseif (words(i)(1:3).eq.'ato') then
 140      continue
          if (stdout) then
              read (*,1440,end=191,err=191)string
          else
              read (1,1440,end=191,err=191)string
          endif
          nline = nline+1
          call untab(string)
          call triml(string)
          if (string(1:1).eq.'-') goto 191
          if ( (string(1:1).eq.' ').or.(string(1:1).eq.'!').or.
     $         (string(1:1).eq.'*').or.(string(1:1).eq.'%').or.
     $         (string(1:1).eq.'#'))
     $                goto 140
          iatom=iatom+1
          if (iatom.gt.iat) then
              ierr = 3
 400          format(' *** Error: Your atoms list exceeds ',i3,
     $                    ', the hardwired limit.')
              write(messg, 400)iat
              call messag(' ')
              call messag(messg)
              call messag('     Reset iat in the source code then ')
              call messag('     recompile ATOMS to handle a larger '//
     $                    'atom list.-')
              goto 191
          endif
          if (expnd) then
              call messag(' *** Error: No expanded atoms list yet.-')
              stop
          else
              call getatm(nline, ierr, string, elemnt(iatom),
     $                    x(iatom), y(iatom), z(iatom),
     $                    tag(iatom))
              if (ierr.eq.2) iatom=iatom-1
          endif
          if (ierr.eq.3) iertot = 1
          ierr = 0
          goto 140
c                                           handle basis
c                                           iatom set to one, one atom
c                                           will be expanded as a
c                                           point and rest will be added
      elseif (words(i)(1:3).eq.'bas') then
          logic(2) = .true.
          iatom  = iatom+1
 155      continue
          if (stdout) then
              read (*,1440,end=191,err=191)string
          else
              read (1,1440,end=191,err=191)string
          endif
          nline = nline+1
          call untab(string)
          call triml(string)
          if (string(1:1).eq.'-') goto 191
          if ( (string(1:1).eq.' ').or.(string(1:1).eq.'!').or.
     $         (string(1:1).eq.'*').or.(string(1:1).eq.'%').or.
     $         (string(1:1).eq.'#'))
     $                goto 155
          ibasis=ibasis+1
          if (ibasis.gt.iat) then
              ierr = 3
 410          format(' *** Error: Your basis list exceeds ',i3,
     $                    ', the hardwired limit.')
              write(messg, 400)iat
              call messag(' ')
              call messag(messg)
              call messag('     Reset iat in the source code then ')
              call messag('     recompile ATOMS to handle a larger '//
     $                    'atom list.-')
              goto 191
          endif
          if (expnd) then
              call messag(' *** Error: No expanded atoms list yet.-')
              stop
          else
              call getatm(nline, ierr, string, elemnt(ibasis),
     $                    x(ibasis), y(ibasis), z(ibasis),
     $                    tag(ibasis))
              if (ierr.eq.2) iatom=iatom-1
          endif
          if (ierr.eq.3) iertot = 1
          ierr = 0
          goto 155
      else
          write(messg, 4020)nline
          call messag(messg)
          iunk  = istrln(words(i))
          messg = '     "'//words(i)(:iunk)//'" is an unknown keyword.-'
          call messag(messg)
          ierr = 2
          i=i+2
      endif
c     if read entire line then read next line else read next word in line
      if (ierr.eq.3) iertot = 1
      ierr = 0
      if (i.ge.nwds) goto 101
      goto 130

c     done reading lines
191   continue
      if (ierr.eq.3) iertot = 1
      if (iertot.ne.0) then
          call messag(' ')
          call messag(' *** Ending early due to faulty input file.-')
          call messag(' ')
          ierr = 1
          stop
      endif
      if (inpnul) then
          logic(1)=.true.
          goto 300
      endif

c--------------------------------------------------------------------
c----- do a few things with the keyword values before leaving -------
c--------------------------------------------------------------------

c  turn off all messages if program compiled for standard in/output
c  also disable dafs, geom.dat, unit.dat, p1.dat outputs
      if (stdout) then
          logic(6)  = .false.
          logic(8)  = .false.
          logic(9)  = .false.
          logic(10) = .false.
          logic(11) = .false.
          logic(14) = .false.
          logic(15) = .false.
          logic(16) = .false.
          logic(17) = .false.
          logic(20) = .false.
          logic(21) = .false.
          logic(22) = .false.
          logic(23) = .false.
          logic(24) = .false.
          logic(25) = .false.
          logic(26) = .false.
      endif

c  iall is the number of atom in the atom or basis list
      iall = iatom
      if (logic(2)) iall = ibasis

c  add shift vector to unique atom coordinates
      do 250 i=1,iall
        x(i) = x(i) + xshift
        y(i) = y(i) + yshift
        z(i) = z(i) + zshift
 250  continue

c  identify space group from supplied symbol
      call groups

c  organize matrices containing site contents
      call dopfix(iat,ndopx,ndpmax,
     $            iall,ndop,tag,doplis,replcd,pclis,
     $            elemnt,dopant,percnt,idop)

c  try to determine the atomic symbol of the core atom
c  it need not be specified if the atom list has one site in it...
      if ((iatom.eq.1).and.(ibasis.le.1).and.(.not.lcore))  then
          core  = elemnt(1)
          iabs  = 1
          call case(test,core)
          lcore = .true.
      endif

c      print*,'initial value of core: ', core, ' iabs=', iabs
c  search through tag list for the specified core...
      co   = core
      call case(test,co)
      do 200 i=1,iall
        tagup=tag(i)
        call case(test,tagup)
        if (tagup.eq.co) then
            icore = icore + 1
            core  = elemnt(i)
            iabs  = i
            lcore = .true.
        endif
 200  continue
 210  continue
c      print*,'after searching tag list: ', core, ' iabs=', iabs

c  if the specified core was not found in the tag list, search the dopants
c  (use the variable name tagup to avoid defining another variable)
      if (iabs.eq.0) then
          do 280 i=1,iall
            do 270 j=2,idop(i)
              co   = core
              call case(test,co)
              tagup=dopant(i,j)
              call case(test,tagup)
              if (tagup.eq.co) then
                  iabs  = i
                  icore = icore+1
              endif
 270        continue
 280      continue
      endif

c----------------------------------------------------------------------
c  choose l3 (=4) or k (=1) edge, the numbers are chosen to suit mucal
      if (edge.eq.'k') then
          iedge = 1
      elseif (edge.eq.'l1') then
          iedge = 2
      elseif (edge.eq.'l2') then
          iedge = 3
      elseif (edge.eq.'l3') then
          iedge = 4
      elseif (edge.eq.' ') then
          if (is2z(core).gt.57) then
              edge  = 'l3'
              iedge = 4
          else
              edge  = 'k'
              iedge = 1
          endif
      endif

c  check if core has been specified, this is to satisfy backwards
c  incompatibility concerns
      if (.not.lcore) then
          call messag(' ')
          call messag(' *** Error:  while reading atom or basis list')
          call messag('     In this version of ATOMS, the central '//
     $                'atom is specified by the keyword "core"')
          call messag('     The fifth column of the atom list is '//
     $                'reserved for the site tag.')
          call messag('Please edit the input file and run atoms '//
     $                'again.-')
          call messag(' ')
          stop
      endif

c  check if 0 sites match core
      if (icore.eq.0) then
          ii = istrln(co)
          call messag(' ')
          call messag(' *** Error:  while reading central atom.')
          call messag('     '//co(:ii)//' matches no sites.-')
          call messag(' ')
          stop
      endif

c  check if 2 or more sites match core
      if (icore.ge.2) then
          ii = istrln(co)
          call messag(' ')
          call messag(' *** Error:  while reading central atom.')
          call messag('     '//co(:ii)//' matches more than one site.-')
          call messag(' ')
          stop
      endif

c  check if outfil or afname will overwrite input file
      if ((outfil.eq.fname).or.(afname.eq.fname).or.
     $            (refile.eq.fname)) then
          call messag(' *** Error: ')
          call messag('     One of your specified output files names '//
     $                'will overwrite the input file.')
          call messag('     This is not allowed.-')
          call messag(' ')
          stop
      endif

 300  continue
      return

c 666  continue
c      stop

c end subroutine atinpt
      end
      subroutine dopfix(iat,ndopx,ndpmax,
     $                  iatom,ndop,tag,doplis,replcd,pclis,
     $                  elemnt,dopant,percnt,idop)

c-----------------------------------------------------------------------
c  inputs:
c    iat:    number of unique crystallographic sites
c    ndop:   total number of dopants
c    tag:    array of site tags
c    doplis: array of all doping atoms
c    replcd: array of all replaced atoms
c    pclis:  array of doping percentages
c    elemnt: array of atoms from atom list
c  output:
c    dopant: iat x ndopx matrix of doping atoms indexed by site
c    percnt: iat x ndopx matrix of doping percentages indexed by site
c    idop:   array indicating number of dopants at each site
c-----------------------------------------------------------------------
      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

c      parameter(iat=50, ndopx=4, ndpmax=10)

      character*2  test
      character*2  doplis(ndpmax), dopant(iat,ndopx), elemnt(iat)
      character*10 tag(iat), tagup, replcd(ndpmax), repup
      dimension    pclis(ndpmax),  percnt(iat,ndopx)
      dimension    idop(iat)

      test = 'ab'

      do 100 i=1,ndop
        repup = replcd(i)
        call case(test,repup)
        do 50 j=1,iatom
          tagup = tag(j)
          call case(test,tagup)
          if (repup.eq.tagup) then
              idop(j) = idop(j) + 1
              dopant(j,idop(j)) = doplis(i)
              percnt(j,idop(j)) = pclis(i)
          endif
 50     continue
 100  continue

      do 200 i=1,iatom
        dopant(i,1) = elemnt(i)
        sumdop = 0
        do 150 j=2,idop(i)
          sumdop = sumdop + percnt(i,j)
 150    continue
        percnt(i,1) = 1. - sumdop
 200  continue

      return
c  end subroutine dopfix
      end

      subroutine getatm(nline, ierr, string, elem, x, y, z, tag)
      implicit real(a-h,o-z)
      implicit integer(i-n)
c--------------------------------------------------------------
c  copyright (c) 1998 Bruce Ravel
c--------------------------------------------------------------
c========================================================================
c  This routine parses a line containing information about a unique atom
c  position in a unit cell and returns the information about that site.
c========================================================================
      parameter(nwdx=20)
      character*2   elem, test, check*10
      character*20  tag*10, words(nwdx)
      character*78  messg
      character*(*) string
      integer       nline, ierr

      test = 'ab'
      do 10 i=1,nwdx
        words(i) = ' '
 10   continue
      nwds = nwdx
      call bwords(string,nwds,words)

      if (nwds.lt.4) goto 666
      check = words(1)
      call lower(check)
      if (check(1:2) .eq. 'nu') goto 20
      if (is2z(words(1)).eq.0) goto 666
      if (istrln(words(1)).gt.2) goto 666
 20   continue
c      if ( (nwds.lt.4).or.(is2z(words(1)).eq.0).or.
c     $            (istrln(words(1)).gt.2) ) goto 666


c  get element symbol, x, y, z, tag in that order
      elem = words(1)
      call case(test,elem)
      call getrea('atomic position', words(2), x, nline, ierr)
      call getrea('atomic position', words(3), y, nline, ierr)
      call getrea('atomic position', words(4), z, nline, ierr)
      if ((nwds.gt.4) .and. ((words(5).ne.'!').and.
     $            (words(5).ne.'%').and.(words(5).ne.'*')) ) then
          tag = words(5)
      else
          tag = words(1)
          call fixsym(tag(1:2))
      endif

      return

 666  continue
      call messag(' ')
      call messag(' ')
 400  format(' *** Warning at line ',i3, ' while reading the atom',
     $            ' or basis list')
      write(messg,400)nline
      call messag(messg)
      call messag(' ')
      call messag('     The atom list is a formatted list, and must '//
     $            'be at the end of')
      call messag('        the input file.')
      call messag('     The first column is a two character '//
     $            'elemental symbol.')
      call messag('     The next three are coordinates in the unit '//
     $            'cell.')
      call messag('     The fifth column is the optional site tag.')
      call messag('     Atoms should be able to continue, but may '//
     $            'not produce')
      call messag('        the correct output.')
      call messag('     You should edit atoms.inp and try again.-')
      call messag(' ')
      ierr = 2

c  end subroutine getatm.f
      end
      subroutine crystl(idebug, ierr)
c=====================================================================
c  atoms module 2:  decode space group, determine unit cell contents
c=====================================================================
c  this module consists of the following subroutines and functions:
c     crystl basfil basort equipt genmul ibravl metric multip syschk
c  and requires
c     case dbglvl istrln messag positn upper
c
c  most variables are passed in common in crystl.h
c  idebug: an integer denoting the debug level, this is interpreted
c          into an array of binary bits which are used as logical flags.
c          multiple debugging features can be enables by specifying a
c          sum of bits
c     0:  disable all debuging function
c     1:  enable positional run-time messages
c  ierr:   output error code (0=no problem, 1=info, 2=warning, 3=error)
c------------------------------------------------------------------------
      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

      parameter(zero=0.e0, one=1.e0)
c       include 'atparm.h'
c-*-fortran-*-
c  These parameters are the variable size declarations for the program
      parameter (iat=50, natx=800, ntitx=9, ndopx=4, ngeomx=natx)
      parameter (neptx=2**11, maxln=natx)
      parameter (nlogx=28, nexafs=13, ndbgx=10)
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:
c
c  iat:    maximum number of unique atom positions
c  natx:   maximum size of atomic cluster
c  ntitx:  maximum number of title lines
c  ndopx:  maximum number of dopants at any site
c  ngeomx: maximum number of lines written to geom.dat
c  neptx:  maximum number of energy points in dafs output files
c  maxln:  maximum number of lines written to feff.inp
c  nlogx:  number of logical parameters in logic array
c  nexafs: number of mcmaster paramters in exafs array
c  ndbgx:  maximum size of debugging code numbers = 2**ndbgx
c------------------------------------------------------------------------
c       include 'crystl.h'
c-*-fortran-*-

c  various parameters used by module crystl
      common /cryint/ iabs, iatom, ibasis, isystm, ispa, iperm, nsites,
     $            ipt(iat), idop(iat), imult(iat)
      save /cryint/

      parameter(nsysm=4, nshwrn=4)
      character*2  dopant(iat,ndopx)
      character*10 spcgrp, tag(iat)
      character*74 shwarn(nshwrn)
      character*77 sysmes(nsysm)
      common /crystr/ shwarn, sysmes, dopant, tag, spcgrp
      save /crystr/

      logical syserr, shift
      common /crylog/ syserr, shift
      save /crylog/

      dimension trmtx(3,3), st(iat,192,3)
      dimension cell(6), x(iat), y(iat), z(iat)
      dimension percnt(iat,ndopx)
      common /cryflt/ trmtx, st, cell, x, y, z, percnt
      save /cryflt/
c -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
c  GLOSSARY:  (* = user input, % = error handling, ! = output needed to
c              construct cluster, the rest are used internally)
c
c * iabs:   index of absorber in unique coordinate list
c * iatom:  >=1 if atoms list is used, else =1
c * ibasis: =1 if basis list is used, else =0
c   isystm: index of crystal system (1..7)=(mono,orth,<not used>,tetr,
c           cubic,hex,triclinic)
c   ispa:   space group index, 1-230 from IXTC
c             0:       not recognized, error in input symbol
c             1-2:     triclinic
c             3-15:    monoclinic
c             16-74:   orthorhombic
c             75-142:  tetragonal
c             143-167: trigonal
c             168-194: hexagonal
c             195-230: cubic
c   iperm:  permutation index for non-standard settings, used in crystl
c             1:     default value -- no permutation necessary
c             1-6:   6 orthorhombic settings (abc, cab, bca, a-cb, ba-c, -cba)
c             11-12: 2 monoclinic settings, z-axis unique, y-axis unique
c             21-22: 2 tetragonal settings, standard and rotated
c ! ipt:    (iat) number of positions of unique atom in unit cell (1..192)
c   imult:  (iat) workspace for calculating multiplicities
c
c % sysmes: 4 line message if space group and axes/angles don't match
c % syserr: true if space group and axes/angles don't match
c % shwarn: 3 line message if space group may require a shift vector
c % shift:  true if space group may require a shift vector
c
c * dopant*2:  (iat,ndopx) matrix with all host and dopant atomic symbols
c * percnt:    (iat,ndopx) matrix with occupancies of hosts and dopants
c * tag*10:    (iat) character tag for each unique site in cell
c * spcgrp*10: space group symbol.  On output it is the short
c              Hermann-Maguin symbol in standard setting.  On input
c              spcgrp can be any short HM, Schoenflies, a number
c              between 1 and 230, or one of a small set of special
c              words (fcc, bcc, etc.).  Other symbol conventions
c              (full HM symbol, Shubnikov, 1935 ITXC, etc.) are not
c              and never will be used.
c
c * x,y,z:  (iat) arrays of fractional coordinates of unique positions
c                 in unit cell
c * cell:   (6) array of a,b,c,alpha,beta,gamma
c
c ! trmtx:  (3,3) transformation matrix between cell-axis and cartesian
c                 bases, see subroutine trans in clustr
c ! st:     (iat,192,3) fractional coordinates of all atoms in unit cell,
c                       first arg refers to unique atom list, second
c                       to position in cell, third is xyz.
c
c           fyi: 192 is the largest possible number of equivalent
c                positions in a cell of any symmetry. see, for example,
c                cubic f m 3 c.
c----------------------------------------------------------------------
      character*1  csymbr
      logical      lbasis
      dimension    ts(3,24), fs(3,3,24)
      integer      idebug, idbg(0:ndbgx)
      character*78 module, string

      lbasis=.false.
      if (ibasis.gt.0) lbasis=.true.
      call dbglvl(idebug,idbg)
c       print*,'idebug=',idebug,' idbg=(',idbg,')'
      module = 'crystl'
      ierr = 0

c  initialize centrosymmetry flag
c  isymce=1 marks a centrosymmetric atom, this is determined in
c  equipt and used elsewhere
      isymce = 0
      ns     = 0
c  initialize the transformation matrices
      do 30 i=1,3
        do 20 j=1,24
          ts(i,j) = zero
          do 10 k=1,3
            fs(i,k,j) = zero
 10       continue
 20     continue
 30   continue
      fs(1,1,1) = one
      fs(2,2,1) = one
      fs(3,3,1) = one

c------------------------------------------------------------
c  reorganize basis list so absorber is at symmetry point in cell
      nat = iatom
      if (lbasis.and.(iabs.ne.1)) then
          string = 'sorting basis'
          if (idbg(0).gt.0) call positn(module, string)
          call basort(iat,ndopx,ibasis,iabs,x,y,z,dopant,tag)
          nat = ibasis
      endif

c------------------------------------------------------------
c  permute atoms to standard setting
      if (iperm.gt.1) call fperm(iat, nat, iperm, cell, x, y, z)

c------------------------------------------------------------
c  calculate the transformation matrix
      string = 'getting matrix'
      if (idbg(0).gt.0) call positn(module, string)
      call metric(cell, trmtx)

c------------------------------------------------------------
c  call the routine that does all the group theory to find
c  individual atom positions from the symmetry properties of
c  the space group and the positions of the unique atoms within
c  the primitive cell
      string = 'entering equipt'
      if (idbg(0).gt.0) call positn(module, string)
      call equipt(isyst, isymce, csymbr, ns, ts, fs, spcgrp)

c------------------------------------------------------------
c  two independent checks on system of crystal
      string = 'checking system'
      if (idbg(0).gt.0) call positn(module, string)
      syserr = .false.
      if (isyst.ne.isystm) then
          call syschk(isyst, isystm, spcgrp, sysmes)
          syserr = .true.
          ierr = 1
      endif

c------------------------------------------------------------
c  get number of bravais lattice for use in multip
      ibravl=ibrav(csymbr)

c------------------------------------------------------------
c  get igen, the multiplicity of the general position then get
c  multiplicities of atom positions,
c  for bases iatom = 1, want to run multip only on first atom in basis
      string = 'calling genmul'
      if (idbg(0).gt.0) call positn(module, string)
      call genmul(ns,isymce,ibravl,igen)
      string = 'callinp multip'
      if (idbg(0).gt.0) call positn(module, string)
      call multip(iatom, ibravl, x, y, z, tag, fs, ts, isymce,
     $            ns, igen, cell, st, ipt, imult, ierr)

c------------------------------------------------------------
c  fill in basis at each point
      if (lbasis) then
          string = 'filling basis'
          if (idbg(0).gt.0) call positn(module, string)
          call basfil(iat, ibasis, x, y, z, ipt, st)
      endif

c------------------------------------------------------------
c  permute atoms back to non-standard setting, but not tetragonal
      if ((iperm.ne.22).and.(iperm.gt.1))
     $            call bperm(iat, nat, iperm, cell, ipt, st)

      return
c  end of module crystl
      end
      subroutine basfil(iat,ibasis,x,y,z,ipt,st)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c----------------------------------------------------------------
c  after the first atom in a basis list has been expanded according
c  to the group summetries, fill in the remainder of the atoms in
c  the basis at each point in the unit cell
c----------------------------------------------------------------
c  input:
c    iat:     parameter set in calling program
c    ibasis:  number of atoms in basis
c    x,y,z:   (iat) fractional coordinates of representative atoms
c  i/o:
c    ipt:     (iat) multiplicity of each site (all equal for a basis)
c    st:      (iat,192,3) fractional coordinates of all atoms in unit cell
c----------------------------------------------------------------
c      parameter (iat=50)

      dimension x(iat), y(iat), z(iat), st(iat,192,3), ipt(iat)

      do 50 i=2,ibasis
        do 40 j=1,ipt(1)
          st(i,j,1) = st(1,j,1) + x(i)
          st(i,j,2) = st(1,j,2) + y(i)
          st(i,j,3) = st(1,j,3) + z(i)
 40     continue
        ipt(i) = ipt(1)
 50   continue

      return
c  end subroutine basfil
      end

      subroutine basort(iat,ndopx,ibasis,iabs,x,y,z,dopant,tag)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c------------------------------------------------------------------
c  want the absorber at the top of the unique atom list in the
c  basis calculation.  this routine puts the absorber on the point
c  of translation symmetry and shifts the coordinates of the rest
c  of the basis accordingly.  it also removes a null site from the
c  list.
c------------------------------------------------------------------
c  input/output:
c    iat,ndopx: parameters set in calling program
c    ibasis: number of sites in basis on input, null site removed
c            on output
c    iabs:   index of absorber in basis list on input, 1 on output
c    x,y,z:  (iat) on input: fractional coordinates of basis as listed
c            in input file, on output:  fractional coordinates
c            shifted to put absorber at point of symmetry
c    dopant*2: (iat,ndopx) list of atomic symbols, reordered on output
c            so that absorber is the first element of the array
c    tag*10: (iat) list of site tags, reordered on output so that
c            absorber is the first element of the array
c------------------------------------------------------------------
c      parameter (iat=50, ndopx=4)
      parameter(ibig=10)
c  this number needs to be bigger than ndopx, 10 should suffice

      dimension    x(iat),y(iat),z(iat)
      character*2  dopant(iat,ndopx), etoss(ibig), el, test
      character*10 tag(iat),ttoss

      test = 'ab'
      if (iabs.eq.0) goto 999

c  on input, first site is the point of symmetry and may be null
c  iabs may be any other point in the list.  want iabs listed first
c  and the null site at the end for easy removal.

c  get vector to shift absorber onto basis point for point transform
      xshift = x(iabs) - x(1)
      yshift = y(iabs) - y(1)
      zshift = z(iabs) - z(1)

c  swap first and absorber and last
      do 10 i=1,ndopx
        etoss(i) = dopant(1,i)
 10   continue
      ttoss = tag(1)
      xtoss = x(1)
      ytoss = y(1)
      ztoss = z(1)

      do 20 i=1,ndopx
        dopant(1,i) = dopant(iabs,i)
 20   continue
      tag(1)    = tag(iabs)
      x(1)      = x(iabs)
      y(1)      = y(iabs)
      z(1)      = z(iabs)

      do 30 i=1,ndopx
        dopant(iabs,i) = dopant(ibasis,i)
 30   continue
      tag(iabs)    = tag(ibasis)
      x(iabs)      = x(ibasis)
      y(iabs)      = y(ibasis)
      z(iabs)      = z(ibasis)

      do 40 i=1,ndopx
        dopant(ibasis,i) = etoss(i)
 40   continue
      tag(ibasis)    = ttoss
      x(ibasis)      = xtoss
      y(ibasis)      = ytoss
      z(ibasis)      = ztoss

c  remove null atom site from list
      el = dopant(ibasis,1)
      call case(test,el)
      if (el.eq.'nu') ibasis = ibasis-1

c  shift coordinates
      do 100 i=1,ibasis
        x(i) = x(i) - xshift
        y(i) = y(i) - yshift
        z(i) = z(i) - zshift
 100  continue

c  reset iabs to reflect the switch just made
      iabs = 1

 999  continue
      return
c  end subroutine basort
      end

      subroutine bperm(iat, nat, iperm, cell, ipt, st)

c  permute monoclinic and orthorhombic back to the non-standard setting

      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

      dimension st(iat,192,3), ipt(iat), cell(6)
      dimension abc(3), perm(3,3)

c       if (iperm.eq.22) then
c           print*,iperm
c           return
c       endif

      do 20 j=1,3
        do 10 i=1,3
          perm(i,j) = 0
 10     continue
        abc(j) = 0.e0
 20   continue

c  set elements fo permutation matrix
c             orthorhombic, cab
      if (iperm.eq.2) then
          perm(2,1) = 1.e0
          perm(3,2) = 1.e0
          perm(1,3) = 1.e0
c             orthorhombic, bca
      elseif (iperm.eq.3) then
          perm(3,1) = 1.e0
          perm(1,2) = 1.e0
          perm(2,3) = 1.e0
c             orthorhombic, a-cb
      elseif (iperm.eq.4) then
          perm(1,1) = 1.e0
          perm(3,2) = 1.e0
          perm(2,3)= -1.e0
c             orthorhombic, ba-c
      elseif (iperm.eq.5) then
          perm(2,1) = 1.e0
          perm(1,2) = 1.e0
          perm(3,3) = -1.e0
c             orthorhombic, -cba
      elseif (iperm.eq.6) then
          perm(3,1) = 1.e0
          perm(2,2) = 1.e0
          perm(1,3) = -1.e0
c             monoclinic, acb(2nd) to abc(1st)
      elseif (iperm.eq.12) then
          perm(1,1) = 1.e0
          perm(2,3) = 1.e0
          perm(3,2) = 1.e0
      endif

c     --- orthorhombic or monoclinic
      if (iperm.lt.20) then
          do 40 i=1,3
            do 30 j=1,3
              abc(i) = abc(i) + cell(j) * abs(perm(i,j))
 30         continue
 40       continue
          do 45 i=1,3
            cell(i) = abc(i)
 45       continue
      endif

c     --- tetragonal
c       else
c           abc(1) = cell(1) * sqrt(2.e0)
c           abc(2) = abc(1)
c           abc(3) = cell(3)

c       print*,'iperm=',iperm
c       print*,'before, after'
c       print*,cell(1),abc(1)
c       print*,cell(2),abc(2)
c       print*,cell(3),abc(3)
c       print*,'angles:',cell(4),cell(5),cell(6)


      do 70 i=1,nat
        do 60 j=1,ipt(i)
          xx = 0
          yy = 0
          zz = 0
c  permute fractional coordinates
          xx = st(i,j,1)*perm(1,1) + st(i,j,2)*perm(1,2) +
     $                st(i,j,3)*perm(1,3)
          yy = st(i,j,1)*perm(2,1) + st(i,j,2)*perm(2,2) +
     $                st(i,j,3)*perm(2,3)
          zz = st(i,j,1)*perm(3,1) + st(i,j,2)*perm(3,2) +
     $                st(i,j,3)*perm(3,3)
c           print*,i,st(i,j,1),xx
c           print*,i,st(i,j,2),yy
c           print*,i,st(i,j,3),zz
          st(i,j,1) = xx
          st(i,j,2) = yy
          st(i,j,3) = zz
 60     continue
 70   continue

      return
c  end subroutine bperm
      end
      subroutine equipt (isystm,isymce,cbravl,ng,tx,fx,spcgrp)
      implicit real (a-h,o-z)
c      implicit double precision (a-h,o-z)
c----------------------------------------------------------------------
c     interprets the  herman-mauguin symbol for space group
c     modified from a program by h.burzlaff, j.appl.cryst, v.15,
c     p.464 (1982)
c  double precision and character variable manipulations: br 8/92
c----------------------------------------------------------------------
c  input:
c    spcgrp*10:  hermann-maguin space group notation
c  output:
c    isystm:     integer denoting bravais lattice, (1..7)=(a,b,r,c,i,f,p)
c    isymce:     centrosymmetry flag, 1=centrosymmetric
c    cbravl*1:   character denoting bravais type, 1st letter in spcgrp
c    ng:         simple multiplicity of lattice type
c    tx(3,24):   fractional coordinates of simply multiple atoms
c    fx(3,3,24): rotational symmetries
c----------------------------------------------------------------------
      parameter (eps=.001e0)
      parameter (one = 1.e0, two = 2.e0, four = 4.e0)
      parameter (half = one/two, quart = one/four)
      integer      e(3,3), sys, isys(7), nul(3)
      dimension    sh(3), te(3)
      dimension    ss(24,3,3), tx(3,24), ts(24,3), fx(3,3,24)
      character*10 spcgrp

      character*1  cspace(10), cbra(7), cbravl, csym(3,4), cbr
      character*1  cm, cn
      character*2  test
      logical      blank

      data (sh(i),  i=1,3)  /0.25e0, 0.25e0, 0.25e0/
      data (nul(i), i=1,3)  /0, 0, 0/
      data (isys(i),i=1,7)  /7, 1, 2, 4, 6, 5, 3/
      data (cbra(i),i=1,7)  /'p','a','b','c','f','i','r'/

c%%%%%
      test  = 'ab'
c --- this stuff was added after running the code through g77 --------
      blank = .false.
      ichar = 0
      nbr   = 0
c --------------------------------------------------------------------
c  the value of test *must* be of the same case as the values of cbra
c-----------------------------------------------------
c  break spcgrp up into its 10 components
c  cspace contains the (up to) 10 characters in the space group designation
      do 9001 i=1,10
        cspace(i) = ' '
9001  continue
      do 9011 i=1,10
        read (spcgrp(i:i),4000) cspace(i)
        call case(test,cspace(i))
9011  continue
4000  format(a1)

c-----------------------------------------------------
c     initialize some things
      nsym = -1
      ng   =  1
      ns   =  1
      do 24 i=1,24
         do 22 j=1,3
           ts(i,j) = 0.e0
           tx(j,i) = 0.e0
           do 20 k=1,3
              e(j,k)    = 0
              ss(i,j,k) = 0.e0
              fx(j,k,i) = 0.e0
 20        continue
 22      continue
 24   continue
      ss(1,1,1) = 1.e0
      ss(1,2,2) = 1.e0
      ss(1,3,3) = 1.e0
      do 32 i=1,3
         te(i) = 0.e0
         do 30 j=1,4
            csym(i,j) = ' '
 30      continue
 32   continue

c-----------------------------------------------------
c  begin interpreting the space group

      do 80 i=1,10
c        reset blank (flag) and ichar each time a blank is found
c        then go on to the next character
         if (cspace(i).eq.' ') then
             blank  = .true.
             ichar  = 0
             goto 80
         endif

c  n.b. this could be cleaned up since spcgrp is trimled before this

c  if bravais lattice type has been found...
         if (nsym.ge.0) goto 70

c  look for bravais lattice type    1..7 = p,a,b,c,f,i,r
c  cbr for internal use, i/cbravl for external use
         do 50 j=1,7
            cbr = cbra(j)
            nbr = j
            if (cspace(i).eq.cbra(j)) goto 60
50       continue

60       continue
         nsym    = 0
         cbravl  = cspace(i)
         goto 80

c  nsym counts the three symbols in the space group after the lattice type.
c  ichar counts the characters in that symbol.  nsym=[1..3],ichar=[1..4].
c  blank=t flags the beginning of a new symbol.
70       continue
         if (blank) nsym  = nsym+1
         ichar            = ichar+1
         csym(nsym,ichar) = cspace(i)
         blank            = .false.
         if (cspace(i).eq.'/') ns = 0
 80   continue

c      print*,'spcgrp=',spcgrp
c      print*,'cspace=',cspace

c-----------------------------------------------------
c  determine the bravais system
c  sys=1...6 ==> (tric,mono,ortho,tetra,hex,cubic)

c  cubes have triad axes.
      do 90 i=1,4
         if (csym(2,i).ne.'3') goto 90
         sys = 6
         goto 145
 90   continue

c  hexagonal cells have hexads, trigonal cells have triads
c  tetragonal cells have tetrads.
      do 110 i=1,4
         if ( (csym(1,i).eq.'3').or.(csym(1,i).eq.'6') ) then
               sys = 5
               goto 145
         endif
         if (csym(1,i).eq.'4') then
             sys = 4
             goto 145
         endif
110   continue

c     choose between triclinic and monoclinic
      if (nsym.eq.1) then
c         monoclinic might only have one symbol, but it's not "1"
          if ( (csym(1,1).eq.'1').or.(csym(1,1).eq.'-') ) then
                sys = 1
                goto 145
          endif
          sys = 2
c         fix csym for point groups 2 and m, but not 2/m
          do 130 i=1,4
             csym(2,i) = csym(1,i)
130       continue
          csym(1,1) = '1'
          csym(3,1) = '1'
      endif

c%#@%#@  this seems to be wrong for group 2/m
c  orthorhombic by default unless monoclinic pt. grp. 2 or m
      sys=3
      if ( (csym(1,1).eq.'1').or.(csym(2,1).eq.'1') ) sys = 2

c---------------------------------------------------------------------
c   leap ahead to the appropriate section of code for the bravais
c   type in question
c   jump to these line numbers: (1000,1100,1200,1300,1400,1500)
c   for appropriate sys, sys=[1..6] according to lattice type:
c   (tric,mono,ortho,tetra,hex,cubic)
145   continue

      if     (sys.eq.1) then
              goto 1000
      elseif (sys.eq.2) then
              goto 1100
      elseif (sys.eq.3) then
              goto 1200
      elseif (sys.eq.4) then
              goto 1300
      elseif (sys.eq.5) then
              goto 1400
      elseif (sys.eq.6) then
              goto 1500
      endif

c---------------------------------------------------------------------
c  triclinic,  described by group containing only the identity (p1)
c              or the identity and a parity operation (p-1)
1000  if (csym(1,1).eq.'-') ns = 0
      goto 760

c---------------------------------------------------------------------
c  monoclinic, primitive (p) or one-face-centered (c).
c              point gps: diad (p2,p21,c2),mirror (pm,cm),glide (pc,cc),
c                         combined (p2/m,p21/m,p2/c,p21/c,c2/m,c2/c)
1100  ng = 2
      ind = 0
      do 180 i=1,3
        if (csym(i,1).ne.'1') ind = i
180   continue
      id = 1
      if (csym(ind,1).eq.'2') id = -1
      do 190 i=1,3
        ss(2,i,i) =  ss(1,i,i)*id
190   continue
      ss(2,ind,ind) = -ss(2,ind,ind)
      do 220 i=1,3
        if ( (csym(i,1).eq.'2').and.(csym(i,2).eq.'1') ) ts(2,i)=0.5
        do 210 j=1,4
          if (csym(i,j).eq.'a') ts(2,1) = 0.5
          if (csym(i,j).eq.'b') ts(2,2) = 0.5
          if (csym(i,j).eq.'c') ts(2,3) = 0.5
          if (csym(i,j).eq.'n') goto 200
          goto 210
 200      continue
          k=i+1
          if (k.gt.3) k = k-3
          ts(2,k)     = 0.5
          ts(2,6-k-i) = 0.5
 210    continue
 220  continue
      goto 760

c---------------------------------------------------------------------
c --- orthorhombic
1200  continue
      ng  = 4
      ic  = 0
      ind = 1

c  if not point group 222...
      if ( (csym(1,1).ne.'2').and.(csym(2,1).ne.'2').and.
     $     (csym(3,1).ne.'2') )  then
            ind = -1
            ns = 0
      endif

      do 240 i=1,3
         id = 1
c        if any diads exist...
         if (csym(i,1).eq.'2') id = -1
         do 230 j=1,3
           ss(1+i,j,j) = ss(1,j,j)*id*ind
230      continue
         ss(1+i,i,i) = -ss(1+i,i,i)
240   continue

      do 282 i=1,3
c        if screw diads exist...
         if ( (csym(i,1).eq.'2').and.(csym(i,2).eq.'1') )
     $          ts(1+i,i) = 0.5
         do 280 j=1,4
c           if a glide plane...; if b glide plane...; if c glide plane...;
c           if diagonal glide plane or 1/4 glide plane
            if  (csym(i,j).eq.'a') ts(1+i,1) = 0.5
            if  (csym(i,j).eq.'b') ts(1+i,2) = 0.5
            if  (csym(i,j).eq.'c') ts(1+i,3) = 0.5
            if ((csym(i,j).eq.'n').or.(csym(i,j).eq.'d')) goto 250
            goto 280

250         continue
            k = i+1
            if (k.gt.3) k = k-3
            if (csym(i,j).eq.'d') goto 260

c           case of not diagonal glide plane
            ts(1+i,k)     = 0.5
            ts(1+i,6-k-i) = 0.5
            goto 280

c           case of 1/4 glide plane
260         ic = 1
c           if no centrosymmetry...
            if (ns.eq.1) goto 270
            ts(1+i,k)     = 0.25
            ts(1+i,6-k-i) = 0.25
            goto 280
 270        ts(1+i,1)     = 0.25
            ts(1+i,2)     = 0.25
            ts(1+i,3)     = 0.25
 280     continue
 282  continue

c
      if (ic.eq.1) goto 760
c     if no centrocymmetry...
      if (ns.eq.1) goto 360

      do 290 i=1,3
         k = 1+i
         if (k.gt.3) k = k-3
         tc = ts(1+k,i)+ts(1+6-i-k,i)
         if (tc.eq.1.0) tc = 0.0
         ts(1+i,i) = tc
290   continue

c----------------------------------------------------------------
c  special treatment of c m m a, c m c a, i m m a
c  ma counts the number of "m's" (ie 2, 1, or 2)
      if ((cbr.eq.'p').or.(cbr.eq.'f')) goto 760
      ma = 0
      do 310 i=1,3
         do 300 j=1,4
           if (csym(i,j).eq.'m') nul(i)=1
300      continue
         ma = ma+nul(i)
310   continue
c     if i m m a...do origin shift
      if ( (cbr.eq.'i').and.(ma.eq.2) ) goto 330
c     weed out remaining orthorhombic c and i space groups
      if ( (ma.eq.0).or.(ma.eq.3).or.(cbr.eq.'i') ) goto 760
c%#%#%#%#%#%#%#%#
      do 320 i=1,3
c  is this a typo?!   what is the purpose of the do loop -- this is equally
c                     mysterious in sexie's equipt.  will deal with this
c                     when time comes.
         if (nul(nbr-1).eq.1) goto 760
         sh(nbr-1) = 0
320   continue

c --- origin shift
 330  do 352 i=1,ng
         do 350 j=1,3
            do 340 k=1,3
               id = 1
               if (j.ne.k) id = 0
               ts(i,j) = ts(i,j) + (id-ss(i,j,k))*sh(k)
340         continue
            if (ts(i,j).gt.1.)  ts(i,j) = ts(i,j)-1
            if (ts(i,j).lt.1.0) ts(i,j) = ts(i,j)+1
350      continue
352   continue
      goto 760

c...where to come for no centrocymmetry
360   ic = 0
      do 370 i=1,3
c                                                  .eq.-1
        if ( abs((ss(1+i,1,1)*ss(1+i,2,2)*ss(1+i,3,3)) +1).lt.eps) ic=1
370   continue
      if (ic.eq.1) goto 410
      tc = ts(2,1)+ts(3,2)+ts(4,3)
      if (abs(tc).lt.eps) goto 760
      do 400 i=1,3
         k = i+1
         if (k.gt.3) k = k-3
         if (tc.gt.0.5) goto 380
         if (abs(ts(1+i,i)).lt.eps) goto 400
         m = i-1
         if (m.eq.0) m = m+3
         ts(1+m,i) = 0.5
         goto 400
 380     if (tc.gt.1.0) goto 390
         if (abs(ts(1+i,i)).gt.eps) goto 400
         l = k+1
         if (l.gt.3) l = l-3
         ts(1+k,l) = 0.5
         ts(1+l,k) = 0.5
         goto 400
 390     ts(1+i,k) = 0.5
 400  continue
      goto 760

410   do 420 i=1,3
        if (abs(ss(1+i,1,1)*ss(1+i,2,2)*ss(1+i,3,3)-1).lt.eps) id=i
420   continue
      do 450 i=1,3
         tc = ts(2,i)+ts(3,i)+ts(4,i)
         if ( (abs(tc).lt.eps).or.(abs(tc-1.0).lt.eps) ) goto 450
c        if mirror and diagonal glide planes exist...
         if ( ((csym(1,1).eq.'m').and.(csym(2,1).eq.'n')) .or.
     $        ((csym(2,1).eq.'m').and.(csym(3,1).eq.'n')) .or.
     $        ((csym(3,1).eq.'m').and.(csym(1,1).eq.'n')) ) goto 440
         do 430 j=1,3
            if (id.eq.j) goto 430
            if (abs(ts(1+j,i)-0.5).lt.eps) goto 430
            ts(1+j,i) = 0.5
 430     continue
         goto 450
 440     k=i-1
         if (k.eq.0) k = k+3
         ts(1+k,i) = 0.5
 450  continue
      goto 760

c---------------------------------------------------------------------
c --- tetragonal


1300  ng = 4
      if (nsym.eq.3) ng=8
      ss(2,1,2) = -1
      ss(2,2,1) =  1
      ss(2,3,3) =  1
      cm        = csym(1,1)
      cn        = csym(1,2)

      do 472 i=1,3
         do 470 j=1,3
c           if enantiomorphous...
            if (cm.eq.'-') ss(2,i,j) = -ss(2,i,j)
 470     continue
 472  continue

      if (cm.eq.'-') goto 500
c     if screw axes of 1/4, 1/2. or 3/4 rotation exist...
      if (cn.eq.'1') ts(2,3) = 0.25
      if (cn.eq.'2') ts(2,3) = 0.5
      if (cn.eq.'3') ts(2,3) = 0.75

c   if horizontal diagonal glide plane or hdgp & # of symbols=3...
      if ( (csym(1,3).eq.'n') .or. ((csym(1,4).eq.'n').and.(nsym.eq.3)))
     $                 ts(2,1) = 0.5

c   if hdgp & # of symbols=1 (p42/n) or 1/4 screw axis & no cent. & i...
      if ( ((csym(1,4).eq.'n').and.(nsym.eq.1)) .or.
     $     ((cn.eq.'1').and.(ns.eq.1).and.(cbr.eq.'i')) ) ts(2,2) = 0.5

c   if 1/4 screw axis & no cent. & i...
      if (  (cn.eq.'1').and.(ns.eq.0).and.(cbr.eq.'i') ) goto 480

c   if secondary 1/4 screw axis or no hdgp & vert. diag. glide plane &
c                                  third mirror plane...
      if ( (csym(2,2).eq.'1') .or.
     $    ((csym(1,4).ne.'n').and.(csym(2,1).eq.'n').and.
     $     (csym(3,1).eq.'m')) ) goto 490

      goto 500

480   continue
      ts(2,1) = 0.25
      ts(2,2) = 0.75
c     if point group 4, -4, or 4/m...
      if (nsym.eq.1) ts(2,1) = 0.75
      if (nsym.eq.1) ts(2,2) = 0.25
      goto 500

490   continue
      ts(2,1)   = 0.5
      ts(2,2)   = 0.5

500   continue
      ss(3,1,1) = -1
      ss(3,2,2) = -1
      ss(3,3,3) =  1
      ts(3,1)   = ss(2,1,2)*ts(2,2)+ts(2,1)
      ts(3,2)   = ss(2,2,1)*ts(2,1)+ts(2,2)
      ts(3,3)   = ss(2,3,3)*ts(2,3)+ts(2,3)

      do 510 i=1,3
        if ( (cbr.eq.'i').and.(abs(ts(3,1)-0.5).lt.eps).and.
     $       (abs(ts(3,2)-0.5).lt.eps).and.(abs(ts(3,3)-0.5).lt.eps) )
     $                   ts(3,i) = 0.0
510   continue

      do 524 i=1,3
         ts(4,i) = ts(2,i)
         do 522 j=1,3
            ts(4,i) = ts(4,i)+ss(2,i,j)*ts(3,j)
            do 520 k=1,3
               ss(4,i,j) = ss(4,i,j)+ss(2,i,k)*ss(3,k,j)
 520        continue
 522     continue
 524  continue

c     if point group 4, -4, or 4/m then done...
      if (nsym.eq.1) goto 760
c     if centrosym...
      if (ns.eq.0) goto 560

      cm = csym(2,1)
      cn = csym(3,1)
c     if no secondary diad exists...
      if ( (cm.ne.'2').and.(cn.ne.'2') ) goto 550
c     if 2 secondary diads exist...
      if ( (cm.eq.'2').and.(cn.eq.'2') ) goto 540
c     if 2nd symbol 1st char = 2 skip to 530..
      if (cm.eq.'2') goto 530

c     if c glide plane or diag. glide plane exists...
      if ( (cm.eq.'c').or.(cm.eq.'n') )  te(3) = 0.5
      e(1,1) = -1
      e(2,2) =  1
      e(3,3) =  1
c     skip to 570 if no diag or b glide plaes exist...
      if ( (cm.ne.'n').and.(cm.ne.'b') ) goto 570
      te(1)  = 0.5
      te(2)  = 0.5
      goto 570

530   continue
      e(1,1) =  1
      e(2,2) = -1
      e(3,3) = -1
c     if c or 1/4 glide plane...
      if (cn.eq.'c') te(3) = 0.5
      if (cn.eq.'d') te(3) = 0.25
      if (cn.eq.'d') te(2) = 0.5
c     if secondary diad is not screw...
      if (csym(2,2).ne.'1') goto 570
      te(1)  = 0.5
      te(2)  = 0.5
      goto 570

540   continue
      e(1,2) =  1
      e(2,1) =  1
      e(3,3) = -1
c     if no secondary screw or i...
      if ( (csym(2,2).ne.' ').or.(cbr.eq.'i').or.(csym(1,2).eq.' ') )
     $     goto 570
c     if principle screw axis is 1/4, 1/2, 3/4...
      if (csym(1,2).eq.'1') te(3) = 0.75
      if (csym(1,2).eq.'2') te(3) = 0.5
      if (csym(1,2).eq.'3') te(3) = 0.25
      goto 570

550   continue
      cm     = csym(2,1)
      e(1,1) = -1
      e(2,2) =  1
      e(3,3) =  1
c     if parallel diag(n), b, or c mirror planes...
      if ( (cm.eq.'c').or.(cm.eq.'n') ) te(3) = 0.5
      if ( (cm.eq.'n').or.(cm.eq.'b') ) te(1) = 0.5
      if ( (cm.eq.'n').or.(cm.eq.'b') ) te(2) = 0.5
      goto 570

560   continue
      e(1,1) = -1
      e(2,2) =  1
      e(3,3) =  1
c     if parallel glide plane...
      if ( (csym(2,1).eq.'c').or.(csym(2,1).eq.'n') ) te(3) = 0.5
      if ( (csym(2,1).eq.'b').or.(csym(2,1).eq.'n') ) te(2) = 0.5
      if ( (csym(2,1).eq.'b').or.(csym(2,1).eq.'n') ) te(1) = 0.5
c     if parallel glide plane...
      if ( (csym(1,3).eq.'n').or.(csym(1,4).eq.'n') ) te(1)= te(1)+0.5
 570  ne = 4
      goto 740

c---------------------------------------------------------------------
c --- hexagonal
1400  continue
c      print*,'starting hex'
      ng = 3
      cm  = csym(1,1)
      cn  = csym(1,2)

c     if centrosymmetric...
      if ( (cm.eq.'-').and.(cn.eq.'3') ) ns = 0
c     if a hexad exists...
      if (cm.eq.'6') goto 610
      ss(2,1,2) = -1
      ss(2,2,1) =  1
      ss(2,2,2) = -1
      ss(2,3,3) =  1

c     if a 1/3, 2/3 screw triad exists...
      if (cn.eq.'1') ts(2,3) = 1.0/3.0
      if (cn.eq.'2') ts(2,3) = 2.0/3.0
      ss(3,1,1) = -1
      ss(3,2,1) = -1
      ss(3,1,2) =  1
      ss(3,3,3) =  1
      ts(3,3)   =  2*ts(2,3)
      if (ts(3,3).ge.1.0) ts(3,3) = ts(3,3)-1

c     if point group -6...
      if ( (nsym.eq.1).and.(cn.ne.'6') ) goto 760
c     if not point group -6m2...
      if (cn.ne.'6') goto 600
      ng = ng+ng
      do 594 i=1,3
         do 592 j=1,3
            do 590 k=1,3
               ss(3+i,j,k) = ss(i,j,k)
               ss(3+i,3,3) = -1
 590        continue
 592     continue
 594  continue

600   continue
      if (nsym.eq.1) goto 760
      if ( (csym(2,1).ne.'c').and.(csym(3,1).ne.'c') ) goto 630
      ts(4,3)   = 0.5
      ts(5,3)   = 0.5
      ts(6,3)   = 0.5
      goto 630

610   continue
      ng        = ng+ng
      ss(2,1,1) =  1
      ss(2,1,2) = -1
      ss(2,2,1) =  1
      ss(2,3,3) =  1
c     if 1/6..5/6 screw hexad...
      if (cn.eq.'1') ts(2,3) = 1.0/6.0
      if (cn.eq.'2') ts(2,3) = 2.0/6.0
      if (cn.eq.'3') ts(2,3) = 3.0/6.0
      if (cn.eq.'4') ts(2,3) = 4.0/6.0
      if (cn.eq.'5') ts(2,3) = 5.0/6.0

      do 626 i=1,4
         do 624 j=1,3
            ts(2+i,j)=ts(2,j)
            do 622 k=1,3
               ts(2+i,j) = ts(2+i,j)+ss(2,j,k)*ts(1+i,k)
               if (ts(2+i,j).gt.1.0) ts(2+i,j) = ts(2+i,j)-1.0
               do 620 l=1,3
                  ss(2+i,j,k) = ss(2+i,j,k)+ss(2,j,l)*ss(1+i,l,k)
 620           continue
 622        continue
 624     continue
 626  continue
      if (nsym.eq.1) goto 760

630   continue
      ng  = ng+ng
      cm  = csym(2,1)
      cn  = csym(3,1)
c     if x,y directions are identity axes...
      if (cm.eq.'1') goto 650
c     if secondary diad exists...
      if (cm.eq.'2') goto 640
      e(1,2) = -1
      e(2,1) = -1
      e(3,3) =  1
c     if c glide plane exists...
      if (cm.eq.'c') te(3)=0.5
      goto 670

640   continue
      e(1,2) =  1
      e(2,1) =  1
      e(3,3) = -1
      te(3)  = 2.0*ts(2,3)

c   group p 31 1 2 and p 32 1 2
      if ( (csym(1,1).eq.'3') .and.
     $    ((csym(1,2).eq.'1').or.(csym(1,2).eq.'2')) ) te(3) = 0
      goto 670

650   continue
c     if secondary diad exists...
      if (cn.eq.'2') goto 660
      e(1,2) = 1
      e(2,1) = 1
      e(3,3) = 1
c     if c glide plane exists...
      if (cn.eq.'c') te(3) = 0.5
      goto 670

660   continue
      e(1,2) = -1
      e(2,1) = -1
      e(3,3) = -1
      te(3)  = 2.0*ts(2,3)
      if (te(3).gt.1.0) te(3) = te(3)-1.0

670   continue
      ne = 6
c     if trigonal (3**) or (-3**)...
      if ( (csym(1,1).eq.'3') .or.
     $    ((csym(1,2).eq.'3').and.(csym(1,1).eq.'-')) ) ne = 3
      goto 740

c---------------------------------------------------------------------
c   cubic
1500  ng = 12
      if (nsym.eq.3) ng = 24
c     if no primary diad or tetrad exists...
      if ( (csym(1,1).ne.'2').and.(csym(1,1).ne.'4').and.
     $     (csym(1,1).ne.'-') ) ns = 0

      do 692 i=1,3
         do 690 j=1,3
            ss(1+i,j,j) =  1
            if (i.eq.j) goto 690
            ss(1+i,j,j) = -1
            if (csym(1,1).eq.'n') ts(1+i,j) = half
            if (csym(1,1).eq.'d') ts(1+i,j) = quart
 690      continue
 692   continue

c     if (no a and no 1/4 glide plane and no 1/4 and no 3/4 screw) or
c     if (face centered)...
      if ( ((csym(1,1).ne.'a').and.(csym(3,1).ne.'d').and.
     $      (csym(1,2).ne.'3').and.(csym(1,2).ne.'1')) .or.
     $      (cbr.eq.'f') )   goto 710

      do 700 i=1,3
         ts(1+i,i) = half
         k = i+1
         if (k.eq.4) k = 1
         ts(1+i,k) = half
700   continue

710   continue
      do 724 i=1,4
         do 722 j=1,3
            do 720 k=1,3
               l = j+1
               if (l.eq.4) l = 1
               m = j-1
               if (m.eq.0) m = 3
               ss(4+i,j,k) = ss(i,l,k)
               ss(8+i,j,k) = ss(i,m,k)
               ts(4+i,j)   = ts(i,l)
               ts(8+i,j)   = ts(i,m)
720         continue
722      continue
724   continue

      if (ng.eq.12) goto 760
      ne     = 12
      e(1,2) = 1
      e(2,1) = 1
      e(3,3) = 1
c     if secondary diad exists...
      if (csym(3,1).eq.'2') e(3,3) = -1
c     if c glide plane exists...
      if (csym(3,1).eq.'c') te(3)  = half

c        if diag glide plane or 2/4 screw diad..
c        if 1/4 glide or 1/4 screw or 3/4 screw..
      do 730 i=1,3
         if ( (csym(3,1).eq.'n').or.(csym(1,2).eq.'2') ) te(i) = half
         if ( (csym(3,1).eq.'d').or.(csym(1,2).eq.'1').or.
     $        (csym(1,2).eq.'3') )                       te(i) = quart
730   continue

c     if 1/4 screw or p...
      if (  (csym(1,2).eq.'1').and.(cbr.eq.'p') )      te(1) = quart * 3
c     if (not 1/4 screw or not i) and (not 3/4 screw or not p)...
      if ( ((csym(1,2).ne.'1').or.(cbr.ne.'i')) .and.
     $     ((csym(1,2).ne.'3').or.(cbr.ne.'p')) )        goto 740

      te(2) = quart * 3
      te(3) = quart * 3

c---------------------------------------------------------------------
 740  do 756 i=1,ne
         do 754 j=1,3
            ts(ne+i,j) = te(j)
            do 752 k=1,3
               ts(ne+i,j) = ts(ne+i,j)+e(j,k)*ts(i,k)
               do 750 l=1,3
                  ss(ne+i,j,k) = ss(ne+i,j,k)+e(j,l)*ss(i,l,k)
750            continue
752         continue
754      continue
756   continue

c---------------------------------------------------------------------
c  just about done.  put all atoms in central cell, finish up fx, mark
c  system and centrosymmetry, and make sure ng>=1.
760   continue

c++++++++++++++++
c  put located atoms back into central cell
      do 812 i=1,ng
         do 810 j=1,3
 770        if (ts(i,j).ge.1.0) then
                ts(i,j) = ts(i,j)-1
                goto 770
            endif
 790        if (ts(i,j).lt.0.0) then
                ts(i,j) = ts(i,j)+1
                goto 790
            endif
 810     continue
 812  continue

c================
c---------------------------------------------------------------------
c  tell the rest of the program what system the lattice is and whether
c  is has centrosymmetry
      isystm = isys(sys)
      if (ns.eq.0) isymce = 1
      if (ns.eq.1) isymce = 0

c---------------------------------------------------------------------
c --- swap ss(i,j,k) to fx(k,i,j)
      do 840 k=1,ng
         do 830 j=1,3
            do 820 i=1,3
              fx(i,j,k) = ss(k,j,i)
820         continue
            tx(j,k)   = ts(k,j)
830      continue
840   continue

c---------------------------------------------------------------------
c  there is always at least one element in a group
c  (identity transformation)
      if (ng.eq.0) ng = 1

      return
c end subroutine equipt
      end

      subroutine fperm(iat, nat, iperm, cell, x, y, z)

      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

      dimension cell(6), x(iat), y(iat), z(iat)
      dimension abc(3), perm(3,3)

      do 20 j=1,3
        do 10 i=1,3
          perm(i,j) = 0
 10     continue
        abc(j) = 0.e0
 20   continue

c  set elements fo permutation matrix
c             orthorhombic, cab
      if (iperm.eq.2) then
          perm(1,2) = 1.e0
          perm(2,3) = 1.e0
          perm(3,1) = 1.e0
c             orthorhombic, bca
      elseif (iperm.eq.3) then
          perm(1,3) = 1.e0
          perm(2,1) = 1.e0
          perm(3,2) = 1.e0
c             orthorhombic, a-cb
      elseif (iperm.eq.4) then
          perm(1,1) = 1.e0
          perm(2,3) = 1.e0
          perm(3,2)= -1.e0
c             orthorhombic, ba-c
      elseif (iperm.eq.5) then
          perm(1,2) = 1.e0
          perm(2,1) = 1.e0
          perm(3,3) = -1.e0
c             orthorhombic, -cba
      elseif (iperm.eq.6) then
          perm(1,3) = 1.e0
          perm(2,2) = 1.e0
          perm(3,1) = -1.e0
c             monoclinic, acb(2nd) to abc(1st)
      elseif (iperm.eq.11) then
          perm(1,1) = 1.e0
          perm(2,3) = 1.e0
          perm(3,2) = 1.e0
c             tetragonal, rotated --> standard
      elseif (iperm.eq.22) then
          perm(1,1) = 1.e0
          perm(1,2) = 1.e0
          perm(2,1) = -1.e0
          perm(2,2) = 1.e0
          perm(3,3) = 1.e0
      endif

c     --- orthorhombic or monoclinic
      if (iperm.lt.20) then
          do 40 i=1,3
            do 30 j=1,3
              abc(i) = abc(i) + cell(j) * abs(perm(i,j))
 30         continue
 40       continue
c     --- tetragonal
      else
          abc(1) = cell(1) / sqrt(2.e0)
          abc(2) = abc(1)
          abc(3) = cell(3)
      endif
c       print*,'iperm=',iperm
c       print*,'before, after'
c       print*,cell(1),abc(1)
c       print*,cell(2),abc(2)
c       print*,cell(3),abc(3)
c       print*,'angles:',cell(4),cell(5),cell(6)

      do 45 i=1,3
        cell(i) = abc(i)
 45   continue

      do 50 i=1,nat
        xx = x(i)*perm(1,1) + y(i)*perm(1,2) + z(i)*perm(1,3)
        yy = x(i)*perm(2,1) + y(i)*perm(2,2) + z(i)*perm(2,3)
        zz = x(i)*perm(3,1) + y(i)*perm(3,2) + z(i)*perm(3,3)
c         print*,i,x(i),xx
c         print*,i,y(i),yy
c         print*,i,z(i),zz
        x(i) = xx
        y(i) = yy
        z(i) = zz
 50   continue

      return
c  end subroutine fperm
      end
      subroutine genmul(ng,isymce,ibravl,igen)
c--------------------------------------------------------------
c  copyright 1993 university of washington         bruce ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c-------------------------------------------------------------------
c  calculates general multiplicity (igen) of unit cell, i.e. the
c  number of times all of the symmetries of a cell generate a new
c  point from any arbitrary point
c-------------------------------------------------------------------
c  input:
c    ng     = multiplicity before bravais translations and
c             centrosymmetry operation
c    isymce = centrosymmetry flag
c    ibravl = bravais lattice type, 1..7 = p i r f a b c
c  output:
c    igen   = general multiplicity, product of simple multiplicity,
c             mult. of bravais lattice, and mult. of centrosymmetry
c-------------------------------------------------------------------
c    icmf   = mult. factor due to centrosymmety
c    ibrmf  = mult. factor due to bravais lattices
c    igen   = total multiplicity
c-------------------------------------------------------------------

      icmf  = 1
      igen  = ng
      ibrmf = 1
c                                          if centrosymmetric...
      if (isymce.eq.1) then
          icmf = 2
          igen  = ng*icmf
      endif
c                                          2=i, body center
c                                          3=r, rhombohedral
c                                          4=f, face center
      if ((ibravl.ge.2).and.(ibravl.le.4)) then
          ibrmf = ibravl
          igen  = ng*icmf*ibrmf
c                                          5..7=a,b,c, single face center
      elseif ((ibravl.ge.5).and.(ibravl.le.7)) then
          ibrmf = 2
          igen   = ng*icmf*ibrmf
      end