%
% Program SYMPWS gives the symmetry adaptation of plane wave orbitals.
%
%
% Section 1
%
%
%
%                      ******************************
%                      *                            *
%                      *                            *
%                      *           SYMPWS           *
%                      *                            *
%                      *        Version 1.0         *
%                      *                            *
%                      *        2006-08-16          *
%                      *                            *
%                      *  Author: Per-Olof Jansson  *
%                      *                            *
%                      *  SYMPWS caculates the      *
%                      *  symmetry projection       *
%                      *  matrices for double space *
%                      *  groups. SYMPWS Matlab     *
%                      *  code is functionally,     *
%                      *  largely based on the      *
%                      *  SYMPW, FORTRAN program.   *
%                      *  SYMPW calculates the      *
%                      *  symmetry projection       *
%                      *  matrices for single space *
%                      *  groups, and is written by *
%                      *  Esko Blokker and Stig     *
%                      *  Flodmark. QCPE 259        *
%                      *                            *
%                      ******************************
%
nopi1 = 1;
nopli1 = 1;
sil = 1;
ksym = 1;
JJJ = 1;
%
% Input
%
read_inputdata = fopen('D4_6h.dat','r');
read_spinm = fopen('spinm','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');
%
% Output
%
fid = fopen('D4_6h_res.m','a');
write_proj = fopen('D4_6h_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:144
  for K1 = 1:2
    for K2 = 1:2
      for K3 = 1:2
        spinmorg(K3,K1,K2,I) = fread(read_spinm, 1, 'float64');
      end
      spinm(K1,K2,I) = spinmorg(1,K1,K2,I) + i*spinmorg(2,K1,K2,I);
    end
  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');
%
% a(1:3,1:3) is the matrix of the unit cell vectors.
% (A1,A2,A3) = (X,Y,Z)*a(1:3,1:3).
%
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');
bia = bi*a;
aib = ai*b;
%
% Section 1.6
%
pgnr = fscanf(read_inputdata,'%3d',1);
%
% 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 group 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^(-4);
for I = 1:2
  for J = 1:36
    npgo(I,J) = fread(read_npgo,1,'int16');
  end  
end
%
% Section 1.7
%
order = npgo(1,pgnr);
first = npgo(2,pgnr);
%
% Section 1.8
%
nge(1:736) = fread(read_nge,736,'int16');
gel(1:order) = nge(first:(first+order-1));
npri(1:100) = fread(read_NP,100,'int16');
if ~steer(20)
%
% Section 1.9
%
  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.10
%
if ~((pgnr >= 16) & (pgnr <= 31))
  K96 = 0;
%
% Section 1.11
%
  fprintf(fid,'The pointgroup (no. %3d) of the crystal is a subgroup of Oh, with element numbers:\n',pgnr);
  for I = 1:96
    for J = 1:96
      mtab(I,J) = fread(read_MOh,1,'int8');
    end
  end
else
%
% Section 1.12
%
  for I = 1:48
    for J = 1:48
      mtab(I,J) = fread(read_MD6h,1,'int8');
    end
  end
  K96 = 96;
  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.13
%
  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.14
%
fprintf(fid,'%2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d %2d\n', gel(1:order));
%
% Section 1.15
%
%
% Load the rotation matrices for the orthogonal coordinate system.
%
for I = 1:144  
  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) + K96);
%
% Section 1.16
%
%
% Transform to reciprocal lattice
%
for I = 1:order
  rg1(1:3,1:3) = bi(1:3,1:3)*rgr(1:3,1:3,I);
  rgr(1:3,1:3,I) = rg1(1:3,1:3)*b(1:3,1:3);
end
%
% Section 1.17
%
number_of_wave_vectors = fscanf(read_inputdata,'%2d',1);
fwrite(write_proj,number_of_wave_vectors, '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
nrec = fscanf(read_inputdata,'%2d',1);
krep = fscanf(read_inputdata,'%2d',1);
while last ~= 1
  steer(11) = 1;
  fprintf(fid,'\n\n');
  fprintf(fid,'Projection matrices for the wave vector ( %7.4f %7.4f %7.4f)\n\n', rk(1:3));
%
% 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. If krep=1, the same reciprocal lattice
% vectors as for the previous wave vector are to be used. nrec
% is the number of reciprocal lattice vectors following this input.
% 
  if wvco == 1
    tsk = rk;
    rk(1:3) = bi(1:3,1:3)*tsk(1:3).';
%
% Wave vector in reciprocal lattice coordinates.
%
  end
  if krep ~= 1
    for I = 1:nrec
      for K1 = 1:3
        rklat(I,K1) = fscanf(read_inputdata,'%g',1);
      end
    end
    fprintf(fid,'Input reciprocal lattice vectors\n');
    for I = 1:nrec
      fprintf(fid,'%2d (%7.4f %7.4f %7.4f)\n', I, rklat(I,1:3));
    end
  end
%
% Section 2
%
% Formation of the point group of the wave vector.
%
  rkli = zeros(order,3);
%
% Test if the k-vector lies within the first Brillouin zone.
%
%
% Section 2.1
%
  ibz = 0;
  ntz = bztest(rk,b);
  if ntz > 0
    fprintf(fid,'\n');
    fprintf(fid,'(%7.4f %7.4f %7.4f), outside Bz\n\n', rk(1:3));
    last = fscanf(read_inputdata,'%2d',1);
    wvco = fscanf(read_inputdata,'%2d',1);
    for I = 1:3
      rk(I) = fscanf(read_inputdata,'%g',1);
    end
    nrec = fscanf(read_inputdata,'%2d',1);
    krep = fscanf(read_inputdata,'%2d',1);
    continue;
  else
    if ntz < 0
      ibz = 1;
      if steer(18)
        fprintf(fid,'\n');
        fprintf(fid,'Wave vector within Bz\n\n');
      end
    end
  end
%
% Section 2.2
%
  if ~((rk(1) == 0) & (rk(2) == 0) & (rk(3) == 0))|(ntz == 0)
    kgord = 1;
    kgel(1) = 1;
    K = 2;
    while K <= order
      rosk(1:3) = rgr(1:3,1:3,K)*rk(1:3).';
      rosk(1:3) = rosk(1:3) - rk(1:3);
      for K1 = 1:3
        if abs(rosk(K1)) <= tsmall
          rlat(K1) = 0;
        else
          if ibz
            break;
          end
          if abs(rosk(K1) - 1) <= tsmall
            rlat(K1) = 1;
          else
            if abs(rosk(K1) + 1) <= tsmall
              rlat(K1) = -1;
            else
              break;
            end
          end
        end
        if K1 == 3
          kgord = kgord + 1;
          kgel(kgord) = K;
          rkli(kgord,1:3) = rlat(1:3);
        end
      end
      K = K + 1;
    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));
%
% Wave vector equal to point group of the space group.
%
  else
    kgord = order;
    mtab2(1:kgord,1:kgord) = mtab(1:kgord,1:kgord);
    kgel(1:kgord) = 1:kgord;
    kkgel(1:kgord) = gel(1:kgord);
    kg = kgord;
    rkli = zeros(kgord,3);
  end
%
% Section 3
%
% Formation of a stable basis.
%
  fprintf(fid,'P(I).rk = rk + K(I), list K(I)\n');
  for I = 1:kgord
    fprintf(fid,'%2d (%7.4f %7.4f %7.4f)\n', I, rkli(I,1:3));
  end
  kmat(1,1:nrec) = 1:nrec;
  J = 1;
  while J <= nrec
    for I = 2:kgord
      K = kgel(I);
      for K1 = 1:3
        roklat(K1) = rgr(K1,1:3,K)*rklat(J,1:3).';
      end
      if ~ibz
        roklat(1:3) = roklat(1:3) + rkli(I,1:3);
      end
      J1 = 1;
      while J1 <= nrec
        K1 = 1;
        while K1 <= 3
          if abs(roklat(K1) - rklat(J1,K1)) > tsmall
            break;
          end
          K1 = K1 + 1;
        end
        if K1 > 3
          kmat(I,J) = J1;
          break;
        end
        J1 = J1 + 1;
      end
      if J1 > nrec
        nrec = J1;
        rklat(nrec,1:3) = roklat(1:3);
        kmat(1,nrec) = nrec;
        kmat(I,J) = nrec;
      end
    end
    J = J + 1;
  end
  fprintf(fid,'The following %3d lattice vectors form a stable basis\n', nrec);
  for I = 1:nrec
    fprintf(fid,'K(%3d) (%7.4f %7.4f %7.4f)\n', I, rklat(I,1:3));
  end
  if steer(17)
    fprintf(fid,'Stable basis transforms as follows\n');
    for I = 1:kgord
      fprintf(fid,'%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', kmat(I,1:nrec));
      fprintf(fid,'\n');
    end
  end
%
% Section 4
%
% Tests for the nonsymmorphic space group.
%
  if ~(steer(20) | ibz)
%
% Test if group of k-vector is symmorphic
%
    ksym = 0;
    I = 1;
    while I <= kgord
      I1 = kgel(I);
      K = 1;
      while K <= 3
        if abs(u(I1,K)) > tsmall
          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');
      end
    else
%
% Section 5
%
% 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.
%
      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).';
        br = T*ar;
        sres = -i*br;
        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);
          for K11 = 1:3
            tsk(K11) = bia(K11,1:3)*til(K,1:3).';
          end
          for K1 = 1:3
            rprod(K1) = rgr(K1,1:3,N3)*tsk(1:3).';
          end
          tsk(1:3) = rprod(1:3);
          for K1 = 1:3
            rprod(K1) = aib(K1,1:3)*tsk(1:3).';
            rprod(K1) = rprod(K1) + til(I,K1);
          end
          ar = rk(1:3)*rprod(1:3).';
          br = T*ar;
          R4 = -i*br;
          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
%
% Section 6
%
% Calculation of the diagonal elements of the irreducible representations.
%
%
% Section 6.1
%
  fprintf(fid,'The point group of the wave vector consists of %3d operators, indexed as nrs:\', kg);
  fprintf(fid,'\n');
  fprintf(fid,'%3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d %3d\n', kkgel(1:kg));
  fprintf(fid,'\n\n');
  if ~(steer(20) | ksym | ibz)
    fprintf(fid,'The factor group Gk/Tk consists of:\n');
    for I = 1:kgord
      fprintf(fid,'%3d point group operator: %3d non-primitive 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
  end
%
% Now we have determined the point group of the wave vector, or (for the
% nonsymmorphic case with wave vector on the BZ-boundary) the factor group
% Gk/Tk. Next we determine all (allowable) irreducible representations
% of this finite group.
%
%
% Section 6.2
%
  [jpdd,allow,ncl,laj] = Irrep(kgord,mtab2,npri,steer,ibz,ksym,nopi1,nopli1,sil);
  if steer(11)
%
% 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
    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);
%
% Section 7
%
% The formation of the projection matrices.
% Here start the loops necessary to form the projection matrix.
%
    nrec2 = 2*nrec;
%
% Section 7.1
%
    for J = 1:ncl
%
% J is the index of the irreducible representation
%
      if allow(J)
        for JD = 1:laj(J)
%
% jpdd(J,JD,I) are the diagonal elements of the irreducible representations.
%
          jpdd(J,JD,1:kgord) = conj(jpdd(J,JD,1:kgord));
          if steer(18)
            fprintf(fid,'jpdd = \n');
            jpddorg(1,1:kgord) = real(jpdd(J,JD,1:kgord));
            jpddorg(2,1:kgord) = imag(jpdd(J,JD,1:kgord));
            fprintf(fid,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n',jpddorg(1:2,1:kgord));
            fprintf(fid,'\n\n');
          end
%
% Form the projection matrix.
%
          jdpk = zeros(nrec2);
          for I = 1:kg
            IR = kkgel(I) + K96;
            for K1 = 1:nrec
              if ~(steer(20) | ksym)
                rin = 0;
                K4 = kmat(I,K1);
                if ~ibz
%
% Case 3, nonsymmorphic Gk. Wave vector on Brillouin zone boundary.
%
                  rksum(1:3) = rklat(K4,1:3) + rk(1:3);
                  rin = til(I,1:3)*rksum(1:3).';
                else
%
% Case 2, nonsymmorphic Gk. Wave vector within Brillouin zone.
%
                  IH = kgel(I);
                  rin = rklat(K4,1:3)*u(IH,1:3).';
                end
                crin = -i*rin*T;
                fac(K1) = exp(crin)*laj(J)*jpdd(J,JD,I)/kg;
              else
%
% Factors for case 2 and case 3
%
% Case 1. Symmorphic Gk.
%
                fac(K1) = laj(J)*jpdd(J,JD,I)/kg;
              end
            end
            for K5 = 1:nrec
              L = kmat(I,K5);
%
% listp(I) = I for I = 1:kg
%
              nn1 = 2*L - 1;
              nn2 = 2*K5 - 1;
              for IP = 1:2
                for IQ = 1:2
                  jdpk(IP+nn1-1,IQ+nn2-1) = jdpk(IP+nn1-1,IQ+nn2-1)...
                                          + fac(K5)*spinm(IP,IQ,IR);
                end
              end
            end
          end
          ndi = nrec2;
%
% Test projection matrix.
%
%
% Section 7.2
%
          jdprod = jdpk*jdpk;
          tmatri = jdprod - jdpk;
          K4 = 1;          
          while K4 <= ndi
            K5 = 1;
            while K5 <= ndi
              if abs(tmatri(K4,K5)) > ttsmall
                break;
              end
              K5 = K5 + 1;
            end
            if K5 <= ndi
              break;
            end
            K4 = K4 + 1;
          end
          if K4 <= ndi
            fprintf(fid,'Error, no projection matrix for J = %3d, JD = %3d, laj = %3d\n\n', J,JD,laj(J));
            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 = jdpk - jdpk';
          K4 = 1;
          while K4 <= ndi
            K5 = 1;
            while K5 <= ndi
              if abs(tmatri(K4,K5)) > ttsmall
                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 7.3
%
%
% We orthonormalise the columns of the projection matrix and store the 
% resulting 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 T-matrix.
% 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
            I = 1;
            I1 = itotal + 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
                      if itotal == ntr
                        break;
                      else
                        I = I + 1;
                        continue
                      end
                    end
                    for NC = I1: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
                        itotal = itotal - 1;
                        break;
                      else
                        tmatri(1:ndi,itotal) = tmatri(1:ndi,itotal) - R4*tmatri(1:ndi,NC);
                        LD1 = 1;
                        while abs(tmatri(LD1,itotal)) < 0.01
%
% If all elements of the column become very small after an orthogonalisation,
% we do not use the column. Renormalisation would give large errors.
%
                          if LD1 == ndi
                            itotal = itotal - 1;
                            break;
                          end
                          LD1 = LD1 + 1;
                        end
                        R1 = sum((abs(tmatri(1:ndi,itotal))).^2);
                        R1 = 1/sqrt(R1);
%
% Renormalisation.
%
                        tmatri(1:ndi,itotal) = tmatri(1:ndi,itotal)*R1;
                      end
                    end
                    if itotal == ntr
                      break;
                    end
                  end
                end
                I = I + 1;
              end
              if I > ndi
                fprintf(fid,'Error, not enough orthonormal columns\n');
                JJJ = 0;
                break;
              end
            end
          else
            itotal = 0;
          end
%
% Section 7.4
%
% Here starts the output section. The projection matrices (T-matrices) are
% printed, after that all of them have been calculated for all values of
% J, JD, ichem and L.
%
          fprintf(fid,'\n\n');
          fprintf(fid,'T-matrix for J = %3d, JD = %3d\n\n', J, JD);
          fwrite(write_proj,itotal, 'int8');
          fwrite(write_proj,ndi, 'int8');
          if itotal ~= 0
            for I = 1:itotal
              fprintf(fid,'%3d\n', I);
              I2 = 0;
              for I1 = 1:2:ndi
                I2 = I2 + 1;
                outp(1,I2) = real(tmatri(I1,I));
                outp(2,I2) = imag(tmatri(I1,I));
                outp(3,I2) = real(tmatri(I1+1,I));
                outp(4,I2) = imag(tmatri(I1+1,I));
                outp(5,I2) = I2;
              end
              fprintf(fid,'((%7.4f,%7.4f)*spin-up + (%7.4f,%7.4f)*spin-down)*K(%3d)+\n', outp(1:5,1:I2));
              fprintf(fid,'\n');
              fwrite(write_proj,outp(1:4,1:I2), 'float64');
            end
          else
            fprintf(fid,'Empty block\n\n\n');
          end
        end
        if JJJ == 0
          break;
        end
      end
    end
    if JJJ == 0
      break;
    end
  end
  last = fscanf(read_inputdata,'%2d',1);
  wvco = fscanf(read_inputdata,'%2d',1);
  for I = 1:3
    rk(I) = fscanf(read_inputdata,'%g',1);
  end
  nrec = fscanf(read_inputdata,'%2d',1);
  krep = fscanf(read_inputdata,'%2d',1);
end
%
% Section 8
%
fprintf(fid,'\n\n');
fprintf(fid,'End of SYMPWS execution');
fclose(fid);
fclose(read_spinm);
fclose(read_npgo);
fclose(read_nge);
fclose(read_NP);
fclose(read_inputdata);
fclose(read_MOh);
fclose(read_MD6h);
fclose(read_rgr3);
fclose(write_proj);
fprintf('SYMPWS has finished executing\n');