-
Notifications
You must be signed in to change notification settings - Fork 1
/
rewrite_S_phi_E.cpp
94 lines (86 loc) · 3.76 KB
/
rewrite_S_phi_E.cpp
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
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cstdlib> // Additions to compile cleanly in 2016...
#include <cstring>
#include "components.h"
#define ZINDO 0
#define HOMO_mon 23 //Counting from zero!
using namespace std;
int main(int argc, char * argv [])
{
if (argc!=4) {
cerr << "***Error*** : Should be called with (in order!):\n"
<< " name of monomer_1 calculations, name of monomer_2 calculations, name of pair calculations\n"
<< " For each monomer and dimer should have a 'log' and a 'pun' file\n";
return (-1);
}
if (ZINDO) { cerr << "Assuming ZINDO calculations\n";}
else { cerr << "Assuming not using ZINDO\n";}
components mon1, mon2, pair;
mon1.readS(argv[1]);
mon1.readMOs(argv[1]);
//mon1.normaliseMOs();
//mon1.checkOrthNorm();
mon2.readS(argv[2]);
mon2.readMOs(argv[2]);
//mon2.normaliseMOs();
pair.readS(argv[3]);
pair.readMOs(argv[3]);
if (mon1.Nbasis!=mon2.Nbasis ) {
cerr << "***WARNING***: Mismatched basis set size for molecula A vs. molecule B. \n";
// return (-2);
}
if (pair.Nbasis!=mon1.Nbasis+mon2.Nbasis) {
cerr << "***WARNING**: The size of the pair basis is not the sum of molecule A + molecule B. This is almost certainly an error. Check your Gaussian jobs! \n";
}
if(!ZINDO) {
//pair.normaliseMOs();
mon1.justPrintS("S_1.txt");
mon1.justPrintMOs("MOs_1.txt");
mon2.justPrintS("S_2.txt");
mon2.justPrintMOs("MOs_2.txt");
pair.justPrintS("S_pair.txt");
pair.justPrintMOs("MOs_pair.txt");
pair.justPrintEvls("Evls_pair.txt");
}
if (ZINDO) {
/********************************************************************************
* If using zindo orbitals are already orthogonal and can simply calculate F_local
********************************************************************************/
double F[pair.Nbasis][pair.Nbasis];
double B[pair.Nbasis][pair.Nbasis]; //B[i][k] = MOs_mons[i][j] * MOs_pair_[j][k]
// = mon1[i][j] * MOs_pair[j][k] for 0 < i,j < N/2
// = mon2[i-N/2][j-N/2] * MOs_pair[j][k] for N/2 < i,j < N
for (int k=0; k<pair.Nbasis; k++) {
for (int i=0; i<mon1.Nbasis; i++) {
B[i][k]=0.0;
for (int j=0; j<mon1.Nbasis; j++) {
B[i][k]+=mon1.MOs[j][i] * pair.MOs[j][k];
}
}
for (int i=mon2.Nbasis; i<pair.Nbasis; i++) {
B[i][k]=0.0;
for (int j=mon2.Nbasis; j<pair.Nbasis; j++) {
B[i][k]+=mon2.MOs[j-mon2.Nbasis][i-mon2.Nbasis] * pair.MOs[j][k];
}
}
}
for (int i=0; i<pair.Nbasis; i++) {
for (int j=0; j<pair.Nbasis; j++) {
F[i][j]=0.0;
for (int k=0; k<pair.Nbasis; k++) {
F[i][j]+=B[i][k] * pair.Evls[k] * B[j][k];
}
}
}
cout << "5-fold HOMO average = "
<< ( F[HOMO_mon][HOMO_mon+mon1.Nbasis] + F[HOMO_mon][HOMO_mon+mon1.Nbasis-1] + F[HOMO_mon][HOMO_mon+mon1.Nbasis-2] + F[HOMO_mon][HOMO_mon+mon1.Nbasis-3] + F[HOMO_mon][HOMO_mon+mon1.Nbasis-4]
+ F[HOMO_mon-1][HOMO_mon+mon1.Nbasis] + F[HOMO_mon-1][HOMO_mon+mon1.Nbasis-1] + F[HOMO_mon-1][HOMO_mon+mon1.Nbasis-2] + F[HOMO_mon-1][HOMO_mon+mon1.Nbasis-3] + F[HOMO_mon-1][HOMO_mon+mon1.Nbasis-4]
+ F[HOMO_mon-2][HOMO_mon+mon1.Nbasis] + F[HOMO_mon-2][HOMO_mon+mon1.Nbasis-1] + F[HOMO_mon-2][HOMO_mon+mon1.Nbasis-2] + F[HOMO_mon-2][HOMO_mon+mon1.Nbasis-3] + F[HOMO_mon-2][HOMO_mon+mon1.Nbasis-4]
+ F[HOMO_mon-3][HOMO_mon+mon1.Nbasis] + F[HOMO_mon-3][HOMO_mon+mon1.Nbasis-1] + F[HOMO_mon-3][HOMO_mon+mon1.Nbasis-2] + F[HOMO_mon-3][HOMO_mon+mon1.Nbasis-3] + F[HOMO_mon-3][HOMO_mon+mon1.Nbasis-4]
+ F[HOMO_mon-4][HOMO_mon+mon1.Nbasis] + F[HOMO_mon-4][HOMO_mon+mon1.Nbasis-1] + F[HOMO_mon-4][HOMO_mon+mon1.Nbasis-2] + F[HOMO_mon-4][HOMO_mon+mon1.Nbasis-3] + F[HOMO_mon-4][HOMO_mon+mon1.Nbasis-4] ) /25.0
<< endl;
}//if(ZINDO)
return 0;
}