      SUBROUTINE MULTF(elemt,kodlp,npos,bpos,xyzpos)
C-----CALCULATES MULTIPLICITIES OF SPECIAL POSITIONS                    MULT0001
C-----and outputs the coordiates of atoms in the unit cell--------
C to unit=unitstruct, the file *.struct

      implicit none      
      include 'arraysize.inc'
      include 'io_files.inc'
      
c...	common blocks
      character*76 COMPND
      character*4 ielmt(16)
      real FOCCU(50,16), FMULT(50,16)
      real FNEUT(16)
      real A,B,C,ALPHA,BETA,GAMMA
      real SYMSY				!! Unused in code
      integer ISYMCE, NCOMPO
      integer NA(16)
      real wl, SL,SH, AGUIN, BGUIN
      real BTEMP(50,16), DELFR(16), DELFI(16)
      integer IMAGE, ns
      character*4 norm,iano
      COMMON /INPUT/ COMPND,IELMT,FOCCU,FMULT,FNEUT,
     &		A,B,C,ALPHA,BETA,GAMMA,
     &		SYMSY,ISYMCE,NCOMPO,NA,WL,SL,SH,NORM,IMAGE,AGUIN,
     &		BGUIN,BTEMP,DELFR,DELFI,IANO,NS

c01Dec      COMMON /INPUT/ COMPND,IELMT,FOCCU,FMULT,FNEUT,A,B,C,ALPHA,BE
c01Dec     1TA,GAMMA,SYMSY,ISYMCE,NCOMPO,NA,WL,SL,SH,NORM,IMAGE,AGUIN,BG
c01Dec     2UIN,BTEMP,DELFR,DELFI,IANO,NS

      integer ibravl, isystm
      common /cell/ ibravl,isystm

      character*2 symlp
      character*1 symbr
      character*4 symwl
      COMMON /INPUTCHAR/ SYMLP,SYMBR,symwl

      double precision x(50,16),y(50,16),z(50,16)
      common /inputdoub/ x,y,z

      double precision ts(3,24),fs(3,3,24)
      COMMON /EQUO/ TS,FS

      integer ICR,IPR,ILU
      COMMON /UNITS/ ICR,IPR,ILU

      real FNEU
      common /neu/ fneu

c...	passed variables
      character*4 elemt(16)
      integer kodlp
      integer npos
      real bpos(200)
      double precision xyzpos(natoms_max,3)
c01Dec      double precision xyzpos(nmax,3)
      
c...	local variables

c01Dec      integer nmax
      double precision xyzunique(96,3)
      double precision xyz(3)
      double precision dif,dif1, anti
      double precision tsb(3)
      double precision b1(3,8)
      integer nbra,ia,ic
C      real b1(3,8)    ! this was (3,9) above ?
C      real xyzunique(96,3)
      integer nunique, natom
      logical uniflag
      integer i, j, nb
      
C ibravl is the Bravais lattice number from function BRA(SYMBR)
C which indexes (blank) P I R F A B C       
C
C nbra is the number of Bravais points
      
      DATA B1(1,1)/.0d0/,B1(2,1)/.0d0/,B1(3,1)/.0d0/,B1(1,2)/.5d0/,B1(2,
     12)/.5d0/,B1(3,2)/.5d0/,B1(1,3)/2.0d0/,B1(2,3)/1.0d0/,B1(3,3)/1.0d0
     2/,B1(1,4)/.5d0/,B1(2,4)/.5d0/,B1(3,4)/.0d0/,B1(1,5)/.0d0/,B1(2,5)/
     3.5d0/,B1(3,5)/.5d0/,B1(1,6)/.5d0/,B1(2,6)/.0d0/,B1(3,6)/.5d0/,B1(1
     4,7)/.5d0/,B1(2,7)/.5d0/,B1(3,7)/.0d0/,B1(1,8)/1.0d0/,B1(2,8)/2.0d0
     5/,B1(3,8)/2.0d0/

      do i=1,3
         b1(i,3)=b1(i,3)/3.0d0
         b1(i,8)=b1(i,8)/3.0d0
      enddo

C B1 as a matrix:  0   1/2   2/3   1/2    0    1/2   1/2   1/3
C                  0   1/2   1/3   1/2   1/2    0    1/2   2/3
C                  0   1/2   1/3    0    1/2   1/2    0    2/3

C start with no atoms
      npos=0
         
      NBRA=1
      IF (IBRAVL.GT.1) NBRA=2       ! not a primitive lattice
      IF (IBRAVL.EQ.3) NBRA=3       ! rhombohedral (oh dear!)
      IF (IBRAVL.EQ.4) NBRA=4       ! face-centred

      write (unitstruct, *) '>> multf():'

      IF (ISYMCE.NE.1) then
c        write(7,25) nbra,
c     &     ' Bravais lattice positions and non-centrosymmetric.'
        write(unitstruct,25) nbra,
     &     ' Bravais lattice positions and non-centrosymmetric.'
      else
c        write(7,25) nbra,
c     &     ' Bravais lattice positions and centrosymmetric.'
        write(unitstruct, 25) nbra,
     &     ' Bravais lattice positions and centrosymmetric.'
      endif  
 25   format(i1,a)

C      write(7,24) 'Atom','  fi  ','  X   ','  Y   ','  Z   '
C      write(7,24) '----',' ---- ','----- ','----- ','----- '
C 24   format(a4,3x,a6,4x,a6,1x,a6,1x,a6)
      
cUnused      nopos=0    ! number of positions which are occupied in the unit cell

C-----LOOP OVER ATOMKINDS                                               MULT0020
      DO 7 IC=1,NCOMPO
        NATOM=NA(IC)
        CALL ANOMA (SYMWL,ELEMT(IC),KODLP)

C-----LOOP OVER ATOMS OF ONE KIND                                       MULT0023
        DO 6 IA=1,NATOM
          nunique=1
C  put in xyz coordinates of the atom, these are now double precision
          XYZUNIQUE(1,1)=X(IA,IC)
          XYZUNIQUE(1,2)=Y(IA,IC)
          XYZUNIQUE(1,3)=Z(IA,IC)
C-----MAKE COORDINATES POSITIVE                                         MULT0028
            DO 1 I=1,3
 1            IF (XYZUNIQUE(1,I).LT.0.0d0) XYZUNIQUE(1,I)=
     &              XYZUNIQUE(1,I)+1.0d0

            FMULT(IA,IC)=0.         ! initialize multiplicity

C-----LOOP OVER EQUIVALENT POSITIONS                                    MULT0032
            DO 4 J=1,NS             ! NS from either SPCGRP or SYMTRY

C-----LOOP OVER BRAVAIS LATTICE POINTS # just calculated                MULT0034
              DO 4 NB=1,NBRA

                ANTI=1.0d0
 2              DIF=0.0d0
C-----CALCULATE ATOM COORDINATES                                        MULT0038
C loop over x,y,z   
C there is an addend of 0.00001 here that I have removed, should be
C okay for double precision
                DO 3 I=1,3
          XYZ(I)=ANTI*(TS(I,J)+FS(1,I,J)*XYZUNIQUE(1,1)+
     &               FS(2,I,J)*XYZUNIQUE(1,2)+
     &               FS(3,I,J)*XYZUNIQUE(1,3)) ! +0.00001
C  calculate the lattice Bravais vectors (TSB)

C  first Bravais lattice point, make TSB 0 vector
                  IF (NB.EQ.1) TSB(I)=B1(I,1)
C  second Bravais lattice point, TSB is equal to IBRAVL column
                  IF (NB.EQ.2) TSB(I)=B1(I,IBRAVL)
C  fourth Bravais lattice point, TSB must be (1/2 0 1/2) (fcc)
                  IF (NB.EQ.4) TSB(I)=B1(I,6)
C  third Bravais lattice point, put in appropriate TSB
                  IF ((NB.EQ.3).AND.(IBRAVL.EQ.3)) TSB(I)=B1(I,8)  ! last of R
                  IF ((NB.EQ.3).AND.(IBRAVL.EQ.4)) TSB(I)=B1(I,5)  ! 3rd of F
C-----ADD BRAVAIS LATTICE VECTORS                                       MULT0047
                  XYZ(I)=XYZ(I)+TSB(I)
C-----MAKE COORDINATES POSITIVE AND LESS THAN 1.                        MULT0049
                  IF (XYZ(I).GE.2.0d0) XYZ(I)=XYZ(I)-1.0d0
                  IF (XYZ(I).GE.1.0d0) XYZ(I)=XYZ(I)-1.0d0
                  IF (XYZ(I).LT.-1.0d0) XYZ(I)=XYZ(I)+1.0d0
                  IF (XYZ(I).LT.0.0d0) XYZ(I)=XYZ(I)+1.0d0
 3              continue

                do i=1,3
                   if ((xyz(i).ge.1.0d0).or.(xyz(i).lt.0.0d0)) 
c     &               write(7,'(a,f6.3,a2,i1,a)') 'Trouble at ',xyz(i),
     &               write(unitstruct,'(a,f6.3,a2,i1,a)') 
     &			'Trouble at ',xyz(i),
     &               '. ',i,'-th coordinate not in unit cell.'
                enddo
C-----CALCULATE NUMBER FMULT OF COINCIDING ATOMS                        MULT0053

C                nsite=nsite+1
C                write(*,17) 'Atom #',ia,' of atom kind #',ic,
C     &           ' at ',nsite,'-th position ',xyz(1),xyz(2),xyz(3)
C 17             format(a,i1,a,i1,a,i3,a,f4.2,1x,f4.2,1x,f4.2)

                uniflag=.true.
                do i=1,nunique
                   DIF=XYZ(1)-XYZUNIQUE(i,1) ! x-coordinates
                   DIF=DABS(DIF)
                   DIF1=XYZ(2)-XYZUNIQUE(i,2) ! y-coordinates
                   DIF1=DABS(DIF1)
                   DIF=DIF+DIF1
                   DIF1=XYZ(3)-XYZUNIQUE(i,3) ! z-coordinates
                   DIF1=DABS(DIF1)
                   DIF=DIF+DIF1
C  increase multiplicity by 1.
                   IF ((DIF.LT..001).and.(i.eq.1))
     &                  FMULT(IA,IC)=FMULT(IA,IC)+1.
                   if (dif.lt..001) uniflag=.false.
                enddo

                if (uniflag) then
                   nunique=nunique+1
                   do i=1,3
                      xyzunique(nunique,i)=xyz(i)
                   enddo
                endif
C-----CHECK CENTROSYMMETRICALLY RELATED ATOMS                           MULT0064
                IF (ISYMCE.NE.1) GO TO 4   ! exit loop if not centrosym
                IF (ANTI.EQ.-1.0d0) GO TO 4
                ANTI=-1.0d0
                GO TO 2                    ! change sign and count again
 4          CONTINUE
            IF (FMULT(IA,IC).NE.0.) GO TO 5
            A=-1.                          ! Problem!
            WRITE (IPR,8) IA,IC
 5          FMULT(IA,IC)=1./FMULT(IA,IC)
c            write(7,18) elemt(ic),'#',ia,
            write(unitstruct,18) elemt(ic),'#',ia,
     &        ' has multiplicity ',int(1./fmult(ia,ic)+.001),
     &        ', scattering length ',fneu,' and occupancy ',
     &          foccu(ia,ic)
 18         format(a3,a1,i1,a,i3,a,f5.3,a,f4.2,1h.)
c            write(7,19) 'There are ',nunique,' positions.'  ! at:'
            write(unitstruct,19) 'There are ',nunique,' positions.'  ! at:'
C  output as two columns to save space
C  divide by 2 and round down
            do i=1,nunique
              if (2*(i/2).ne.i) then
c                write(7,20) (xyzunique(i,j),j=1,3)
                write(unitstruct,20) (xyzunique(i,j),j=1,3)
              else
c                write(7,21) (xyzunique(i,j),j=1,3)
                write(unitstruct,21) (xyzunique(i,j),j=1,3)
              endif
               bpos(npos+i)=fneu*foccu(ia,ic)
               do j=1,3
                 xyzpos(npos+i,j)=xyzunique(i,j)
               enddo
            enddo
c            if (2*(nunique/2).ne.nunique) write(7,*)           
            if (2*(nunique/2).ne.nunique) write(unitstruct,*) 
 19         format(a,i3,a)
 20         format(f7.5,1x,f7.5,1x,f7.5,$)
 21         format(7x,f7.5,1x,f7.5,1x,f7.5)            
C 20         format(f5.3,1x,f5.3,1x,f5.3,$)
C 21         format(7x,f5.3,1x,f5.3,1x,f5.3)            
            npos=npos+nunique
 6        CONTINUE
 7      CONTINUE
        write(unitstruct,*)
c        write(7,*)

      RETURN
C                                                                       MULT0077
 8    FORMAT (22H WRONG MULTIPLICITY OF,I5,16H TH ATOM OF TYPE,I5,43H CH
     1ECK COORDINATES AND EQUIVALENT POSITIONS)
      END
