PVOX-readOuTput.m函数

function [data, wfn, D_vlid, W_valid] = readOutput(outName, datName)
%
%READOUTPUT
%
%       Runs parsec parser and gathers data
%

D_valid = false;
W_valid = [0,0,0,0,0];
data = [];
wfn = [];

% Open parsec.out
try
    h = waitbar(0,'Opening parsec.out');

    [str,maxsize,endian] = computer;

    fName = fopen('pTemp/fname.dat','w');

    fwrite(fName,outName);

    fclose(fName);
    waitbar(0.1,h);

    if strncmp(str,'PCWIN',5)
        try   % Parsec Abridged Output
            [s,w] = dos('del /Q pTemp\main.dat');
            !"./pTools/PARSEC_PARSER_ABRIDGED.py"
        catch % Parsec Un-Abridged Output
            !"./pTools/PARSEC_PARSER.py"
        end     
    else
        [s,w] = unix('rm -f pTemp/main.dat');
  try
    [s,w] = unix('python ./pTools/PARSEC_PARSER_ABRIGED.py');
  catch
    [s,w]=unix('python ./pTools/PARSEC_PARSER.py');
    if (s ~= 0)
      fclose('all');
      close(h);
      return;
    end
    data.abridged = 0;
   end
  if (s ~= 0)
    [s,w]=unix('python ./pTools/PARSEC_PARSER.py');
    if (s ~= 0)
      fclose('all');
                close(h);
                return;
            end
            data.abridged = 0;
        else
            data.abridged = 1;
        end
        [s,w] = unix('chmod 660 pTemp/*');
    end
    waitbar(0.2,h);

    inMain = load('pTemp/main.dat');
    data.timeDiag = inMain(:,1);
    data.FermiLevel = inMain(:,2);
    data.timeHart = inMain(:,3);
    data.ChargeDmin = inMain(:,4);
    data.ChargeDmax = inMain(:,5);
    data.Eig_Energy = inMain(:,6);
    data.H_Energy = inMain(:,7);
    data.Int_Vxc_rho = inMain(:,8);
    data.Exc = inMain(:,9);
    data.E_I_Energy = inMain(:,10);
    data.I_I_Energy = inMain(:,11);
    data.Enew_Eold = inMain(:,12);
    data.Tot_Energy = inMain(:,13);
    data.Energy_atom = inMain(:,14);
    data.SREpot = inMain(:,15);
    data.ChrgWghtPot = inMain(:,16);

    fOptions = fopen('pTemp/misc.dat');
    opt = fgetl(fOptions);
    data.num_states = str2num(strtok(opt));
    opt = fgetl(fOptions);
    data.num_atoms = str2num(strtok(opt));

    iter = size(data.timeDiag);
    data.num_iterations = iter(1);

    waitbar(0.3,h);

    inEigs = load('pTemp/eigenvalues.dat');
    m = data.num_states;
    n = size(inEigs(:,1));
    n = n(1);
    n = n/m;
    data.State = reshape(inEigs(:,1),m,n);
    data.EigenvalueRY = reshape(inEigs(:,2),m,n);
    data.EigenvalueEV = reshape(inEigs(:,3),m,n);
    data.Occup = reshape(inEigs(:,4),m,n);
    data.Repr = reshape(inEigs(:,5),m,n);

    inForces = load('pTemp/forces.dat');
    forceIon = reshape(inForces(1,:),3,data.num_atoms);
    data.forceIon = transpose(forceIon);
    forceSelfC = reshape(inForces(2,:),3,data.num_atoms);
    data.forceSelfC = transpose(forceSelfC);
    forceCoreC = reshape(inForces(3,:),3,data.num_atoms);
    data.forceCoreC = transpose(forceCoreC);

    inForces2 = load('pTemp/forces2.dat');
    forceLocal = reshape(inForces2(1,:),6,data.num_atoms);
    data.forceLocal = transpose(forceLocal);
    forceNonLcl = reshape(inForces2(2,:),6,data.num_atoms);
    data.forceNonLcl = transpose(forceNonLcl);
    forceTotal = reshape(inForces2(3,:),6,data.num_atoms);
    data.forceTotal = transpose(forceTotal);
    forceTotal2 = reshape(inForces2(4,:),6,data.num_atoms);
    data.forceTotal2 = transpose(forceTotal2);

    inForces3 = load('pTemp/forces3.dat');
    data.dipole = reshape(inForces3(1,:),1,3);
    data.netForce = reshape(inForces3(2,:),1,3);
    data.netTorque = reshape(inForces3(3,:),1,3);

    D_valid = true;
end

waitbar(0.3,h);

% OPEN WFN.DAT
try
    waitbar(0.4,h,'Opening wfn.dat   ');

    [wfn, pot, rho, wavedata, XYZ] = wfnread(datName, h);

    % wfn.dat read to completion
    W_valid(1) = 1;

    waitbar(0.6,h,'Processing wfn.dat');

    p = floor( (2*wfn.radius-1)/wfn.step + 1 );

    o = floor(wfn.radius/wfn.step)+1;

    wfn.pot = zeros( p,p,p );
    wfn.rho = zeros( p,p,p );

    % Determine degree of symmetry
    sym = 1;
    sym = round (wfn.ndim / wfn.num_wedges);

    % Initialize sparse matracies
    wfn.pot = zeros( p,p,p );
    wfn.rho = zeros( p,p,p );
    wfn.wfn = zeros(p,p,p,wfn.num_states);

    % Read in data and reflect if necessary
    if (sym == 1)
        % Known symmetry found
        W_valid(2) = 1;

        % Read pot and rho
        for i=(1:wfn.num_wedges)
            wfn.pot( XYZ(i,2)+o, XYZ(i,1)+o,XYZ(i,3)+o ) = pot(i);
            wfn.rho( XYZ(i,2)+o, XYZ(i,1)+o,XYZ(i,3)+o ) = abs(rho(i));
        end
        % pot and rho read correctly
        W_valid(3) = 1;
        waitbar(0.65,h);

        % Read wfn's
        for i = (1:wfn.num_states)
            for j = (1:wfn.num_wedges)
                wfn.wfn( XYZ(j,2)+o, XYZ(j,1)+o,XYZ(j,3)+o, i) = abs(wavedata(j,i));
            end
        end
        % wfn's read correctly
        W_valid(4) = 1;

    elseif (sym == 4)
        % Check if data is in the correct regions
        a = min(XYZ(:,1));
        b = min(XYZ(:,2));
        if (a+b == 0)
            % Known symmetry found
            W_valid(2) = 1;

            % Read and reflect pot and rho
            for i=(1:wfn.num_wedges)
                wfn.pot( XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot( XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);

                wfn.rho( XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho( XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
            end
            % pot and rho read correctly
            W_valid(3) = 1;       
            waitbar(0.65,h);

            % Read and reflect wfn's
            for i = (1:wfn.num_states)
                for j = (1:wfn.num_wedges)
                    wfn.wfn( XYZ(j,2)+o, XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn( XYZ(j,2)+o,-XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o, XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o,-XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                end
            end
            % wfn's read correctly
            W_valid(4) = 1;
        end

    elseif (sym == 8)
        % Check if data is in the correct regions
        a = min(XYZ(:,1));
        b = min(XYZ(:,2));
        c = min(XYZ(:,3));
        if (a+b+c == 0)
            % Known symmetry found
            W_valid(2) = 1;

            % Read and reflect pot and rho
            for i=(1:wfn.num_wedges)
                wfn.pot( XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot( XYZ(i,2)+o, XYZ(i,1)+o,-XYZ(i,3)+o ) = pot(i);
                wfn.pot( XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot( XYZ(i,2)+o,-XYZ(i,1)+o,-XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o, XYZ(i,1)+o,-XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = pot(i);
                wfn.pot(-XYZ(i,2)+o,-XYZ(i,1)+o,-XYZ(i,3)+o ) = pot(i);

                wfn.rho( XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho( XYZ(i,2)+o, XYZ(i,1)+o,-XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho( XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho( XYZ(i,2)+o,-XYZ(i,1)+o,-XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o, XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o, XYZ(i,1)+o,-XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o,-XYZ(i,1)+o, XYZ(i,3)+o ) = abs(rho(i));
                wfn.rho(-XYZ(i,2)+o,-XYZ(i,1)+o,-XYZ(i,3)+o ) = abs(rho(i));
            end
            % pot and rho read correctly
            W_valid(3) = 1;
            waitbar(0.65,h);

            % Read and reflect wfn's
            for i = (1:wfn.num_states)
                for j = (1:wfn.num_wedges)
                    wfn.wfn( XYZ(j,2)+o, XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn( XYZ(j,2)+o, XYZ(j,1)+o,-XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn( XYZ(j,2)+o,-XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn( XYZ(j,2)+o,-XYZ(j,1)+o,-XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o, XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o, XYZ(j,1)+o,-XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o,-XYZ(j,1)+o, XYZ(j,3)+o, i) = abs(wavedata(j,i));
                    wfn.wfn(-XYZ(j,2)+o,-XYZ(j,1)+o,-XYZ(j,3)+o, i) = abs(wavedata(j,i));
                end
            end
            % wfn's read correctly
            W_valid(4) = 1;
        end
    end

    wfn.wfn_min = min(min(min(min(wfn.wfn))));
    wfn.wfn_max = max(max(max(max(wfn.wfn))));
    waitbar(0.85,h);

    % Compute number of z slices
    [trash, trash, wfn.n, trash] = size(wfn.wfn);

    q = ( -1*wfn.radius:wfn.step: (wfn.radius-wfn.step) );
    [wfn.X,wfn.Y,wfn.Z] = meshgrid( (wfn.shift(1)*wfn.step+q), (wfn.shift(2)*wfn.step+q), (wfn.shift(3)*wfn.step+q) );
    % meshgrid created correctly
    W_valid(5) = 1;
end

% Read in Atom Data
try
    waitbar(0.9,h,'Reading Atom Data ');

    % atom data
    data.numTypes = 0;
    i = 1;
    fAtom = fopen('pTemp/atom.dat');
    input = fgetl(fAtom);
    count = 0;
    while input ~= -1
        data.numTypes = data.numTypes + 1;
        inputd = sscanf(mat2str(input), '%d %s');

        % get the name of the atom
        if ( size(inputd) < 3 )
            b = char(inputd(2));
            b(2) = ' ';
        else
            b = char(inputd(2:size(inputd)));
        end

        for ( j = 1:inputd(1) )
            data.atomName(count+j,:) = b;
        end

        data.typeNum(data.numTypes) = inputd(1);
        data.typeName(data.numTypes,:) = b;
        data.typeSize(data.numTypes) = 1;
        switch(b)
            case 'C '
                data.typeColor(data.numTypes,:) = [.3 .3 .3];
            case 'O '
                data.typeColor(data.numTypes,:) = [.938 0 0];
            case 'N '
                data.typeColor(data.numTypes,:) = [.368 .705 1];
            case 'H '
                data.typeColor(data.numTypes,:) = [1 1 1];
            case 'S '
                data.typeColor(data.numTypes,:) = [1 1 0];
            case 'P '
                data.typeColor(data.numTypes,:) = [1 .446 0];
            case 'Cl'
                data.typeColor(data.numTypes,:) = [0 1 0];
            case 'Ca'
                data.typeColor(data.numTypes,:) = [.4 .4 .4];
            otherwise
                data.typeColor(data.numTypes,:) = [.491 0 .385];
        end
        count = count + inputd(1);
        input = fgetl(fAtom);
    end
catch
    D_valid = false;
end

fclose('all');

waitbar(1,h);

close(h);