1 // ENG1_QM_MOPAC.CPP
2
3 // Copyright (C) 2001 Tommi Hassinen.
4
5 // This package is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU General Public License as published by
7 // the Free Software Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9
10 // This package is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU General Public License for more details.
14
15 // You should have received a copy of the GNU General Public License
16 // along with this package; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 /*################################################################################################*/
20
21 #include "libghemicalconfig2.h"
22 #include "eng1_qm_mopac.h"
23
24 #include "local_i18n.h"
25
26 #ifdef ENABLE_MOPAC7
27
28 #include <cstring>
29 #include <iomanip>
30 #include <sstream>
31 using namespace std;
32
33 #include <mopac7/libmopac7.h>
34
35 /*################################################################################################*/
36
37 eng1_qm_mopac * eng1_qm_mopac::mopac_lock = NULL;
38
eng1_qm_mopac(setup * p1,i32u p2,i32u mode)39 eng1_qm_mopac::eng1_qm_mopac(setup * p1, i32u p2, i32u mode) : engine(p1, p2), eng1_qm(p1, p2)
40 {
41 if (mopac_lock != NULL)
42 {
43 // the main MOPAC locking mechanism is now in model::CreateDefaultEngine().
44 // if the locking fails in this stage, take it as a serious error and shutdown...
45
46 assertion_failed(__FILE__, __LINE__, "mopac_lock != NULL");
47 }
48 else mopac_lock = this;
49
50 if (GetSetup()->GetModel()->GetConstD_count() > 0)
51 {
52 GetSetup()->GetModel()->WarningMessage(CONSTRAINTS_NOT_SUPPORTED);
53 }
54
55 char fn_FOR005[256] = "FOR005";
56 if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
57
58 cout << _("writing MOPAC-input file ") << fn_FOR005 << endl;
59
60 ofstream file;
61 file.open(fn_FOR005, ios::out);
62
63 file << "XYZ NOLOG ";
64 file << "SCFCRT=0.000001 "; // slightly less strict than default, but seems to stabilize geomopt.
65
66 switch (mode)
67 {
68 case MOPAC_MINDO3:
69 file << "MINDO ";
70 break;
71
72 case MOPAC_AM1:
73 file << "AM1 ";
74 break;
75
76 case MOPAC_PM3:
77 file << "PM3 ";
78 break;
79
80 // if we don't write anything here, the default MNDO method will be used...
81 // if we don't write anything here, the default MNDO method will be used...
82 // if we don't write anything here, the default MNDO method will be used...
83 };
84
85 // one should use UHF for transition state searching???
86 // but it will break the energy level diagrams...
87 // and unfortunately is quite slow...
88
89 //file << "UHF ";
90 //file << "TRIPLET ";
91 //file << "PULAY ";
92
93 file << "CHARGE=" << GetSetup()->GetModel()->GetQMTotalCharge() << " ";
94 file << "GEO-OK MMOK ";
95 file << endl;
96
97 file << "an automatically generated MOPAC input file." << endl << endl;
98
99 // add one dummy atom, to avoid program crash with molecules that have 3 first atoms linearly arranged.
100 // using this dummy atom won't have any other effect, except it cheats a test at GETGEO.F around line 357.
101 // removing this dummy atom and the test at GETGEO.F seems to produce exactly similar results...
102
103 // this problem can be studied with a molecule CH2=C=CH2 that can be drawn in either "right" way (carbons
104 // are the first three atoms) or "wrong" way... produces slightly higher energies and different MOs.
105
106 // HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!
107 // HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!
108 // HOW TO AVOID THIS PROBLEM: draw the molecule in 1-3-2 order instead of 1-2-3 order!!!
109
110 ////////////////////////////////////////////////////////////////////////////////
111 //file << "XX +1.00 1 +1.00 1 +1.00 1" << endl; // this line can be commented out...
112 ////////////////////////////////////////////////////////////////////////////////
113
114 atom ** atmtab = GetSetup()->GetQMAtoms();
115 for (i32s index = 0;index < GetSetup()->GetQMAtomCount();index++)
116 {
117 const fGL * defcrd = atmtab[index]->GetCRD(0);
118 const fGL cf = 10.0; // conversion factor for [nm] -> [�]
119
120 file << atmtab[index]->el.GetSymbol() << "\t";
121 file << setprecision(6) << setw(12) << (defcrd[0] * cf) << " 1 ";
122 file << setprecision(6) << setw(12) << (defcrd[1] * cf) << " 1 ";
123 file << setprecision(6) << setw(12) << (defcrd[2] * cf) << " 1 ";
124 file << endl;
125 }
126
127 file << endl; // add an empty line to mark end of this input???
128 file.close();
129
130 // the MOPAC output is now directed to console; that's because at the beginning of LM7START()
131 // at m7MAIN.c the unit 6 file is not opened. to send output to a logfile, open the unit 6 file.
132
133 lm7start_();
134 lm7_ini_full_xyz();
135 }
136
~eng1_qm_mopac(void)137 eng1_qm_mopac::~eng1_qm_mopac(void)
138 {
139 if (mopac_lock != this) return; // LOCK!!!
140
141 lm7stop_();
142
143 char fn_FOR005[256] = "FOR005";
144 if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
145
146 char fn_SHUTDOWN[256] = "SHUTDOWN";
147 if (getenv("SHUTDOWN") != NULL) strcpy(fn_SHUTDOWN, getenv("SHUTDOWN"));
148
149 cout << _("removing intermediate files... ");
150
151 static ostringstream cs_FOR005;
152 cs_FOR005 << "rm " << fn_FOR005 << ends;
153 system(cs_FOR005.str().c_str());
154
155 static ostringstream cs_SHUTDOWN;
156 cs_SHUTDOWN << "rm " << fn_SHUTDOWN << ends;
157 system(cs_SHUTDOWN.str().c_str());
158
159 cout << "OK!" << endl;
160 mopac_lock = NULL;
161 }
162
FixMeLaterTSS(void)163 void eng1_qm_mopac::FixMeLaterTSS(void)
164 {
165 // this is a temporary fix for some VERY bad problems in TSS...
166 // this is a temporary fix for some VERY bad problems in TSS...
167 // this is a temporary fix for some VERY bad problems in TSS...
168
169 lm7stop_();
170
171 char fn_FOR005[256] = "FOR005";
172 if (getenv("FOR005") != NULL) strcpy(fn_FOR005, getenv("FOR005"));
173 char fn_SHUTDOWN[256] = "SHUTDOWN";
174 if (getenv("SHUTDOWN") != NULL) strcpy(fn_SHUTDOWN, getenv("SHUTDOWN"));
175
176 static ostringstream cs_FOR005;
177 cs_FOR005 << "rm " << fn_FOR005 << ends;
178 system(cs_FOR005.str().c_str());
179 static ostringstream cs_SHUTDOWN;
180 cs_SHUTDOWN << "rm " << fn_SHUTDOWN << ends;
181 system(cs_SHUTDOWN.str().c_str());
182
183 mopac_lock = NULL;
184 }
185
GetElectronCount(void)186 i32s eng1_qm_mopac::GetElectronCount(void)
187 {
188 return lm7_get_electron_count();
189 }
190
GetOrbitalCount(void)191 i32s eng1_qm_mopac::GetOrbitalCount(void)
192 {
193 return lm7_get_orbital_count();
194 }
195
GetOrbitalEnergy(i32s orb_i)196 f64 eng1_qm_mopac::GetOrbitalEnergy(i32s orb_i)
197 {
198 return lm7_get_orbital_energy(orb_i);
199 }
200
Compute(i32u p1,bool)201 void eng1_qm_mopac::Compute(i32u p1, bool)
202 {
203 if (mopac_lock != this) return; // LOCK!!!
204
205 double e; double hf; int i;
206
207 for (i = 0;i < lm7_get_atom_count();i++)
208 {
209 lm7_set_atom_crd(i, & crd[l2g_qm[i] * 3]);
210 }
211
212 if (p1 == 0) // energy was requested...
213 {
214 lm7_call_compfg(& e, & hf, 0);
215 }
216 else if (p1 == 1) // energy and gradient was requested...
217 {
218 lm7_call_compfg(& e, & hf, 1);
219
220 for (i = 0;i < lm7_get_atom_count();i++)
221 {
222 lm7_get_atom_grad(i, & d1[l2g_qm[i] * 3]);
223 }
224 }
225 else // can't calculate higher derivatives just yet...
226 {
227 assertion_failed(__FILE__, __LINE__, "not_implemented");
228 }
229
230 //cout << "energy = " << e << " kJ/mol" << endl;
231 //cout << "heat of formation = " << hf << " kcal/mol" << endl;
232
233 energy = e;
234
235 // the rest is for transition state search only!!!
236 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
237
238 if (tss_ref_str == NULL) return;
239
240 // this is currently implemented for MOPAC only; there is no problems in
241 // making it for MPQC as well, but since semi-empirical methods are quicker
242 // this seems to be enough. ONE CAN USE AB INITIO IN STATIONARY STATE SEARCH!
243
244 tss_delta_ene = 0.0;
245 for (i = 0;i < lm7_get_atom_count();i++)
246 {
247 f64 t1a[3]; f64 t1b = 0.0;
248 for (i32s n5 = 0;n5 < 3;n5++)
249 {
250 f64 t9a = crd[l2g_qm[i] * 3 + n5];
251 f64 t9b = tss_ref_str[i * 3 + n5];
252
253 t1a[n5] = t9a - t9b;
254 t1b += t1a[n5] * t1a[n5];
255 }
256
257 f64 t1c = sqrt(t1b);
258 for (i32s n5 = 0;n5 < 3;n5++)
259 {
260 f64 t9a = t1a[n5] / t1c;
261
262 t1a[n5] = +t9a;
263 }
264
265 // f = a(x)^2
266 // df/dx = 2a(x)
267
268 f64 t2a = tss_force_const * t1b;
269 energy += t2a;
270
271 // df/fa = (x)^2
272
273 tss_delta_ene += t1b;
274
275 if (p1 > 0)
276 {
277 f64 t2b = 2.0 * tss_force_const * t1c;
278
279 for (i32s n5 = 0;n5 < 3;n5++)
280 {
281 f64 t2c = t1a[n5] * t2b;
282
283 // d1[i * 3 + n5] += t2c;
284 d1[l2g_qm[i] * 3 + n5] += t2c;
285 }
286 }
287 }
288 }
289
SetupPlotting(void)290 void eng1_qm_mopac::SetupPlotting(void)
291 {
292 lm7iniplt_();
293 }
294
GetESP(fGL * pp,fGL * dd)295 fGL eng1_qm_mopac::GetESP(fGL * pp, fGL * dd)
296 {
297 if (mopac_lock != this) return 0.0; // LOCK!!!
298
299 const double cf1 = 10.0; // conversion factor for [nm] -> [�]
300 double tmpcrd[3] = { pp[0] * cf1, pp[1] * cf1, pp[2] * cf1 };
301
302 lm7_set_num_potesp(1);
303 lm7_set_crd_potesp(0, tmpcrd);
304
305 getesp_();
306
307 const double cf2 = 2625.5; // conversion factor for [Hartree] -> [kJ/mol]
308 fGL value = lm7_get_val_potesp(0) * cf2;
309
310 if (dd != NULL) // numerical gradient...
311 {
312 fGL old;
313 const fGL delta = 0.0001;
314
315 old = pp[0]; pp[0] += delta;
316 dd[0] = (GetESP(pp, NULL) - value) / delta;
317 pp[0] = old;
318
319 old = pp[1]; pp[1] += delta;
320 dd[1] = (GetESP(pp, NULL) - value) / delta;
321 pp[1] = old;
322
323 old = pp[2]; pp[2] += delta;
324 dd[2] = (GetESP(pp, NULL) - value) / delta;
325 pp[2] = old;
326 }
327
328 return value;
329 }
330
GetElDens(fGL * pp,fGL * dd)331 fGL eng1_qm_mopac::GetElDens(fGL * pp, fGL * dd)
332 {
333 if (mopac_lock != this) return 0.0; // LOCK!!!
334
335 const double cf = 18.897162; // conversion factor for [nm] -> [bohr]
336 double tmpcrd[3] = { pp[0] * cf, pp[1] * cf, pp[2] * cf };
337
338 lm7_set_num_potesp(1);
339 lm7_set_crd_potesp(0, tmpcrd);
340
341 geteldens_();
342
343 fGL value = lm7_get_val_potesp(0);
344
345 if (dd != NULL) // numerical gradient...
346 {
347 fGL old;
348 const fGL delta = 0.0001;
349
350 old = pp[0]; pp[0] += delta;
351 dd[0] = (GetElDens(pp, NULL) - value) / delta;
352 pp[0] = old;
353
354 old = pp[1]; pp[1] += delta;
355 dd[1] = (GetElDens(pp, NULL) - value) / delta;
356 pp[1] = old;
357
358 old = pp[2]; pp[2] += delta;
359 dd[2] = (GetElDens(pp, NULL) - value) / delta;
360 pp[2] = old;
361 }
362
363 return value;
364 }
365
GetOrbital(fGL * pp,fGL * dd)366 fGL eng1_qm_mopac::GetOrbital(fGL * pp, fGL * dd)
367 {
368 if (mopac_lock != this) return 0.0; // LOCK!!!
369
370 const double cf = 18.897162; // conversion factor for [nm] -> [bohr]
371 double tmpcrd[3] = { pp[0] * cf, pp[1] * cf, pp[2] * cf };
372
373 lm7_set_num_potesp(1);
374 lm7_set_crd_potesp(0, tmpcrd);
375
376 lm7_set_plots_orbital_index(GetSetup()->GetModel()->qm_current_orbital + 1);
377
378 getorb_();
379
380 fGL value = lm7_get_val_potesp(0);
381
382 if (dd != NULL) // numerical gradient...
383 {
384 fGL old;
385 const fGL delta = 0.0001;
386
387 old = pp[0]; pp[0] += delta;
388 dd[0] = (GetOrbital(pp, NULL) - value) / delta;
389 pp[0] = old;
390
391 old = pp[1]; pp[1] += delta;
392 dd[1] = (GetOrbital(pp, NULL) - value) / delta;
393 pp[1] = old;
394
395 old = pp[2]; pp[2] += delta;
396 dd[2] = (GetOrbital(pp, NULL) - value) / delta;
397 pp[2] = old;
398 }
399
400 return value;
401 }
402
GetOrbDens(fGL * pp,fGL * dd)403 fGL eng1_qm_mopac::GetOrbDens(fGL * pp, fGL * dd)
404 {
405 if (mopac_lock != this) return 0.0; // LOCK!!!
406
407 // this orbital density plot is closely related to the orbital plot;
408 // we just square the orbital function and multiply it with 2 (assuming 2 electrons per orbital, true for RHF).
409
410 fGL orb = GetOrbital(pp, dd);
411 fGL value = 2.0 * orb * orb;
412
413 if (dd != NULL) // numerical gradient...
414 {
415 fGL old;
416 const fGL delta = 0.0001;
417
418 old = pp[0]; pp[0] += delta;
419 dd[0] = (GetOrbDens(pp, NULL) - value) / delta;
420 pp[0] = old;
421
422 old = pp[1]; pp[1] += delta;
423 dd[1] = (GetOrbDens(pp, NULL) - value) / delta;
424 pp[1] = old;
425
426 old = pp[2]; pp[2] += delta;
427 dd[2] = (GetOrbDens(pp, NULL) - value) / delta;
428 pp[2] = old;
429 }
430
431 return value;
432 }
433
434 /*################################################################################################*/
435
436 #endif // ENABLE_MOPAC7
437
438 // eof
439