-
Notifications
You must be signed in to change notification settings - Fork 1
/
basisread.m
94 lines (79 loc) · 2.93 KB
/
basisread.m
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
% Reads in a basis set definition file (in Gaussian94 format)
% and returns a structure containing all basis set details.
%
% Input:
% BasisSetName basis set name, e.g. '6-31G'
%
% Output:
% BasisSet cell array with one entry per element
% BasisSet{iElement} is a structure array containing the fields
% .shelltype type of shell ('S', 'P', 'SP', 'D', etc)
% .coeffs contraction coefficients
% .exponents radial exponents of primitives
%
% The basis set files are expected to be in a subfolder called
% basissets. The file names are expected to have the extension
% .basis. For instance, the 6-31G basis should be in the file
% basissets/6-31G.basis
function BasisSet = basisread(BasisSetName)
basissetfilename = ['basissets/' BasisSetName '.basis'];
% Read all lines from file into a single cell array (one cell per line)
f = fopen(basissetfilename);
if f<0
error('Could not open file ''%s''.',basissetfilename);
end
allLines = textscan(f,'%s','Delimiter','\n','Whitespace','');
allLines = allLines{1};
fclose(f);
% Remove comment lines (starting with !) and empty lines
for iLine = 1:numel(allLines)
rmv(iLine) = isempty(allLines{iLine}) || (allLines{iLine}(1)=='!');
end
allLines(rmv) = [];
% Identify lines delimiting atoms (****)
for iLine = 1:numel(allLines)
delimLines(iLine) = (allLines{iLine}(1)=='*');
end
delimLines = find(delimLines);
nEntries = numel(delimLines)-1; % number of entries (i.e. atoms) in the file
for iEntry = 1:nEntries
% Extract lines for this atom
lines = allLines(delimLines(iEntry)+1:delimLines(iEntry+1)-1);
% Get element number
ElementNo = elementsymb2no(lines{1}(1:2));
% Determine number of shells
shellDelim = 0;
for iLine = 2:numel(lines)
shellDelim(iLine) = lines{iLine}(1)~=' ';
end
shellDelim(end+1) = true;
shellDelim = find(shellDelim);
nShells = numel(shellDelim)-1;
% Read in each shell in turn
for iShell = 1:nShells
% Get shell type (S, P, SP, D, etc) from shell header line
shellHeader = lines{shellDelim(iShell)};
newShell.shelltype = strtrim(shellHeader(1:2));
% Get exponents and contraction coefficients
stri = char(lines(shellDelim(iShell)+1:shellDelim(iShell+1)-1));
numbers = str2num(stri); %#ok
exponents = numbers(:,1).';
coeffs = numbers(:,2:end).';
% Store exponents and contraction coefficients in output structure
newShell.exponents = exponents;
newShell.coeffs = coeffs;
BasisSet{ElementNo}(iShell) = newShell;
end
end
function elno = elementsymb2no(symbol)
if numel(symbol)==1, symbol(2) = ' '; end
ElementSymbols = [...
'H He'...
'LiBeB C N O F Ne'...
'NaMgAlSiP S ClAr'...
'K CaScTiV CrMnFeCoNiCuZnGaGeAsSeBrKr'...
'RbSrY ZrNbMoTcRuRhPdAgCdInSnSbTeI Xe'...
'CsBaLaCePrNdPmSmEuGdTbDyHoerTmYbLuHfTaW ReOsIrPtAuHgTlPbBiPoAtRn'...
'FrRaAcThPaU NpPuAmCmBkCfEsFmMdNoLrRfDbSgBhHsMtDsRgCnNhFlMcLvTsOg'];
idx = strfind(ElementSymbols,symbol);
elno = (idx+1)/2;