/
mocalc.m
148 lines (127 loc) · 5.48 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
%% milestone 5 part (d)
%% (also milestone 4 parts (c)-(e))
% defines function mocalc
% documentation (from stefan's notes):
%
% out = mocalc(atoms,xyz_a0,totalcharge,settings)
%
% 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 (a_0^(-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, 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)
% calculates the Coulumb and exchange matrices J and K using the given values for P and the electron-electron repulsion integrals ERI
function [J K] = calculateJK(P, ERI, nBasis)
J = zeros(nBasis, nBasis);
K = zeros(nBasis, nBasis);
for mu = 1:nBasis
for nu = 1:nBasis
for kappa = 1:nBasis
for lambda = 1:nBasis
J(mu,nu) = J(mu,nu) + P(kappa, lambda)*ERI(mu, nu, lambda, kappa);
K(mu,nu) = K(mu,nu) + (1/2)*P(kappa, lambda)*ERI(mu, kappa, lambda, nu);
end
end
end
end
end
% build the basis and record its length
basis = buildbasis(atoms, xyz_a0, basisread(settings.basisset));
nBasis = length(basis);
% N = number of unique spatial orbitals required for RHF
N = (sum(atoms) - totalcharge)/2;
% compute quantities that are fixed (do not get updated in the SCF loop below)
S = int_overlap(basis);
T = int_kinenergy(basis);
Vne = int_attraction(atoms, xyz_a0, basis);
ERI = int_repulsion(basis);
% initialize P (and epsilon) to all zeros
% this way the first iteration of the loop below drops two-electron terms from the Fock matrix, as required
P = zeros(nBasis, nBasis);
epsilon = zeros(nBasis,nBasis);
if settings.method == 'RKS'
grid = molecular_grid(atoms, xyz_a0, settings.nRadialPoints, settings.nAngularPoints);
end
% this is the SCF loop; it iterates until the convergence criteria are met, whence the "break" statement below is executed to terminate the loop
while true
% calculate J and K using the current density matrix P and the (fixed) values of ERI and nBasis
[J K] = calculateJK(P, ERI, nBasis);
% if we are doing DFT, compute exchange-correlation quantities and calculate F appropriately
if settings.method == 'RKS'
[Vxc Exc rhoInt] = int_xc(basis, P, grid, settings.ExchFunctional, settings.CorrFunctional);
F = T + Vne + J + Vxc;
% otherwise we are doing HF; calculate F using J and K
else
F = T + Vne + J - K;
end
% solve the Roothan-Hall equations for the current F
[C epsilonNew] = eig(F,S);
% ensure C and epsilonNew are sorted by increasing energy
[eps I] = sort(diag(epsilonNew));
epsilonNew = epsilonNew(I,I);
C = C(:,I);
% ensure C is properly normalized
normalizingConstants = diag(sqrt(C'*S*C));
for k = 1:length(basis)
C(:,k) = C(:,k)/normalizingConstants(k);
end
% calculate the new density matrix PNew
PNew = 2*C(:,1:N)*C(:,1:N)';
% check if we have met the convergence criteria
if max(max(abs(P - PNew))) < settings.tolDensity && max(max(abs(epsilon - epsilonNew))) < settings.tolEnergy
% if we have, break out of the update loop
break
end
% if not, store the new values of P and epsilon and run another iteration
P = PNew;
epsilon = epsilonNew;
end
% populate the fields of the output structure
out.basis = basis;
out.S = S;
out.T = T;
out.Vne = Vne;
out.ERI = ERI;
out.epsilon = diag(epsilon);
out.C = C;
out.P = P;
out.J = J;
out.K = K;
if settings.method == 'RKS'
out.E0 = sum(sum(transpose(P).*(T + Vne + J/2))) + Exc;
out.K = zeros(nBasis, nBasis);
out.Vxc = Vxc;
out.Exc = Exc;
out.rhoInt = rhoInt;
else
out.E0 = sum(sum(transpose(P).*(T + Vne + (J - K)/2)));
end
out.Etot = out.E0 + nucnucrepulsion(atoms, xyz_a0);
end