実装
本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。
最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……
#include <math.h>
#include <stdio.h>
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);
double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
V1111, V2111, V2121, V2211, V2221, V2222;
double COEF[3][3] = {{1.0, 0.678914, 0.444635},
{0.0, 0.430129, 0.535328},
{0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
{0.0, 0.851819, 0.405771},
{0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;
double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];
int main(void) {
int IOP = 2;
int N = 3;
double R = 1.4632;
double ZETA1 = 2.0925;
double ZETA2 = 1.24;
double ZA = 2.0;
double ZB = 1.0;
HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
double ZB) {
if (IOP != 0) {
printf(" STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
}
INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
double ZB) {
double R2 = R * R;
for (int i = 0; i < N; ++i) {
A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
}
S12 = 0.0;
T11 = 0.0;
T12 = 0.0;
T22 = 0.0;
V11A = 0.0;
V12A = 0.0;
V22A = 0.0;
V11B = 0.0;
V12B = 0.0;
V22B = 0.0;
V1111 = 0.0;
V2111 = 0.0;
V2121 = 0.0;
V2211 = 0.0;
V2221 = 0.0;
V2222 = 0.0;
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
double RAP = A2[j] * R / (A1[i] + A2[j]);
double RAP2 = pow(RAP, 2.0);
double RBP2 = pow(R - RAP, 2.0);
S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
}
}
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
for (int k = 0; k < N; ++k) {
for (int l = 0; l < N; ++l) {
double RAP = A2[i] * R / (A2[i] + A1[j]);
double RAQ = A2[k] * R / (A2[k] + A1[l]);
double RBQ = R - RAQ;
double RPQ = RAP - RAQ;
double RAP2 = RAP * RAP;
double RBQ2 = RBQ * RBQ;
double RPQ2 = RPQ * RPQ;
V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0)
* D1[i] * D1[j] * D1[k] * D1[l];
V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
* D2[i] * D1[j] * D1[k] * D1[l];
V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
* D2[i] * D1[j] * D2[k] * D1[l];
V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
* D2[i] * D2[j] * D1[k] * D1[l];
V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
* D2[i] * D2[j] * D2[k] * D1[l];
V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
D2[j] * D2[k] * D2[l];
}
}
}
}
if (IOP != 0) {
printf(" R ZETA1 ZETA2 S12 T11\n\n");
printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
printf(" T12 T22 V11A V12A V22A\n\n");
printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
printf(" V11B V12B V22B V1111 V2111\n\n");
printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
V2111);
printf(" V2121 V2211 V2221 V2222\n\n");
printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
}
}
double calcS(double A, double B, double RAB2) {
return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}
double calcT(double A, double B, double RAB2) {
return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}
double F0(double ARG) {
double ret;
if (ARG < 1.0e-6) {
ret = 1.0 - ARG / 3.0;
} else {
ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
}
return ret;
}
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
return -v * ZC;
}
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
double RPQ2) {
return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
* F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
* exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
double ZB) {
H[0][0] = T11 + V11A + V11B;
H[0][1] = T12 + V12A + V12B;
H[1][0] = H[0][1];
H[1][1] = T22 + V22A + V22B;
S[0][0] = 1.0;
S[0][1] = S12;
S[1][0] = S[0][1];
S[1][1] = 1.0;
X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
X[1][0] = X[0][0];
X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
X[1][1] = -X[0][1];
XT[0][0] = X[0][0];
XT[0][1] = X[1][0];
XT[1][0] = X[0][1];
XT[1][1] = X[1][1];
TT[0][0][0][0] = V1111;
TT[1][0][0][0] = V2111;
TT[0][1][0][0] = V2111;
TT[0][0][1][0] = V2111;
TT[0][0][0][1] = V2111;
TT[1][0][1][0] = V2121;
TT[0][1][1][0] = V2121;
TT[1][0][0][1] = V2121;
TT[0][1][0][1] = V2121;
TT[1][1][0][0] = V2211;
TT[0][0][1][1] = V2211;
TT[1][1][1][0] = V2221;
TT[1][1][0][1] = V2221;
TT[1][0][1][1] = V2221;
TT[0][1][1][1] = V2221;
TT[1][1][1][1] = V2222;
if (IOP != 0) {
MATOUT(S, 2, 2, 2, 2, "S ");
MATOUT(X, 2, 2, 2, 2, "X ");
MATOUT(H, 2, 2, 2, 2, "H ");
printf("\n");
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
printf(" (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
}
}
}
}
}
}
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
double ZB) {
double CRIT = 1.0e-4;
int MAXIT = 25;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
P[i][j] = 0.0;
}
}
if (IOP == 2) {
MATOUT(P, 2, 2, 2, 2, "P ");
}
int ITER = 0;
double EN;
for (;;) {
ITER++;
if (IOP == 2) {
printf("\n START OF ITERATION NUMBER = %2d\n", ITER);
}
FORMG();
if (IOP == 2) {
MATOUT(G, 2, 2, 2, 2, "G ");
}
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
F[i][j] = H[i][j] + G[i][j];
}
}
EN = 0.0;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
}
}
if (IOP == 2) {
MATOUT(F, 2, 2, 2, 2, "F ");
printf("\n\n\n ELECTRONIC ENERGY = %20.12lf\n", EN);
}
MULT(F, X, G, 2, 2);
MULT(XT, G, FPRIME, 2, 2);
DIAG(FPRIME, CPRIME, E);
MULT(X, CPRIME, C, 2, 2);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
OLDP[i][j] = P[i][j];
P[i][j] = 0.0;
for (int k = 0; k < 1; ++k) {
P[i][j] += 2.0 * C[i][k] * C[j][k];
}
}
}
if (IOP == 2) {
MATOUT(FPRIME, 2, 2, 2, 2, "F' ");
MATOUT(CPRIME, 2, 2, 2, 2, "C' ");
MATOUT(E, 2, 2, 2, 2, "E ");
MATOUT(C, 2, 2, 2, 2, "C ");
MATOUT(P, 2, 2, 2, 2, "P ");
}
double DELTA = 0.0;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
}
}
DELTA = sqrt(DELTA / 4.0);
if (IOP != 0) {
printf("\n DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
}
if (DELTA < CRIT) break;
if (ITER < MAXIT) continue;
printf(" NO CONVERGENCE IN SCF");
return;
}
double ENT = EN + ZA * ZB / R;
if (IOP != 0) {
printf("\n\n CALCULATION CONVERGED"
"\n\n ELECTORONIC ENERGY = %20.12lf"
"\n\n TOTAL ENERGY = %20.12f", EN, ENT);
}
if (IOP == 1) {
MATOUT(G, 2, 2, 2, 2, "G ");
MATOUT(F, 2, 2, 2, 2, "F ");
MATOUT(E, 2, 2, 2, 2, "E ");
MATOUT(C, 2, 2, 2, 2, "C ");
MATOUT(P, 2, 2, 2, 2, "P ");
}
MULT(P, S, OLDP, 2, 2);
if (IOP != 0) {
MATOUT(OLDP, 2, 2, 2, 2, "PS ");
}
}
void FORMG(void) {
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
G[i][j] = 0.0;
for (int k = 0; k < 2; ++k) {
for (int l = 0; l < 2; ++l) {
G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
}
}
}
}
}
void DIAG(double F[][2], double C[][2], double E[][2]) {
double theta;
if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
} else {
theta = PI / 4.0;
}
C[0][0] = cos(theta);
C[1][0] = sin(theta);
C[0][1] = sin(theta);
C[1][1] = -cos(theta);
E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
+ F[0][1] * sin(2.0 * theta);
E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
- F[0][1] * sin(2.0 * theta);
E[1][0] = 0.0;
E[0][1] = 0.0;
if (E[1][1] > E[0][0]) return;
double temp = E[1][1];
E[1][1] = E[0][0];
E[0][0] = temp;
temp = C[0][1];
C[0][1] = C[0][0];
C[0][0] = temp;
temp = C[1][1];
C[1][1] = C[1][0];
C[1][0] = temp;
}
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
for (int i = 0; i < M; ++i) {
for (int j = 0; j < M; ++j) {
C[i][j] = 0.0;
for (int k = 0; k < M; ++k) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
}
void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
int IHIGH = 0;
for (;;) {
int LOW = IHIGH + 1;
IHIGH = IHIGH + 5;
IHIGH = (IHIGH < N) ? IHIGH : N;
printf("\n\n\n THE %s ARRAY\n", LABEL);
printf(" ");
for (int i = LOW; i <= IHIGH; ++i) {
printf(" %3d ", i);
}
printf("\n");
for (int i = 1; i <= M; ++i) {
printf("%10d ", i);
for (int j = LOW; j <= IHIGH; ++j) {
printf(" %18.10lf", A[i - 1][j - 1]);
}
printf("\n");
}
if (N <= IHIGH) break;
}
}
最終的な出力は
CALCULATION CONVERGED
ELECTORONIC ENERGY = -4.227529304014
TOTAL ENERGY = -2.860662163500
THE PS ARRAY
1 2
1 1.5296370381 1.1199277981
2 0.6424383868 0.4703629619でした。参照したFortranコードの出力は
CALCULATION CONVERGED
ELECTRONIC ENERGY = -0.422752913203D+01
TOTAL ENERGY = -0.286066199152D+01
THE PS ARRAY
1 2
1 0.1529637121D+01 0.1119927796D+01
2 0.6424383099D+00 0.4703628793D+00でした。
(a)は目的の化合物（図ではIbuprofen)からはじめて分子をばらしていき、全てが既知の入手可能な分子（図では赤で示されている）にまで還元できたら逆合成が完了するというコンセプトを示しています。
(b)は(a)で用いられた既知の反応を示しています。
(c)は(a)の結果得られた反応経路から実際に目的の化合物を合成する過程を示しています。
(d)がこの論文の中心となるア
