EnergyCalcsMatlab/parseAmpacFromText.m at master · clegaspi/EnergyCalcsMatlab · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
function res = parseAmpacFromText(arc, out)

    ts = textscan(arc,'%s');
    arctext = ts{1};

    ts = textscan(out,'%s');
    outtext = ts{1};

    % AMPAC changed the column labels in version 10 so it broke Dave's
    % code. Here's a workaround.
    cartoff = 7;
    netoff = 13;
    ver = find(ismember(outtext,'Version')==1);
    vernum = regexp(outtext{ver+1},'(?<v1>\d+).(?<v2>\d+).(?<v3>\d+)','names');
    if (str2double(vernum.v1) > 9)  % If higher than version 9, change the offset
        cartoff = 6;
        netoff = 11;
    end

    Hf = find(ismember(arctext,'FORMATION')==1);
    [nfound,~] = size(Hf);
    if (nfound ~= 1)
       throw(MException('parseAmpac:TextParseError',...
           ['parseAmpac: found FORMATION ',num2str(nfound),' times']));
    end
    res.Hf = str2double(arctext{Hf+2});

    % Read in cartesian coordinates
    t1 = find(ismember(outtext,'CARTESIAN')==1);
    % take last occurence of the word cartesian
    t1 = t1(size(t1,1),1);
    t1 = t1 + cartoff;
    done = 0;
    iatom = 1;
    while (~done)
       if (str2num(outtext{t1}) == iatom)
          element{iatom} = outtext{t1+1};
          r(1,iatom) = str2double( outtext{t1+2} );
          r(2,iatom) = str2double( outtext{t1+3} );
          r(3,iatom) = str2double( outtext{t1+4} );
          t1 = t1 + 5;
          iatom = iatom + 1;
       else
          done = 1;
          iatom = iatom -1;
       end
    end
    res.natom = iatom;
    res.r = r;
    res.element = element;

    %%
    t1 = find(ismember(outtext,'NET')==1);
    if (size(t1,1) > 0)
       % take last occurence of the word 'NET'
       t1 = t1(size(t1,1),1);
       t1 = t1 + netoff;
       done = 0;
       iatom = 1;
       while (~done)
          if (str2num(outtext{t1}) == iatom)
             c_element{iatom} = outtext{t1+1};
             charge(iatom) = str2double( outtext{t1+2} );
             t1 = t1 + 4;
             iatom = iatom + 1;
          else
             done = 1;
             iatom = iatom -1;
          end
       end
       res.rho = charge;
    end

    %% Read in charges
    t1 = find(ismember(outtext,'ELECTROSTATIC')==1);
    if (size(t1,1) > 1)
       % take last occurence of the word electrostatic
       t1 = t1(size(t1,1),1);
       t1 = t1 + 7;     % This may be version dependent...if it breaks, you know why
       totalCharge = str2double( outtext{t1} );
       t1 = t1 + 7;
       done = 0;
       iatom = 1;
       while (~done)
          if (str2num(outtext{t1}) == iatom)
             c_element{iatom} = outtext{t1+1};
             charge(iatom) = str2double( outtext{t1+2} );
             t1 = t1 + 3;
             iatom = iatom + 1;
          else
             done = 1;
             iatom = iatom -1;
          end
       end
       res.esp = charge;
    end
end