Siriius/GFN2-xTB
0
1#include <iostream>
2#include "ulysses/src/GFN.hpp"
3#include "ulysses/src/math/SolverPackage.hpp"
4#include "ulysses/src/Gas.hpp"
5
6int main(int argc, char** argv) {
7
8 //check the files gfn2-xtb_optg.cpp and pm6-corrected.cpp to extend the options here
9
10 //arguments
11 // 0 exe
12 // 1 geometry
13 // 2 charge
14 // 3 Telec = 300
15 // 4 solvation?
16 // 5 solvent name
17 // 6 optimise geometry?
18 // 7 name for new geometry file
19 // 8 thermo?
20 // 9 energy threshold for geometry optimization = 1.0e-6
21 // 10 gradient threshold for geometry optimization = 1.0e-3
22 // 11 calculate density?
23 // 12 name for density file
24 // 13 electronic reactivity indices?
25 // 14 orbital reactivity indices?
26 // 15 Koopman IP?
27 // 16 IP?
28 // 17 EA?
29 // 18 electronativity?
30 // 19 hardness?
31
32 //parameters passed as argument
33 char *p;
34 int charge = strtol(argv[2],&p,10);
35 char *q;
36 double Telec = strtod(argv[3],&q);
37 char *r;
38 int solvation = strtol(argv[4],&r,10);
39 char *v;
40 int optgeom = strtol(argv[6],&v,10);
41 char *s;
42 int thermo = strtol(argv[8],&s,10);
43 char *t;
44 double energy_threshold = strtod(argv[9],&t);
45 char *u;
46 double gradient_threshold = strtod(argv[10],&u);
47 char *w;
48 int calcdensity = strtol(argv[11],&w,10);
49 char *z1;
50 int elecrx = strtol(argv[13],&z1,10);
51 char *z2;
52 int orbrx = strtol(argv[14],&z2,10);
53 char *z3;
54 int koopman = strtol(argv[15],&z3,10);
55 char *z4;
56 int ip = strtol(argv[16],&z4,10);
57 char *z5;
58 int ea = strtol(argv[17],&z5,10);
59 char *z6;
60 int electronegativity = strtol(argv[18],&z6,10);
61 char *z7;
62 int hardness = strtol(argv[19],&z7,10);
63
64 //system declaration
65 std::cout << "running " << argv[1] << "\n";
66 std::cout << "charge = " << charge << std::endl;
67 std::cout << "T electron = " << Telec << std::endl;
68
69 //allocate molecules
70 Molecule Mol1(argv[1],charge,1,"C1");
71
72 //define method and basis set
73 BSet basis(Mol1,"gfn2");
74 GFN2 electron(basis,Mol1);
75 electron.setElectronTemp(Telec);
76
77 //use ALPB solvation?
78 if (solvation > 0) {
79 electron.setSolvent(argv[5]);
80 // "water"
81 // "acetone"
82 // "acetonitrile"
83 // "aniline"
84 // "benzaldehyde"
85 // "benzene"
86 // "dichloromethane"
87 // "chloroform"
88 // "carbon disulfide"
89 // "dioxane"
90 // "dmf"
91 // "dmso"
92 // "ethanol"
93 // "diethyl ether"
94 // "ethyl acetate"
95 // "furane"
96 // "hexadecane"
97 // "hexane"
98 // "methanol"
99 // "nitromethane"
100 // "octanol"
101 // "phenol"
102 // "thf"
103 // "toluene"
104 // "water"
105 // "octanol wet"
106 }
107
108 electron.Calculate(0);
109
110 //optimise geometry?
111 if (optgeom > 0) {
112 BFGSd solve(4,6);
113 SolverOpt(electron,solve,4,0,energy_threshold,gradient_threshold);
114 Molecule Mol2 = electron.Component();
115 Mol2.WriteXYZ(argv[7]);
116 }
117
118 electron.Calculate(1);
119
120 std::cout << std::setprecision(7) << "\n";
121 //perform thermodynamics?
122 if (thermo > 0) {
123 //get vibrations
124 std::vector<double> all_vibrations = electron.CalcVibrFrequencies();
125 int nvibrations = 6;
126 std::vector<double> vibrations(all_vibrations.size() - nvibrations); //CalcVibrFrequencies returns also translation and rotation modes; these must be removed
127 for (size_t idvibr = 0; idvibr < vibrations.size(); ++idvibr) {
128 vibrations[idvibr] = all_vibrations[idvibr + nvibrations];
129 }
130 std::cout << ">all vibrational frequencies" << std::endl;
131 for (size_t idvibr = 0; idvibr < all_vibrations.size(); ++idvibr) {
132 std::cout << all_vibrations[idvibr] << std::endl;
133 }
134 std::cout << "<all vibrational frequencies" << std::endl;
135 //get electronic energies
136 std::vector<double> Eel;
137 Eel.push_back(electron.getEnergy(1)); //the one means that the D3H4X correction is applied to the total energy; use 0 if you want non-corrected energies
138 //get the degeneracy of ground state
139 std::vector<double> gel(1,1.0);
140
141 //get the eigenvalues of inertia matrix
142 std::vector<double> inertia = electron.Component().InertiaEigenvalues();
143
144 double T = 298.15;
145 bool grimmecorrection = true;
146 double numbermolecules = NA; //1 mol
147 double volume = 0.0224;
148 PBlRRlHOE IdealGas(T,argv[1],inertia,vibrations,Eel,gel,charge,1,"C1","0",grimmecorrection,numbermolecules,volume);
149
150 //loop over temperatures and print out
151 double temperature = 100.0; //K
152 std::cout << ">Thermodynamics" << std::endl;
153 for (size_t idx = 0; idx < 2201; ++idx) {
154 IdealGas.changeT(temperature);
155 std::cout << temperature << ";" << IdealGas.S() << ";" << IdealGas.H() << ";" << IdealGas.G() << ";" << IdealGas.U() << ";" << IdealGas.A() << ";" << IdealGas.CP() << ";" << IdealGas.CV() << std::endl;
156 temperature += 0.5;
157 }
158 std::cout << "<Thermodynamics" << std::endl;
159 }
160
161 //get the density?
162 if (calcdensity > 0) {
163 electron.ElectronicDensity(argv[12]);
164 }
165
166 //get charges and polarisabilities
167 std::vector<size_t> atoms = Mol1.Atoms();
168 std::vector<double> AtmCharge = electron.getQAtoms();
169 std::vector<double> polarizabilities;
170 electron.AtomicPolarizabilities(polarizabilities,AtmCharge);
171 std::cout << ">atom;charge;pol\n";
172 for (size_t idx = 0; idx < atoms.size(); ++idx) {
173 std::cout << atoms[idx] << ";";
174 std::cout << AtmCharge[idx] << ";" << polarizabilities[idx] << "\n";
175 }
176 std::cout << "<atom;charge;pol\n";
177 double polbity = 0.0;
178 electron.TotalPolarizability(polbity,AtmCharge);
179 std::cout << " Total Polarizability " << polbity << "\n";
180
181 //additional properties
182 matrixE RxData(1,1);
183 if (elecrx > 0) {
184 electron.ReactivityIndices(RxData,false);
185 std::cout << ">Electronic Reactivity indices" << std::endl;
186 RxData.Print(4);
187 std::cout << "<Electronic Reactivity indices" << std::endl;
188 }
189
190 if (orbrx > 0) {
191 electron.ReactivityIndices(RxData,true);
192 std::cout << ">Orbital Reactivity indices" << std::endl;
193 RxData.Print(4);
194 std::cout << "<Orbital Reactivity indices" << std::endl;
195 }
196
197 if (koopman > 0) {std::cout << "Ionization Potential (Koopman): " << electron.IonizationPotential(true)*au2eV << " eV" << std::endl;}
198 if (ip > 0) {std::cout << "Ionization Potential (Definition): " << electron.IonizationPotential(false)*au2eV << " eV" << std::endl;}
199 if (ea > 0) {std::cout << "Electron Affinity (Definition): " << electron.ElectronAffinity()*au2eV << " eV" << std::endl;}
200
201 if ((electronegativity > 0)||(hardness > 0)) {
202 double chi;
203 double eta;
204 electron.HSABdata(chi,eta);
205 std::cout << "Electronegativity: " << chi*au2eV << " eV" << std::endl;
206 std::cout << "Hardness: " << eta*au2eV << " eV" << std::endl;
207 }
208
209 return 0;
210}
211