%
% Program SYMPRJ gives the symmetry adaptation of orbitals.
%
%
%
%
%                      ******************************
%                      *                            *
%                      *                            *
%                      *          SYMPRJ            *
%                      *                            *
%                      *        Version 1.0         *
%                      *                            *
%                      *        2008-02-26          *
%                      *                            *
%                      *  Author: Per-Olof Jansson  *
%                      *                            *
%                      *  SYMPRJ caculates the      *
%                      *  symmetry projection       *
%                      *  matrices for single space *
%                      *  groups. SYMPRJ Matlab     *
%                      *  code is functionally,     *
%                      *  largely based on the      *
%                      *  SYMPRJ, FORTRAN program.  *
%                      *  SYMPRJ calculates the     *
%                      *  symmetry projection       *
%                      *  matrices for single space *
%                      *  groups, and is written by *
%                      *  Esko Blokker and Stig     *
%                      *  Flodmark. QCPE 259        *
%                      *                            *
%                      ******************************
%
%
%
% Section 1
%
nopi1 = 1;
nopli1 = 1;
sil = 1;
ksym = 1;
ntz = 0;
%
% Input
%
read_inputdata = fopen('InSb.dat','r');
read_npgo = fopen('npgo','r');
read_nge = fopen('nge','r');
read_NP = fopen('NP','r');
read_MOh = fopen('MOh','r');
read_MD6h = fopen('MD6h','r');
read_rgr3 = fopen('rgr3','r');
read_ldrmm = fopen('ldrmm','r');
%
% Output
%
fid = fopen('InSb_res.m','a');
write_proj = fopen('InSb_proj','a');
%
% Section 1.1
%
steer = fscanf(read_inputdata,'%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d',[1,20]);
%
% Section 1.2
%
a = fscanf(read_inputdata,'%g %g %g',[3,3]);
a = a.';
for I = 1:72
  for K1 = 1:84
    for K2 = 1:2
      ldrmmorg(K2,I,K1) = fread(read_ldrmm, 1,'float64');
    end
    ldrmm(I,K1) = ldrmmorg(1,I,K1) + i*ldrmmorg(2,I,K1);
  end
end
%
% Section 1.3
%
fprintf(fid,'Input, crystal with unit cell vectors (X, Y, Z)*\n');
for I = 1:3
  fprintf(fid,'(%9.6f %9.6f %9.6f)\n', a(I,1:3));
end
fprintf(fid,'= (A1, A2, A3)\n\n');
T = 2*pi;
%
% Section 1.4
%
b = (inv(a)).';
bi = a.';
ai = b.';
%
% Section 1.5
%
fprintf(fid,'Reciprocal unit cell vectors are (X, Y, Z)*\n');
for I = 1:3
  fprintf(fid,'(%9.6f %9.6f %9.6f)\n', b(I,1:3));
end
fprintf(fid,'= (B1, B2, B3)\n\n');
%
% Section 1.6
%
pgnr = fscanf(read_inputdata,'%3d',1);
nel = fscanf(read_inputdata,'%3d',1);
%
% lmax is the maximum value of orbital quantum number L.
%
% This version of SYMPRJ has the option; different lmax for different atoms.
%
for I = 1:nel
  lmax(I) = fscanf(read_inputdata,'%3d', 1);
end
%
% Section 1.7
%
fprintf(fid,'The unit cell contains %3d chemical elements\n', nel);
%
% pgnr is the index number of the point group, 1=C1, 2=S2, 3 = C2,
% 4=Cs, 5=C2h, 6=D2, 7=C2v, 8=D2h, 9=C4, 10=S4, 11=C4h, 12=D4, 13=C4v,
% 14=D2d, 15=D4h, 16=C3, 17=S6, 18=D3 (for space groups 149, 151, 153),
% 19=D3 (for space groups 150, 152, 154, 155), 20 = C3v (for space
% groups 156, 158, 160, 161), 21 = C3v (for space groups 157, 159),
% 22 = D3d (for space groups 164, 165, 166, 167), 23 = D3d
% (for space groups 162, 163) 24 = C6, 25 = C3h, 26= C6h,
% 27 = D6, 28 = C6v, 29 = D3h (for space groups 187, 188), 30 = D3h
% (for space groups 189, 190), 31 = D6h, 32 = T, 33 = Th, 34 = O,
% 35 = Td, 36 = Oh
%
tsmall = 1*10^(-5);
ttsmall = 1*10^(-6);
%
% Section 1.8
%
for I = 1:nel
  nat(I) = fscanf(read_inputdata,'%3d', 1);
end
%
% nat(I) is the number of atoms of chemical elements I, per unit cell.
%
%
% Section 1.9
%
for I = 1:nel
  K = nat(I);
  for J = 1:K
    atco = fscanf(read_inputdata,'%3d',1);
    r(I,J,1:3) = fscanf(read_inputdata,'%g %g %g',3);
%
% atco=1 means cartesian coordinates, atco=0 means lattice coordinates.
%
    if atco ~= 1
      tsk(1:3) = r(I,J,1:3);
      r(I,J,1:3) = a(1:3,1:3)*tsk(1:3).';
    end
  end
end
%
% Section 1.10
%
for I = 1:nel
  K = nat(I);
  fprintf(fid,'Of chemical element no %3d there are %3d atoms, with positions:\n\n', I, K);
  for J = 1:K
    fprintf(fid,'Atom %3d (%7.4f %7.4f %7.4f)\n', J, r(I,J,1:3));
  end
end
fprintf(fid,'All atom positions are printed in cartesian coordinates\nAll wave vectors will be printed as coordinates of the\nreciprocal lattice vectors b1, b2, b3, divided by 2*pi\n\n'); 
%
% r(I,J,1:3) are the coordinates of atom J of chemical element I,
% expressed in cartesian coordinates.
%
for I = 1:2
  for J = 1:36
    npgo(I,J) = fread(read_npgo,1,'int16');
  end  
end
%
% Section 1.11
%
order = npgo(1,pgnr);
first = npgo(2,pgnr);
%
% Section 1.12
%
nge(1:368) = fread(read_nge,368,'int16');
gel(1:order) = nge(first:(first+order-1));
npri(1:100) = fread(read_NP,100,'int16');
if ~steer(20)
%
% Section 1.13
%
  for I = 1:order
%
% The nonprimitive translations.
%
    uco = fscanf(read_inputdata,'%2d',1);
    for K1 = 1:3
      u(K1,I) = fscanf(read_inputdata,'%g',1);
    end
%
% uco=1 means cartesian coordinates, uco=0 means lattice coordinates.
%
    if uco ~= 0
      tsk(1:3) = u(1:3,I);
      u(1:3,I) = ai(1:3,1:3)*tsk(1:3).';
%
% Nonprimitive translations are calculated in lattice coordinates.
%
    end
  end
  u = u.';
else
  u = zeros(order,3);
end
%
% Section 1.14
%
if ~((pgnr >= 16) & (pgnr <= 31))
%
% Section 1.15
%
  K48 = 0;
  fprintf(fid,'The pointgroup (no. %3d) of the crystal is a subgroup of Oh, with element numbers:\n',pgnr);
  for I = 1:48
    for J = 1:48
      mtab(I,J) = fread(read_MOh,1,'int8');
    end
  end
else
%
% Section 1.16
%
  for I = 1:24
    for J = 1:24
      mtab(I,J) = fread(read_MD6h,1,'int8');
    end
  end
  K48 = 48;
  fprintf(fid,'The pointgroup (no. %3d) of the crystal is a subgroup of D6h, with element numbers:\n',pgnr);
end 
if (pgnr ~= 31) & (pgnr ~= 36)
%
% Section 1.17
%
  mtab(1:order,1:order) = mtab(gel(1:order),gel(1:order));
  inver(gel(1:order)) = 1:order;
  mtab(1:order,1:order) = inver(mtab(1:order,1:order));
end
%
% inver gives the numbering of the elements in the specific point
% group. Example: in the group 5=C2h, element 28 is the 4th
% element of the group, so inver(28)=4.
%
%
% Section 1.18
%
fprintf(fid,'%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', gel(1:order));
fprintf(fid,'\n\n');
fprintf(fid,'The maximum value for the orbital quantum number L is set to: %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', lmax(1:nel));
%
% Section 1.19
%
%
% Load the rotation matrices for the orthogonal coordinate system.
%
for I = 1:72  
  for K1 = 1:3
    rgr3(K1,1:3,I) = fread(read_rgr3,3,'float64');
  end
end
rgr(1:3,1:3,1:order) = rgr3(1:3,1:3,gel(1:order) + K48);
%
% Section 1.20
%
number_of_wave_vectors = fscanf(read_inputdata,'%2d',1);
fwrite(write_proj, number_of_wave_vectors, 'int8');
fwrite(write_proj, nel, 'int8');
fwrite(write_proj, nat(1:nel), 'int8');
fwrite(write_proj, lmax(1:nel), 'int8');
last = fscanf(read_inputdata,'%2d',1);
wvco = fscanf(read_inputdata,'%2d',1);
for I = 1:3
  rk(I) = fscanf(read_inputdata,'%g',1);
end
nfacto = fscanf(read_inputdata,'%2d',1);
steer(11) = 1;
while last ~= 1
%
% last=1 for the last k-vector in the input, otherwise last=0.
% wvco=1 means cartesian coordinates; wvco=0 means reciprocal 
% lattice coordinates. rk(1:3) is a k-vector, expressed in the
% reciprocal unit vectors. nfacto is not zero, if the user wants
% to make calculations for other k-vectors with the same
% direction, but different length. Then nfacto is the number
% of such vectors to be treated after the first k-vector in that
% direction.
%
  if wvco == 1
    tsk = rk;
    rk(1:3) = bi(1:3,1:3)*tsk(1:3).';
%
% Wave vector in reciprocal lattice coordinates.
%
  end
  if nfacto > 0
    factor = fscanf(read_inputdata,'%g %g %g %g %g %g %g %g', [1,nfacto]);
  end
  ark(1:3) = rk(1:3);
%
% Section 1.21
%
  factor(1) = 1;
%
% Section 1.22
%
  IV = 1;
  while (IV <= 1) | (IV <= nfacto)
%
% Section 1.23
%
    if IV <= 2
%
% Section 2
%
% Formation of the point group of the wave vector.
%
%
% Section 2.1
%
      srk(1:3) = ark(1:3)*factor(IV);
      rk(1:3) = srk(1:3)*T;
      if ~((rk(1) == 0) & (rk(2) == 0) & (rk(3) == 0))
%
% Section 2.3
%
        kgord = 1;
        kgel(1) = 1;
        brk(1:3) = b(1:3,1:3)*rk(1:3).';
%
% brk is the wave vector expressed in cartesian coordinates
%
        for K = 2:order
          robrk(1:3) = rgr(1:3,1:3,K)*brk(1:3).';
          rosk(1:3) = bi(1:3,1:3)*robrk(1:3).';
%
% rosk is the rotated wave vector in reciprocal lattice coordinates.
%
          rosk(1:3) = rosk(1:3) - rk(1:3);
%
% This difference should be a big k-vector, if the element belongs to the
% group of the k-vector.
%
          K1 = 1;
          while K1 <= 3
            TT = abs(rosk(K1));
            TTT = abs(TT-T);
            if (TT > tsmall) & (TTT > tsmall)
%
% If this condition is fulfilled, rosk is not a big k-vector.
%
              break;
            end
            K1 = K1 + 1;
          end
          if K1 > 3
            kgord = kgord + 1;
            kgel(kgord) = K;
          end
        end
        kkgel(1:kgord) = gel(kgel(1:kgord));
        kg = kgord;
%
% The Kth element of the group of the wave vector is the kkgel(K)th
% element of the complete set of point group operators,
%
        mtab2(1:kgord,1:order) = mtab(kgel(1:kgord),1:order);
        mtab2(1:kgord,1:kgord) = mtab2(1:kgord,kgel(1:kgord));
        inverk(kgel(1:kgord)) = 1:kgord;
        mtab2(1:kgord,1:kgord) = inverk(mtab2(1:kgord,1:kgord));
%
% Section 3
%
% Tests for the nonsymmorphic space group.
%

        if ~steer(20)
%
% Section 3.1
%
%
% Test if the k-vector lies within the first Brillouin zone.
%
          ibz = 0;
          ntz = bztest(srk,b);
          if ntz > 0
%
% k-vector outside Bz
%
            fprintf(fid,'(%7.4f %7.4f %7.4f) outside the Bz\n\n', rk(1:3));
            break;
          end
          if ntz < 0
            ibz = 1;
            if steer(18)
              fprintf(fid,'Nonsymmorphic, but within Bz\n\n');
            end
          end
%
% Section 3.2
%
%
% Test if the group of the k-vector is symmorphic.
%
          ksym = 0;
          I = 1;
          while I <= kgord
            I1 = kgel(I);
            K = 1;
            while K <= 3
              if u(I1,K) ~= 0
                break;
              end
              K = K + 1;
            end
            if K <= 3
              break;
            end
            I = I + 1;
          end
          if I > kgord
            ksym = 1;
            if steer(18)
              fprintf(fid,'Nonsymmorphic, but symmorph Gk\n\n');
            end
          end
            
%
% Section 4
%
% Formation of the factor group Gk/Tk
%
% Here we treat the case of a nonsymmorphic group of the k-vector, with the
% k-vector on one of the sides of the Brillouin zone boundaries. First, initiate
% the tables used to construct the factor group Gk/Tk. k2gord will be the order 
% of this group. listp(I) contains the index of the pointgroup operator of the
% Ith element of this group. Several elements may have the same pointgroup
% operator, so that they differ only in the translation part til(I,K), K = 1:3
% and the corresponding exponential part sil(I). nopi(I) is the number of
% elements with the same Ith pointgroup operator. The Kth operator, which has
% I as its pointgroup operator is group element number nopli(I,K) in the factor
% group. nr(I) is the number of elements in the Ith row of the group
% multiplication table mtab3(I,K), which has been calculated. At the end all
% nr(I) = k2gord.
%
          if (~ibz) & (I <= kgord)
            k2gord = kgord;
            for I = 1:k2gord
              listp(I) = I;
              K1 = kgel(I);
              til(I,1:3) = u(K1,1:3);
              nopi(I) = 1;
              nr(I) = 1;
              nopli(I,1) = I;
              mtab3(1,I) = I;
              mtab3(I,1) = I;
              ar = til(I,1:3)*rk(1:3).';
              sres = -i*ar;
              sil(I) = exp(sres);
            end
            nr(1) = k2gord;
            I = 2;
            while I <= k2gord
              K = nr(I) + 1;
              N1 = listp(I);
              N3 = kgel(N1);
              while K <= k2gord
                N2 = listp(K);
                index1 = mtab2(N1,N2);
                tsk(1:3) = a(1:3,1:3)*til(K,1:3).';
%
% Transform to Cartesian coordinates.
%
                rprod(1:3) = rgr(1:3,1:3,N3)*tsk(1:3).';
                tsk(1:3) = rprod(1:3);
                for K1 = 1:3
                  rprod(K1) = ai(K1,1:3)*tsk(1:3).';
                  rprod(K1) = rprod(K1) + til(I,K1);
                end
                ar = rk(1:3)*rprod(1:3).';
                R4 = -i*ar;
                sres = exp(R4);
                NT = nopi(index1);
                K1 = 1;
                while K1 <= NT
                  indk1 = nopli(index1,K1);
                  if abs(sil(indk1) - sres) <= 5*10^(-2)
                    mtab3(I,K) = indk1;
                    nr(I) = nr(I) + 1;
                    break;
                  end
                  K1 = K1 + 1;
                end
                if K1 > NT
                  k2gord = k2gord + 1;   
                  mtab3(I,K) = k2gord;
                  nr(I) = nr(I) + 1;
                  mtab3(1,k2gord) = k2gord;
                  nr(1) = k2gord;
                  mtab3(k2gord,1) = k2gord;
                  nr(k2gord) = 1;
                  listp(k2gord) = index1;
                  til(k2gord,1:3) = rprod(1:3);
                  sil(k2gord) = sres;
                  nopi(index1) = nopi(index1) + 1;
                  K2 = nopi(index1);
                  nopli(index1,K2) = k2gord;
                  if I ~= 2
                    I = 2;
                    K = nr(I);
                    N1 = listp(I);
                    N3 = kgel(N1);
                  end
                end
                K = K + 1;
              end
              I = I + 1;
            end
            kgord = k2gord;
            nopi1 = nopi(1);
            nopli1(1:nopi1) = nopli(1,1:nopi1);
            llistp(1:kgord) = kkgel(listp(1:kgord));
            mtab2(1:kgord,1:kgord) = mtab3(1:kgord,1:kgord);
          end
        end
      else
%
% Section 2.2
%
%
% Section for the case of zero wave vector. Then the point group of the wave
% vector is equal to the point group of the space group.
%
        kgord = order;
        mtab2(1:kgord,1:kgord) = mtab(1:kgord,1:kgord);
        kgel(1:kgord) = 1:kgord;
        kkgel(1:kgord) = gel(1:kgord); 
        ibz = 1;
        kg = kgord;
      end
%
% Section 5
%
% Wave vectors in the same direction. Here we treat the case with the k-vector 
% in the same direction as the previous one.
%
    else
      srk(1:3) = ark(1:3)*factor(IV);
      rk(1:3) = srk(1:3)*T;
    end
%
% Section 6
%
% Calculation of the diagonal elements of the irreducible representations.
%
%
% Section 6.1
%
    fprintf(fid,'\n\n');
    fprintf(fid,'Projection matrices for the wave vector (%7.4f %7.4f %7.4f)\n', srk(1:3));
    fprintf(fid,'The pointgroup of the wave vector consists of %3d operarators, indexed as numbers:\n', kg);
    fprintf(fid,'%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', kkgel(1:kg));
    if ~(steer(20) | ksym | ibz)
      fprintf(fid,'\n\n');
      fprintf(fid,'The factor group Gk/Tk consists of\n') ;
      for I = 1:kgord
        fprintf(fid,'%3d pointgroup operator: %3d  nonprimitive translation: (%7.4f %7.4f %7.4f), exp = (%7.4f %7.4f)\n', I, llistp(I), til(I,1:3), real(sil(I)), imag(sil(I)));
      end
      fprintf(fid,'\n');
    end
%
% Now we have determined the point group of the wave vector, or (for 
% the nonsymmorphic case with wave vector the BZ-boundary) the factor
% group Gk/Tk. Next we determine all (allowable) irreducible
% representations of this finite group.
%
    if (IV <= 2) | (steer(20) | ksym | ibz)
%
% Section 6.2
%
      [jpdd,allow,ncl,laj] = Irrep(kgord,mtab2,npri,steer,ibz,ksym,nopi1,nopli1,sil);
      if ~steer(11)
        break;
      end
%
% Section 6.3
%
      nup = 0;
      nip = 0;
      for I = 1:ncl
        nip = nip + laj(I);
        if allow(I)
          nup = nup + laj(I);
        end
      end
      if nip ~= nup
        nal = 0;
        for I = 1:ncl
          if allow(I)
            nal = nal + 1;
            nalr(nal) = I;
          end
        end
        fprintf(fid, '\n\n');
        fprintf(fid,'Allowable irreps for J = %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', nalr(1:nal));
      end
      nblock = 0;
      for ilmax = 1:nel
        nblock = nblock + lmax(ilmax) + 1;
      end
      fprintf(fid,'\n');
      fprintf(fid, 'There are %3d projection matrices for this wave vector, with index J for the irreps and index JD for the diagonal elements of the irrep.\n\n', nup);
      fprintf(fid, 'Each projection matrix is blockdiagonalised with respect to the indices ichem = 1, %2d , and L = %1d %1d %1d %1d %1d %1d %1d %1d %1d %1d, so each projection matrix consists of %3d subblocks of sub-projection matrices along the main diagonal. These are independently orthonormalised into sub-T-matrices, from which the complete T-matrix can be constructed.\n\n', nel, lmax(1:nel), nblock);                   
%
% Section 7
%
% Formation of the summation sets G(MU,NU). nup is the number of projection
% matrices for this wave vector. We form np(K,I,J), which is the number
% of elements of the group of the k-vector (groupk), for which the difference
% rgr*(coordinate of atom I) - (coordinate of atom J) of chemical element K,
% is a lattice vector. npl(1:nel,1:nat(I1),1:nat(I1),1:np(1:nel,1:nat(I1),
% 1:nat(I1))) gives the indices of these pointgroup operators.
% nvec(1:nel,1:nat(I1),1:nat(I1),1:np(1:nel,1:nat(I1),1:nat(I1)),1:3)
% gives the corresponding lattice vectors.
%
      I1 = 0;
      while I1 < nel
        I1 = I1 + 1;
%
% I1 is here the index of the chemical element
%
        JJ = nat(I1);
%
% JJ is the number of atoms of chemical element I1.
%
        I3 = 0;
        while I3 < JJ
          I3 = I3 + 1;
%
% I3 is the atom row index.
%
          np(I1,I3,1:JJ) = zeros(1,JJ);
          I4 = 0;
          while I4 < kgord
            I4 = I4 + 1;
%
% I4 is the index of elements of groupk
%
            I10 = I4;
%
% Inverse point group operator.
%
            if ~(steer(20) | ksym | ibz)
              I5 = listp(I10);
              I5 = kgel(I5);
              trac(1:3) = a(1:3,1:3)*til(I4,1:3).';
%
% Nonprimitive translations in Cartesian coordinates.
%
            else
              I11 = kgel(I4);
              I5 = kgel(I10);
              trac(1:3) = a(1:3,1:3)*u(I11,1:3).';
            end
            for I7 = 1:3
              for I8 = 1:3
%
% Space group transformation of the atom position vector.
%
                trac(I7) = trac(I7) + rgr(I7,I8,I5)*r(I1,I3,I8);
              end
            end
            I9 = 1;
            while I9 <= JJ
%
% I9 is the atom column index.
%
              for I12 = 1:3
                dif(I12) = trac(I12) - r(I1,I9,I12);
              end
              difi(1:3) = ai(1:3,1:3)*dif(1:3).';
%
% Transform to lattice coordinates. Test if difi is a lattice vector.
%
              for I10 = 1:3
                isign = 1;
                if difi(I10) < 0
                  isign = -1;
                end
                D = isign*difi(I10);
                ifd = fix(D);
                if (D - ifd) >= 0.98
                  ifd = ifd + 1;
                end
                if abs(D-ifd) > 0.02
                  break;
                end
                difi(I10) = isign*ifd;
              end
              if abs(D-ifd) <= 0.02
                np(I1,I3,I9) = np(I1,I3,I9) + 1;
                K = np(I1,I3,I9);
                npl(I1,I3,I9,K) = I4;
                nvec(I1,I3,I9,K,1:3) = difi(1:3);
                if steer(18)
                  fprintf(fid,'np(%3d %3d %3d) = %3d npl = %3d nvec = (%7.4f %7.4f %7.4f)\n',I1, I3, I9, K, I4, difi(1:3));
                end
                break;
              end
              I9 = I9 + 1;
            end
            if I9 > JJ
              fprintf(fid,'Wrong space group %3d %3d %3d %3d\n',I1, I3, I4 ,I5);
              break;
            end
          end
          if I9 > JJ
            break;
          end      
        end
        if I9 > JJ
          break;
        end   
      end
    end
    if I9 > JJ
      break;
    end   
%
% Section 8
%
% The formation of the projection matrices. Here starts the loops necessary to
% form the projection matrix.
%
%
% Section 8.1
%
    J = 0;
    N31 = -1;
    while J < ncl
      J = J + 1;          
%
% J is the index of the irreducible representation
%
      if allow(J)
        for JD = 1:laj(J)
          fprintf(fid, '\n\n');
          fprintf(fid,'                             T-matrix for J = %3d, JD = %3d\n\n', J, JD);
          nb = 1;
%
% JD is the index of the diagonal elements of the irreducible representation J. 
%
          for I = 1:kgord
            jpdd(J,JD,I) = conj(jpdd(J,JD,I));
          end
          if steer(18)
            fprintf(fid,'jpdd = \n');
            for I = 1:kgord
              jpddorg(1,JD,I) = real(jpdd(J,JD,I));
              jpddorg(2,JD,I) = imag(jpdd(J,JD,I));
            end
            fprintf(fid,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n',jpddorg(1:2,JD,1:kgord));
            fprintf(fid,'\n\n');
          end
%
% ichem is the index of the chemical element.
%
          for ichem = 1:nel
            lsqsum = 0;
            for L = 0:lmax(ichem)
%
% L is the orbital quantum number.
%
              if L > 0
                lsqsum = lsqsum + (2*L - 1)^2;
              end
              N = 2*L + 1;
              for mu1 = 1:nat(ichem)
%
% mu1 is the atom index, row index in the projection matrix.
%
                for mu2 = 1:nat(ichem)
%
% mu2 is the atom index, column index in the projection matrix.
%
                  ncoset = np(ichem,mu1,mu2);
%
% ncoset is the order of the little coset.
%
                  if ncoset ~= 0
                    for I = 1:ncoset
                      for J1 = 1:3 
                        nrn(I,J1) = nvec(ichem,mu1,mu2,I,J1);
                      end
                    end
                  end
%
% Section 8.2
%
                  for M1 = -L:L
                    N1 = M1 + L + 1;
                    NN1 = N1 + (mu2 - 1)*N;
                    for M2 = -L:L
                      N2 = M2 + L + 1;
                      NN2 = N2 + (mu1 - 1)*N;
                      jdpk(NN1,NN2) = 0;
                      if ncoset ~= 0
                        N3 = lsqsum + (N1 - 1)*N + N2;
                        if N3 ~= N31
                          N31 = N3;
                          ldmm(1:72) = ldrmm(1:72,N3);
                          if steer(18)
                            fprintf(fid,'ldmm\n');
                            for I = 1:72
                              fprintf(fid,'(%7.4f %7.4f)\n',real(ldmm(I)), imag(ldrmm(I)));
                            end
                          end
                        end
                        for I = 1:ncoset
                          K3 = npl(ichem,mu1,mu2,I);
                          if  ~(steer(20) | ksym | ibz)
                            K = kkgel(listp(K3));
                            R5 = 0;
                          else
                            K = kkgel(K3);
                            KI = kgel(K3);
                            R5 = i*(u(KI,1:3)*rk(1:3).');
                          end
                          K = K + K48;
                          R4 = -i*(nrn(I,1:3)*rk(1:3).');
                          R4 = R4 + R5;
                          ep = exp(R4);
%
% Section 8.3
%
                          R4 = jpdd(J,JD,K3)*ep*ldmm(K);
                          jdpk(NN1,NN2) = jdpk(NN1,NN2) + R4;
                        end
                      end
                    end
                  end
                end
              end
              ndi = N*nat(ichem);
              rh = laj(J)/kgord;
              jdpk(1:ndi,1:ndi) = jdpk(1:ndi,1:ndi)*rh;
%
% Section 8.4
%
%
% Test projection matrix
%
              jdprod(1:ndi,1:ndi) = jdpk(1:ndi,1:ndi)*jdpk(1:ndi,1:ndi);
              tmatri(1:ndi,1:ndi) = jdprod(1:ndi,1:ndi) - jdpk(1:ndi,1:ndi);
              K4 = 1;
              while K4 <= ndi
                K5 = 1;
                while K5 <= ndi
                  if abs(tmatri(K4,K5)) > tsmall
                    break;
                  end
                  K5 = K5 + 1;
                end
                if K5 <= ndi
                  break;
                end
                K4 = K4 + 1;
              end
%
% Error, no projection matrix for J.
%                     
              if K4 <= ndi
                fprintf(fid,'Error, no projection matrix for J = %3d, JD = %3d, laj = %3d, ichem = %3d, L = %3d, ndi = %3d\n\n', J,JD,laj(J),ichem,L,ndi);
                jdpkorg(1,1:ndi,1:ndi) = real(jdpk(1:ndi,1:ndi));
                jdpkorg(2,1:ndi,1:ndi) = imag(jdpk(1:ndi,1:ndi));
                for K4 = 1:ndi
                  fprintf(fid,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n', jdpkorg(1:2,K4,1:ndi));
                end
                fprintf(fid,'\n');
                fprintf(fid,'jdpk*jdpk - jdpk\n');
                tmatriorg(1,1:ndi,1:ndi) = real(tmatri(1:ndi,1:ndi));
                tmatriorg(2,1:ndi,1:ndi) = imag(tmatri(1:ndi,1:ndi));
                for K4 = 1:ndi
                  fprintf(fid,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n', tmatriorg(1:2,K4,1:ndi));
                end
                fprintf(fid,'\n\n');
              end
              tmatri(1:ndi,1:ndi) = jdpk(1:ndi,1:ndi) - jdpk(1:ndi,1:ndi)';
              K4 = 1;
              while K4 <= ndi
                K5 = 1;
                while K5 <= ndi
                  if abs(tmatri(K4,K5)) > tsmall
                    break;
                  end
                  K5 = K5 + 1;
                end
                if K5 <= ndi
                  break;
                end
                K4 = K4 + 1;
              end
              if K4 <= ndi
                fprintf(fid,'jdpk - hermite conjugate(jdpk)\n');
                fprintf(fid,'Error, no projection matrix for J = %3d, JD = %3d, laj = %3d\n', J,JD,laj(J));
                fprintf(fid,'\n');
                tmatriorg(1,1:ndi,1:ndi) = real(tmatri(1:ndi,1:ndi));
                tmatriorg(2,1:ndi,1:ndi) = imag(tmatri(1:ndi,1:ndi));
                for K4 = 1:ndi
                  fprintf(fid,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n', tmatriorg(1:2,K4,1:ndi));
                end
                fprintf(fid,'\n\n');
              end
%
% Section 8.5
%
%
% We orthonormalise the submatrix (fixed L, ichem) and store the resulting
% sub-T-matrix.
%
              ptrace = sum(diag(jdpk(1:ndi,1:ndi)));
              rntr = real(ptrace);
              ntr = fix(rntr);
              if (rntr - ntr) > 0.5
                ntr = ntr + 1;
              end
              if ntr ~= 0
%
% ntr is the trace, the number of linearly independent columns in the submatrix.
% First we select columns with the diagonal term equal to 1, since these are
% automatically orthonormal to all other columns. Skip all columns with diagonal
% term equal to 0, since these are zero columns.
%
                nspec = zeros(1,ndi);
%
% nspec(I) registrates these special columns.
%
                itotal = 0;
                for I = 1:ndi
                  if (abs(jdpk(I,I) - 1))^2 <= ttsmall
                    itotal = itotal + 1;
                    for N = 1:ndi
                      tmatri(N,itotal) = 0;
                    end
                    tmatri(I,itotal) = 1;
                    nspec(I) = 1;
                  else
                    if (abs(jdpk(I,I)))^2 <= ttsmall
                      nspec(I) = -1;
                    end
                  end
                end
                I1 = itotal + 1;
                I = 1;
                if itotal ~= ntr
                  while I <= ndi
                    if nspec(I) == 0
                      ito = itotal;
                      itotal = itotal + 1;
                      sumtot = sum((abs(jdpk(1:ndi,I))).^2);
                      if sumtot >= ttsmall
                        sumtot = 1/sqrt(sumtot);
                        tmatri(1:ndi,itotal) = jdpk(1:ndi,I)*sumtot;
                        if itotal ~= I1
                          NC = I1;
                          while NC <= ito
%
% Orthonormalisation against previously found orthonormal columns, which have
% not one element equal to 1.
%
                            R4 = tmatri(1:ndi,NC)'*tmatri(1:ndi,itotal);
                            if abs((abs(R4))^2 - 1) > 0.01
                              tmatri(1:ndi,itotal) = tmatri(1:ndi,itotal) - R4*tmatri(1:ndi,NC);
                              R1 = 0;
                              LD1 = 1;
                              while LD1 <= ndi
                                if abs(tmatri(LD1,itotal)) >= 0.01
%
% If all elements of the column become very small after an orthonormalisation, 
% we do not use the column. Renormalisation would give large errors.
% 
                                  break;
                                end
                                LD1 = LD1 + 1;
                              end
                              if LD1 <= ndi
                                R1 = sum((abs(tmatri(1:ndi,itotal))).^2);
                                R1 = 1/sqrt(R1);
%
% Renormalisation.
%
                                tmatri(1:ndi,itotal) = tmatri(1:ndi,itotal)*R1;
                              else
                                itotal = itotal - 1;
                                break;
                              end
                            else
                              itotal = itotal - 1;
                              break;
                            end
                            NC = NC + 1;
                          end
                          if NC > ito
                            if itotal == ntr
                              break;
                            end
                          end
                        else
                          if itotal == ntr
                            break;
                          end
                        end
                      end
                    end
                    I = I + 1;
                  end
                  if itotal ~= ntr
                    fprintf(fid,'Error, not enough orthonormal columns');
                    break;
                  end
                end
              else
                itotal = 0;
              end
%
% Section 8.6
%
% Here starts the output section. The projection matrices (T-matrices) are
% printed. Thereafter all of them have been calculated for all values of J,
% JD, ichem, L.
%
              fprintf(fid,'\n');
              fwrite(write_proj,ichem, 'int8');
              fwrite(write_proj,L, 'int8');
              fwrite(write_proj,itotal, 'int8');
              fwrite(write_proj,ndi, 'int8');
              if itotal ~= 0
                for I = 1:itotal
                  for I1 = 1:ndi
                    tmatriorg(1,I1,I) = real(tmatri(I1,I));
                    tmatriorg(2,I1,I) = imag(tmatri(I1,I));
                  end
                end
                fprintf(fid, 'Subblock %3d sub-T-matrix for ichem = %3d, L = %3d, %3d atoms per unit cell\n', nb, ichem, L, nat(ichem));
                fprintf(fid, '%3d columns, %3d rows, the rows are indexed as follows, (atom 1 M = -L,.,L) (atom 2 M = -L,.,L) etc.\n', itotal, ndi);
                for I = 1:itotal
                  fprintf(fid, 'Column %3d\n', I);
                  fprintf(fid, '(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n', tmatriorg(1:2,1:ndi,I));
                  fprintf(fid,'\n');
                  fwrite(write_proj,tmatriorg(1:2,1:ndi,I), 'float64');
                end
                fprintf(fid,'\n');
                nb = nb + 1;
              else
                fprintf(fid, 'Subblock, %3d sub-T-matrix for ichem = %3d, L = %3d, %3d atoms per unit cell is empty\n\n', nb, ichem, L, nat(ichem));
                nb = nb + 1;
              end
            end
          end
        end
      end
    end 
    IV = IV + 1;
  end
  if I9 > JJ
    break;
  end
  if ~steer(20)
    if ntz > 0
      break;
    end
  end
  steer(11) = 1;
  last = fscanf(read_inputdata,'%2d',1);
  wvco = fscanf(read_inputdata,'%2d',1);
  for I = 1:3
    rk(I) = fscanf(read_inputdata,'%g',1);
  end
  nfacto = fscanf(read_inputdata,'%2d',1);
end
%
% Section 9
%
fprintf(fid, '\n\n');
fprintf(fid, 'SYMPRJ has finished executing\n');
fclose(fid);
fclose(read_npgo);
fclose(read_nge);
fclose(read_NP);
fclose(read_inputdata);
fclose(read_MOh);
fclose(read_MD6h);
fclose(read_rgr3);
fclose(read_ldrmm);
fclose(write_proj);
fprintf('SYMPRJ has finished executing\n');