-
Notifications
You must be signed in to change notification settings - Fork 1
/
mocalc.m
210 lines (179 loc) · 8.14 KB
/
mocalc.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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
% Contains the SCF algorithm.
%
% Input:
% atoms list of element numbers (array with K elements);
% e.g. [6 8] for CO
% xyz_a0 K×3 array of Cartesian coordinates of nuclei, in bohr
% totalcharge total charge of the molecule, in units of elementary charge
% settings a structure that contains several fields:
% .basisset string specifying the basis set, e.g. '6-31G', 'cc-pVDZ', 'STO-3G'
% .tolEnergy SCF convergence tolerance for the energy (hartrees)
% .tolDensity SCF convergence tolerance for the density (a0^3)
% .method string specifying the method, either 'RHF' for restricted Hartree-Fock
% or 'RKS' for restricted Kohn-Sham DFT
% .ExchFunctional string specifying the exchange functional, 'Slater'
% .CorrFunctional string specifying the correlation functional, 'VWN3' or 'VWN5'
% .nRadialPoints number of radial points for the integration grid
% .nAngularPoints number of angular points for the integration grid
%
% Output:
% out a structure that contains several fields:
% .basis list of basis functions, as generated by buildbasis
% .S overlap matrix (M×M)
% .T kinetic energy matrix (M×M)
% .Vne electron-nuclear attraction matrix (M×M)
% .J matrix of Coulomb integrals (M×M)
% .K matrix of exchange integrals (M×M)
% .ERI 4D array of electron-electron repulsion integrals (M×M×M×M)
% .epsilon MO energies (1×M), in hartrees, in ascending order,
% consisting of both occupied and virtual orbitals
% .C MO coefficient matrix (M×M), of occupied and virtual
% orbitals, sorted in ascending order of orbital energy
% .P density matrix (M×M)
% .E0 electronic ground-state energy of the molecule, in hartrees
% .Etot total ground-state energy (including nuclear-nuclear repulsion;
% but without the vibrational zero-point energy), in hartrees
% .Exc exchange-correlation energy, in hartrees
% .Vxc matrix of exchange-correlation integrals (MxM), in hartrees
% .rhoInt integral of electron density over all 3D space
function out = mocalc(atoms,xyz_a0,totalcharge,settings)
% Cleaning the inputs.
BasisSetName = settings.basisset;
E_tol = settings.tolEnergy; % tolerance energy change (hartrees)
P_tol = settings.tolDensity; % tolerance density element change ((a0)^-3)
method = settings.method;
% Defining some variables that will be used throughout.
out = {}; % initializing the return structure.
N = sum(atoms) - totalcharge; % total # of e in molecule.
% Building the outputs.
bdef = basisread(BasisSetName);
basis = buildbasis(atoms,xyz_a0,bdef);
S = int_overlap(basis); % symmetrized, diagonals = 1
T = int_kinenergy(basis);
Vne = int_attraction(atoms,xyz_a0,basis);
Vnn = nucnucrepulsion(atoms,xyz_a0);
ERI = int_repulsion(basis);
M = numel(basis); % number of basis functions.
out.basis = basis;
out.S = S;
out.T = T;
out.Vne = Vne;
out.ERI = ERI;
% SCF Loop.
counter = 0;
converged = false;
while ~converged
counter = counter + 1;
if method == 'RHF'
if counter == 1
F = T + Vne; % initial approx F matrix (ignore e-e rep).
else
F = T + Vne + J - K; % MxM F matrix for SCF iterations after
% starting density matrix has been estimated.
end
[C,epsi] = eig(F,S); % Matlab's general eigenproblem solver to
% approx MO coeff and energy matrices.
diagonals = diag(epsi); % diagonal energy matrix -> column vect.
[sorted,idx] = sort(diagonals);
epsi = sorted; % Mx1 matrix of energy vals in ascending order.
for k = 1:M
norms(k) = sqrt(C(:,k).'*S*C(:,k));
C(:,k) = C(:,k)./norms(k); % normalizing C
end
C = C(:,idx); % C sorted based on sorted e values.
C_occ = C(:,1:N/2); % reducing C to occupied orbitals only
P = 2*(C_occ*C_occ.'); % estimation of starting density matrix.
J = zeros(M); % initializing Coulomb matrix.
K = zeros(M); % initializing Exchange matrix.
for mu = 1:M
for nu = 1:M
for kap = 1:M
for lam = 1:M
J(mu,nu) = J(mu,nu) + P(kap,lam)*ERI(mu,nu,lam,kap); % calc J elements
K(mu,nu) = K(mu,nu) + 0.5*P(kap,lam)*ERI(mu,kap,lam,nu); % calc K elements
end
end
end
end
E0 = 0;
for u = 1:M
for v = 1:M
E0 = E0 + P(u,v)*(T(u,v) + Vne(u,v) + 0.5*(J(u,v) - K(u,v)));
end
end
if counter > 1
E_change = E0 - E0_prev;
P_change = max(abs((P(:) - P_prev(:))));
if E_change < E_tol && P_change < P_tol
converged = true;
end
end
E0_prev = E0;
P_prev = P;
elseif method == 'RKS'
% Adding field to settings structure for DFT;RKS method.
ExchFunctional = settings.ExchFunctional;
CorrFunctional = settings.CorrFunctional;
nRadialPoints = settings.nRadialPoints;
nAngularPoints = settings.nAngularPoints;
grid = molecular_grid(atoms,xyz_a0,nRadialPoints,nAngularPoints);
K = 0;
if counter == 1
F = T + Vne; % initial approx F matrix (ignore e-e rep).
else
F = T + Vne + J + Vxc; % MxM F matrix for SCF iterations after
% starting density matrix has been estimated.
end
[C,epsi] = eig(F,S); % Matlab's general eigenproblem solver to
% approx MO coeff and energy matrices.
diagonals = diag(epsi); % diagonal energy matrix -> column vect.
[sorted,idx] = sort(diagonals);
epsi = sorted; % Mx1 matrix of energy vals in ascending order.
for k = 1:M
norms(k) = sqrt(C(:,k).'*S*C(:,k));
C(:,k) = C(:,k)./norms(k); % normalizing C
end
C = C(:,idx); % C sorted based on sorted e values.
C_occ = C(:,1:N/2); % reducing C to occupied orbitals only
P = 2*(C_occ*C_occ.'); % estimation of starting density matrix.
J = zeros(M); % initializing Coulomb matrix.
% K = zeros(M); % initializing Exchange matrix.
for mu = 1:M
for nu = 1:M
for kap = 1:M
for lam = 1:M
J(mu,nu) = J(mu,nu) + P(kap,lam)*ERI(mu,nu,lam,kap); % calc J elements
% K(mu,nu) = K(mu,nu) + P(kap,lam)*ERI(mu,kap,lam,nu); % calc K elements
end
end
end
end
[Vxc,Exc,rhoInt] = int_xc(basis,P,grid,ExchFunctional,CorrFunctional);
E0 = Exc;
for u = 1:M
for v = 1:M
E0 = E0 + P(u,v)*(T(u,v) + Vne(u,v) + 0.5*J(u,v));
end
end
if counter > 1
E_change = E0 - E0_prev;
P_change = max(abs((P(:) - P_prev(:))));
if E_change < E_tol && P_change < P_tol
converged = true;
end
end
E0_prev = E0;
P_prev = P;
out.Exc = Exc;
out.Vxc = Vxc;
out.rhoInt = rhoInt;
end
end
out.J = J;
out.K = K;
out.epsilon = epsi;
out.C = C;
out.P = P;
out.E0 = E0;
out.Etot = E0 + Vnn;
end