function MatrixU()
% forms Uq or Us

% define foldername
modelfolder = getmodel;

% read model files
InfTable = zeros(1,9);
filename = strcat(modelfolder,'InfTable.txt');
fid = fopen(filename,'r');
InfTable = fscanf(fid,'%i',8);
fclose(fid);
nVibrCoords = InfTable(7);
nBlocks = InfTable(8);

Uq = zeros(nVibrCoords,nVibrCoords);
filename = strcat(modelfolder,'UqSparse.txt');
load(filename);
mu = size(UqSparse, 1);
for i   = 1:mu
    iu = UqSparse(i,1);
    ju = UqSparse(i,2);
    u = UqSparse(i,3);
    Uq(iu, ju) = u;
    Uq(ju, iu) = u;
end % i
clear UqSparse iu ju u mu

for i = 1:nVibrCoords
    if Uq(i,i) == 0
        s = sprintf('i=%i Uq(i,i)=0 ',i);
        disp(s)
        return
    end    
end % i

if nBlocks == 1
    filename = strcat(modelfolder,'Uq.txt');
    fid = fopen(filename,'w');
    nUq3 = nVibrCoords*(nVibrCoords+1)/2; % elements in Uq upper triangle
    Uq3 = zeros(1,nUq3); % Uq upper triangle
    k = 0;
    for i = 1:nVibrCoords
        for j = i:nVibrCoords
            k = k + 1;
            Uq3(k) = Uq(i,j);
        end % j    
    end % i
    fprintf(fid,'  %e  %e  %e  %e \n',Uq3);
    fclose(fid);
    return
end % nBlocks == 1

% Then form Us
C = zeros(nVibrCoords,nVibrCoords);
nC = zeros(1,nBlocks);
filename = strcat(modelfolder,'Cmatrix.txt');
fid = fopen(filename,'rt');
s = fgetl(fid);
nC = fscanf(fid,'%i',nBlocks);
irow = 0;
for iblock = 1:nBlocks
    nrows = nC(iblock);
    for i = 1:nrows
        irow = irow + 1;
        ncols = fscanf(fid,'%i',1);
        for jcol = 1:ncols
            j = fscanf(fid,'%i',1);
            c = fscanf(fid,'%f',1);
            C(irow, j) = c;
        end % jcol
    end % i
end % iblock
fclose(fid);
Us = C*Uq*C';

filename = strcat(modelfolder,'Us.txt');
fid = fopen(filename,'w');
i1 = 0;
i2 = 0;
for iblock = 1:nBlocks
    nUs3 = nC(iblock)*(nC(iblock)+1)/2; % elements in Us upper triangle
    Us3 = zeros(1,nUs3); % Us upper triangle
    k = 0;
    i1 = i2 + 1;
    i2 = i2 + nC(iblock);
    for i = i1:i2
        for j= i:nVibrCoords
            if j <= i2
                k = k + 1;
                Us3(k) = Us(i,j);
            else
                usij = Us(i,j);
                if abs(usij) > 0.0001
                    s = sprintf('i=%i j=%i Us(i,j)=%f',i,j,usij);
                    disp(s)
                    return
               end % if
            end % if
        end % j
    end % i
    fprintf(fid,'  %e  %e  %e  %e \n',Us3);
    fprintf(fid,' \n');
end % iblock
fclose(fid);

% subfunctions ---------------------
function modelfolder = getmodel()
% gets model name in LevInfinite

filename = 'C:\LevInfinite\Config\Molecules.ini';
fid = fopen(filename,'rt');
molfolder = fgetl(fid);
fclose(fid);

filename = 'C:\LevInfinite\Config\ToCalculate.ini';
fid = fopen(filename,'rt');
modelname = fgetl(fid);
fclose(fid);

modelfolder = strcat(molfolder,modelname);
modelfolder = strcat(modelfolder,'\');

