14 alpha = zeros<rowvec>(4);
30 int pIndex = p % nOrbitalsPerNuclei;
31 int qIndex = q % nOrbitalsPerNuclei;
32 int rIndex = r % nOrbitalsPerNuclei;
33 int sIndex = s % nOrbitalsPerNuclei;
35 int Rp = p / nOrbitalsPerNuclei;
36 int Rq = q / nOrbitalsPerNuclei;
37 int Rr = r / nOrbitalsPerNuclei;
38 int Rs = s / nOrbitalsPerNuclei;
41 double B = alpha[rIndex] + alpha[sIndex];
43 rowvec Ra = (alpha[pIndex]*
R.row(Rp) + alpha[qIndex]*
R.row(Rq))/A;
44 rowvec Rb = (alpha[rIndex]*
R.row(Rr) + alpha[sIndex]*
R.row(Rs))/B;
47 double t = (A*B/(A + B))*
dot(Ra-Rb,Ra-Rb);
54 return 1/sqrt(
dot(
R.row(0) -
R.row(1),
R.row(0)-
R.row(1)));
65 double f = 1.0/arg * erf(arg) *sqrt(acos(-1))/2.0;
79 int pIndex = p % nOrbitalsPerNuclei;
80 int qIndex = q % nOrbitalsPerNuclei;
82 double ap =
alpha(pIndex);
83 double aq =
alpha(qIndex);
85 double factor = ap*aq/(ap+aq);
87 int RpIndex = p / nOrbitalsPerNuclei;
88 int RqIndex = q / nOrbitalsPerNuclei;
90 rowvec Rp =
R.row(RpIndex);
91 rowvec Rq =
R.row(RqIndex);
93 double Rpq =
dot(Rp-Rq,Rp-Rq);
94 double expTerm = exp(-factor*Rpq);
95 double kin = 0.5*factor*(6-4*factor*Rpq)*pow(acos(-1)/(ap+aq),3.0/2.0)*expTerm;
108 int pIndex = p % nOrbitalsPerNuclei;
109 int qIndex = q % nOrbitalsPerNuclei;
111 double ap =
alpha(pIndex);
112 double aq =
alpha(qIndex);
114 double factor = 1.0/(ap+aq);
116 int Rp = p / nOrbitalsPerNuclei;
117 int Rq = q / nOrbitalsPerNuclei;
120 double Rpq =
dot(
R.row(Rp)-
R.row(Rq),
R.row(Rp)-
R.row(Rq));
121 double expTerm = exp(-ap*aq*factor*Rpq);
123 rowvec Rmc = (ap*
R.row(Rp) + aq*
R.row(Rq))*factor;
129 double nucAtt = -2*Z*factor*acos(-1)*expTerm*(F0p+F0q);
135 int pIndex = p % nOrbitalsPerNuclei;
136 int qIndex = q % nOrbitalsPerNuclei;
138 double ap =
alpha(pIndex);
139 double aq =
alpha(qIndex);
141 double factor = 1.0/(ap+aq);
143 int RpIndex = p / nOrbitalsPerNuclei;
144 int RqIndex = q / nOrbitalsPerNuclei;
146 rowvec Rp =
R.row(RpIndex);
147 rowvec Rq =
R.row(RqIndex);
149 double Rpq =
dot(Rp-Rq,Rp-Rq);
150 double expTerm = exp(-ap*aq*factor*Rpq);
151 double overlap = pow(acos(-1)*factor,3.0/2.0)*expTerm;