1 // SEARCH.CPP
2
3 // Copyright (C) 1998 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 "search.h"
23
24 #include "utility.h"
25 #include "eng1_qm_mopac.h"
26
27 #include "local_i18n.h"
28 #include "notice.h"
29
30 #include <sstream>
31 using namespace std;
32
33 #ifdef WIN32
34 #include <time.h>
35 #endif // WIN32
36
37 /*################################################################################################*/
38
random_search(model * p1,i32s p2,i32s p3,i32s p4,i32s p5,i32s p6)39 random_search::random_search(model * p1, i32s p2, i32s p3, i32s p4, i32s p5, i32s p6)
40 {
41 mdl = p1;
42 molnum = p2;
43 in_crdset = p3;
44 out_crdset = p4;
45
46 cycles = p5;
47 optsteps = p6;
48
49 if (!mdl->IsGroupsClean()) mdl->UpdateGroups(); // intcrd requires this!!!
50 if (!mdl->IsGroupsSorted()) mdl->SortGroups(); // intcrd requires this!!!
51
52 ic = new intcrd((* mdl), molnum, in_crdset);
53
54 eng = mdl->GetCurrentSetup()->GetCurrentEngine();
55 go = NULL;
56
57 counter1 = 0;
58 counter2 = NOT_DEFINED;
59
60 if (!ic->GetVariableCount())
61 {
62 mdl->ErrorMessage(_("ERROR: no rotatable bonds!!!"));
63 counter1 = cycles; // skip the search...
64 }
65
66 CopyCRD(mdl, eng, in_crdset); CopyCRD(eng, mdl, out_crdset);
67 eng->Compute(0); min_energy = eng->energy;
68
69 last_step = NOT_DEFINED;
70 last_E = NOT_DEFINED;
71
72 srand(time(NULL));
73 }
74
~random_search(void)75 random_search::~random_search(void)
76 {
77 if (go != NULL) delete go;
78 delete ic;
79
80 // do not delete eng since it was obtained from setup->GetCurrentEngine()!
81 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
82 }
83
TakeStep(void)84 i32s random_search::TakeStep(void)
85 {
86 last_step = NOT_DEFINED;
87 last_E = NOT_DEFINED;
88
89 if (counter1 < cycles)
90 {
91 if (counter2 == NOT_DEFINED) // start a new cycle...
92 {
93 counter1++; counter2 = 0;
94
95 fGL rotprob = 1.0 / sqrt((fGL) ic->GetVariableCount());
96 for (i32s n1 = 0;n1 < ic->GetVariableCount();n1++)
97 {
98 fGL random;
99
100 random = (fGL) rand() / (fGL) RAND_MAX;
101 if (random > rotprob) continue;
102
103 random = (fGL) rand() / (fGL) RAND_MAX;
104 ic->SetVariable(n1, 2.0 * M_PI * random);
105 }
106
107 ic->UpdateCartesian();
108 mdl->CenterCRDSet(in_crdset, true);
109 CopyCRD(mdl, eng, in_crdset);
110
111 if (go != NULL) delete go;
112 go = new geomopt(eng, 50, 0.005, 10.0);
113 }
114
115 for (i32s n1 = 0;n1 < UPDATEFRQ;n1++) // optimize...
116 {
117 counter2++;
118 go->TakeCGStep(conjugate_gradient::Newton2An);
119 }
120
121 CopyCRD(eng, mdl, in_crdset);
122
123 i32s retval = counter2;
124 if (counter2 >= optsteps)
125 {
126 eng->Compute(0);
127 if (eng->energy < min_energy)
128 {
129 CopyCRD(eng, mdl, out_crdset);
130 min_energy = eng->energy;
131 }
132
133 ostringstream str1;
134 str1 << _("step ") << (counter1 + 1) << "/" << cycles << _(" energy = ") << eng->energy << " kJ/mol" << endl << ends;
135 mdl->PrintToLog(str1.str().c_str());
136
137 last_step = (counter1 + 1);
138 last_E = eng->energy;
139
140 counter2 = NOT_DEFINED;
141 }
142
143 return retval;
144 }
145 else return -1;
146 }
147
148 /*################################################################################################*/
149
systematic_search(model * p1,i32s p2,i32s p3,i32s p4,i32s p5,i32s p6)150 systematic_search::systematic_search(model * p1, i32s p2, i32s p3, i32s p4, i32s p5, i32s p6)
151 {
152 mdl = p1;
153 molnum = p2;
154 in_crdset = p3;
155 out_crdset = p4;
156
157 divisions = p5;
158 optsteps = p6;
159
160 if (!mdl->IsGroupsClean()) mdl->UpdateGroups(); // intcrd requires this!!!
161 if (!mdl->IsGroupsSorted()) mdl->SortGroups(); // intcrd requires this!!!
162
163 ic = new intcrd((* mdl), molnum, in_crdset);
164
165 eng = mdl->GetCurrentSetup()->GetCurrentEngine();
166 go = NULL;
167
168 nvar = ic->GetVariableCount();
169 if (nvar != 0)
170 {
171 counter1 = new i32s[nvar];
172 for (i32s n1 = 0;n1 < nvar;n1++)
173 {
174 counter1[n1] = 0;
175 }
176 }
177 else
178 {
179 mdl->ErrorMessage(_("ERROR: no rotatable bonds!!!"));
180 counter1 = NULL; // skip the search...
181 }
182
183 counter2 = NOT_DEFINED;
184
185 CopyCRD(mdl, eng, in_crdset); CopyCRD(eng, mdl, out_crdset);
186 eng->Compute(0); min_energy = eng->energy;
187 }
188
~systematic_search(void)189 systematic_search::~systematic_search(void)
190 {
191 if (counter1 != NULL) delete[] counter1; // this should never happen...
192
193 if (go != NULL) delete go;
194 delete ic;
195
196 // do not delete eng since it was obtained from setup->GetCurrentEngine()!
197 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
198 }
199
TakeStep(void)200 i32s systematic_search::TakeStep(void)
201 {
202 if (counter1 != NULL)
203 {
204 if (counter2 == NOT_DEFINED) // start a new cycle...
205 {
206 counter2 = 0;
207
208 i32s n1; bool overflow = false;
209 for (n1 = 0;n1 < nvar;n1++)
210 {
211 overflow = false;
212
213 counter1[n1]++;
214 if (counter1[n1] >= divisions)
215 {
216 counter1[n1] = 0;
217 overflow = true;
218 }
219
220 fGL value = (fGL) counter1[n1] / (fGL) divisions;
221 ic->SetVariable(n1, 2.0 * M_PI * value);
222
223 if (!overflow) break;
224 }
225
226 if (overflow && n1 == nvar) // if overflow happened in the last counter, the search is ready!!!
227 {
228 delete[] counter1;
229 counter1 = NULL;
230 }
231
232 ic->UpdateCartesian();
233 mdl->CenterCRDSet(in_crdset, true);
234 CopyCRD(mdl, eng, in_crdset);
235
236 if (go != NULL) delete go;
237 go = new geomopt(eng, 50, 0.005, 10.0);
238 }
239
240 for (i32s n1 = 0;n1 < UPDATEFRQ && counter2 < optsteps;n1++) // optimize...
241 {
242 counter2++;
243 go->TakeCGStep(conjugate_gradient::Newton2An);
244 }
245
246 CopyCRD(eng, mdl, in_crdset);
247
248 i32s retval = counter2;
249 if (counter2 >= optsteps)
250 {
251 eng->Compute(0);
252 if (eng->energy < min_energy)
253 {
254 CopyCRD(eng, mdl, out_crdset);
255 min_energy = eng->energy;
256 }
257
258 if (counter1 != NULL)
259 {
260 stringstream str1;
261 str1 << _("step "); for (int v = 0;v < nvar;v++) { str1 << (char) ('A' + counter1[(nvar - 1) - v]); }
262 str1 << (" energy = ") << eng->energy << " kJ/mol" << endl << ends;
263 mdl->PrintToLog(str1.str().c_str());
264 }
265
266 counter2 = NOT_DEFINED;
267 }
268
269 return retval;
270 }
271 else return -1;
272 }
273
274 /*################################################################################################*/
275
monte_carlo_search(model * p1,i32s p2,i32s p3,i32s p4,i32s p5,i32s p6,i32s p7)276 monte_carlo_search::monte_carlo_search(model * p1, i32s p2, i32s p3, i32s p4, i32s p5, i32s p6, i32s p7)
277 {
278 mdl = p1;
279 molnum = p2;
280 in_crdset = p3;
281 out_crdset = p4;
282
283 n_init_steps = p5;
284 n_simul_steps = p6;
285 optsteps = p7;
286
287 // OPTIMIZE WITH NO NONBONDED TERMS TO GET AN IDEAL SYMMETRIC STARTING GEOMETRY?!?!?
288 // OPTIMIZE WITH NO NONBONDED TERMS TO GET AN IDEAL SYMMETRIC STARTING GEOMETRY?!?!?
289 // OPTIMIZE WITH NO NONBONDED TERMS TO GET AN IDEAL SYMMETRIC STARTING GEOMETRY?!?!?
290
291 if (!mdl->IsGroupsClean()) mdl->UpdateGroups(); // intcrd requires this!!!
292 if (!mdl->IsGroupsSorted()) mdl->SortGroups(); // intcrd requires this!!!
293
294 ic = new intcrd((* mdl), molnum, in_crdset);
295
296 eng = mdl->GetCurrentSetup()->GetCurrentEngine();
297 go = NULL;
298
299 counter1 = -n_init_steps;
300 counter2 = NOT_DEFINED;
301
302 if (!ic->GetVariableCount())
303 {
304 mdl->ErrorMessage(_("ERROR: no rotatable bonds!!!"));
305 counter1 = n_simul_steps; // skip the search...
306 }
307
308 nvar = ic->GetVariableCount();
309 curr_ic1 = new f64[nvar]; curr_ic2 = new f64[nvar];
310 for (i32s n1 = 0;n1 < nvar;n1++)
311 {
312 curr_ic1[n1] = ic->GetVariable(n1);
313 }
314
315 CopyCRD(mdl, eng, in_crdset); CopyCRD(eng, mdl, out_crdset);
316 eng->Compute(0); curr_energy = min_energy = eng->energy;
317
318 srand(time(NULL));
319 }
320
~monte_carlo_search(void)321 monte_carlo_search::~monte_carlo_search(void)
322 {
323 delete[] curr_ic1;
324 delete[] curr_ic2;
325
326 if (go != NULL) delete go;
327 delete ic;
328
329 // do not delete eng since it was obtained from setup->GetCurrentEngine()!
330 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
331 }
332
TakeStep(void)333 i32s monte_carlo_search::TakeStep(void)
334 {
335 if (counter1 < n_simul_steps)
336 {
337 if (counter2 == NOT_DEFINED) // start a new cycle...
338 {
339 counter1++; counter2 = 0;
340
341 const fGL rotprob = 1.0 / sqrt((fGL) nvar);
342 for (i32s n1 = 0;n1 < nvar;n1++)
343 {
344 fGL value = curr_ic1[n1];
345 curr_ic2[n1] = value;
346
347 fGL random;
348
349 random = (fGL) rand() / (fGL) RAND_MAX;
350 if (random > rotprob) continue;
351
352 random = (fGL) rand() / (fGL) RAND_MAX;
353 curr_ic2[n1] = 2.0 * M_PI * random;
354 }
355
356 for (i32s n1 = 0;n1 < nvar;n1++)
357 {
358 ic->SetVariable(n1, curr_ic2[n1]);
359 }
360
361 ic->UpdateCartesian();
362 mdl->CenterCRDSet(in_crdset, true);
363 CopyCRD(mdl, eng, in_crdset);
364
365 if (go != NULL) delete go;
366 go = new geomopt(eng, 50, 0.005, 10.0);
367 }
368
369 for (i32s n1 = 0;n1 < UPDATEFRQ && counter2 < optsteps;n1++) // optimize...
370 {
371 counter2++;
372 go->TakeCGStep(conjugate_gradient::Newton2An);
373 }
374
375 CopyCRD(eng, mdl, in_crdset);
376
377 i32s retval = counter2;
378 if (counter2 >= optsteps)
379 {
380 bool accept = false;
381
382 eng->Compute(0);
383 if (eng->energy < curr_energy) accept = true;
384
385 f64 tmp1 = eng->energy - curr_energy;
386 f64 tmp2 = 1000.0 * tmp1 / (8.31451 * 300.0);
387
388 if (!accept && ((fGL) rand() / (fGL) RAND_MAX < exp(-tmp2))) accept = true;
389 cout << counter1 << " " << eng->energy << " " << curr_energy << _(" TESTVALUE = ") << exp(-tmp2) << endl;
390
391 if (accept)
392 {
393 for (i32s n1 = 0;n1 < nvar;n1++)
394 {
395 curr_ic1[n1] = curr_ic2[n1];
396 }
397
398 curr_energy = eng->energy;
399
400 stringstream str1;
401 str1 << _("step ") << (counter1 + 1) << "/" << n_simul_steps << _(" energy = ") << eng->energy << " kJ/mol" << endl << ends;
402 mdl->PrintToLog(str1.str().c_str());
403
404 // set the weight???
405 // set the weight???
406 // set the weight???
407 }
408 else
409 {
410 counter1--;
411
412 // increase the weighting factor???
413 // increase the weighting factor???
414 // increase the weighting factor???
415 }
416
417 counter2 = NOT_DEFINED;
418
419 if (eng->energy < min_energy)
420 {
421 CopyCRD(eng, mdl, out_crdset);
422 min_energy = eng->energy;
423 }
424 }
425
426 return retval;
427 }
428 else return -1;
429 }
430
431 /*################################################################################################*/
432 /*################################################################################################*/
433 /*################################################################################################*/
434
435 // Szabo A, Ostlund N : "Modern Quantum Chemistry : introduction to advanced electronic
436 // structure theory" McGraw-Hill, 1989
437
transition_state_search(model * p1,f64 p2,f64 p3)438 transition_state_search::transition_state_search(model * p1, f64 p2, f64 p3)
439 {
440 mdl = p1;
441 deltae = p2;
442 fc[0] = fc[1] = p3;
443
444 init_failed = true;
445 target[0] = target[1] = NULL;
446
447 setup * su = mdl->GetCurrentSetup();
448 suQM = dynamic_cast<setup1_qm *>(su);
449
450 engQM = NULL;
451
452 bool bad_setup = false;
453 if (suQM == NULL) bad_setup = true;
454 else
455 {
456 if (suQM->GetCurrEngIndex() < 0) bad_setup = true; // valid values (MOPAC) are currently from 0 to 3!!!
457 if (suQM->GetCurrEngIndex() > 3) bad_setup = true; // valid values (MOPAC) are currently from 0 to 3!!!
458 }
459
460 if (bad_setup)
461 {
462 mdl->ErrorMessage(_("Must use an all-QM/MOPAC setup!\nPlease see the User's Manual."));
463 return;
464 }
465
466 if (mdl->GetAtomCount() % 2 != 0)
467 {
468 mdl->ErrorMessage(_("Atom count must be even!\nPlease see the User's Manual."));
469 return;
470 }
471
472 i32u half = mdl->GetAtomCount() / 2;
473 atom ** atab1 = new ref_atom[half];
474 atom ** atab2 = new ref_atom[half];
475
476 i32u counter1 = 0; i32u counter2 = 0;
477 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
478 {
479 if (counter1 < half)
480 {
481 atab1[counter1] = & (* it1);
482 counter1++;
483 }
484 else
485 {
486 atab2[counter2] = & (* it1);
487 counter2++;
488 }
489 }
490
491 bool identical = true;
492 for (i32u n1 = 0;n1 < half;n1++)
493 {
494 if (atab1[n1]->el.GetAtomicNumber() != atab2[n1]->el.GetAtomicNumber()) identical = false;
495 }
496
497 if (!identical)
498 {
499 delete[] atab1; atab1 = NULL;
500 delete[] atab2; atab2 = NULL;
501
502 mdl->ErrorMessage(_("No proper pair of reactants/products found!\nPlease see the User's Manual."));
503 return;
504 }
505
506 init_failed = false; // ok, it seems that we can go on...
507
508 // copy the coordinates...
509
510 i32s new_csets1 = 2 - mdl->GetCRDSetCount();
511 if (new_csets1 > 0) mdl->PushCRDSets(new_csets1);
512
513 for (i32s n1 = 0;n1 < mdl->GetCRDSetCount();n1++)
514 {
515 mdl->SetCRDSetVisible(n1, true); // debug...
516 }
517
518 for (i32u n1 = 0;n1 < half;n1++)
519 {
520 const fGL * crd = atab2[n1]->GetCRD(0);
521 atab1[n1]->SetCRD(1, crd[0], crd[1], crd[2]);
522 }
523
524 // handle the bonds that are not identical in reactants/products; collect patoms.
525 // TODO : add an option that all atoms are written into patoms??? but keep the current simple system!
526
527 vector<bond *> bvect_r; vector<bond *> bvect_p;
528 for (iter_bl it1 = mdl->GetBondsBegin();it1 != mdl->GetBondsEnd();it1++)
529 {
530 bool bonded_to_r = false;
531 for (i32u n1 = 0;n1 < half;n1++)
532 {
533 if ((* it1).atmr[0] == atab1[n1]) bonded_to_r = true;
534 if ((* it1).atmr[1] == atab1[n1]) bonded_to_r = true;
535 }
536
537 bool bonded_to_p = false;
538 for (i32u n1 = 0;n1 < half;n1++)
539 {
540 if ((* it1).atmr[0] == atab2[n1]) bonded_to_p = true;
541 if ((* it1).atmr[1] == atab2[n1]) bonded_to_p = true;
542 }
543
544 if (bonded_to_r && bonded_to_p) continue; // no such cases should exist, ignored...
545 else
546 {
547 i32s ind1[2] = { NOT_DEFINED, NOT_DEFINED };
548 for (i32u n1 = 0;n1 < half;n1++)
549 {
550 if ((* it1).atmr[0] == atab1[n1]) ind1[0] = n1;
551 if ((* it1).atmr[1] == atab1[n1]) ind1[1] = n1;
552
553 if ((* it1).atmr[0] == atab2[n1]) ind1[0] = n1 + half;
554 if ((* it1).atmr[1] == atab2[n1]) ind1[1] = n1 + half;
555 }
556
557 if (ind1[0] == NOT_DEFINED || ind1[1] == NOT_DEFINED)
558 {
559 assertion_failed(__FILE__, __LINE__, "tss ind1 search failed.");
560 }
561
562 bool has_equivalent = false;
563 for (iter_bl it2 = mdl->GetBondsBegin();it2 != mdl->GetBondsEnd();it2++)
564 {
565 if (it1 == it2) continue;
566
567 i32s ind2[2] = { NOT_DEFINED, NOT_DEFINED };
568 for (i32u n1 = 0;n1 < half;n1++)
569 {
570 if ((* it2).atmr[0] == atab1[n1]) ind2[0] = n1;
571 if ((* it2).atmr[1] == atab1[n1]) ind2[1] = n1;
572
573 if ((* it2).atmr[0] == atab2[n1]) ind2[0] = n1 + half;
574 if ((* it2).atmr[1] == atab2[n1]) ind2[1] = n1 + half;
575 }
576
577 if (ind2[0] == NOT_DEFINED || ind2[1] == NOT_DEFINED)
578 {
579 assertion_failed(__FILE__, __LINE__, "tss ind2 search failed.");
580 }
581
582 if (ind1[0] == ind2[0] + (i32s) half && ind1[1] == ind2[1] + (i32s) half) has_equivalent = true;
583 if (ind1[0] + (i32s) half == ind2[0] && ind1[1] + (i32s) half == ind2[1]) has_equivalent = true;
584 if (ind1[0] == ind2[1] + (i32s) half && ind1[1] == ind2[0] + (i32s) half) has_equivalent = true;
585 if (ind1[0] + (i32s) half == ind2[1] && ind1[1] + (i32s) half == ind2[0]) has_equivalent = true;
586
587 if (has_equivalent && (* it1).bt.GetValue() != (* it2).bt.GetValue()) // for graphics only...
588 {
589 (* it1).bt = bondtype('C');
590 }
591
592 if (has_equivalent) break; // just a speed optimization...
593 }
594
595 if (!has_equivalent)
596 {
597 if (bonded_to_r && !bonded_to_p) bvect_r.push_back(& (* it1));
598 else if (!bonded_to_r && bonded_to_p) bvect_p.push_back(& (* it1));
599 else
600 {
601 assertion_failed(__FILE__, __LINE__, "a fatal error in tss!");
602 }
603 }
604 }
605 }
606
607 for (i32u n1 = 0;n1 < bvect_r.size();n1++)
608 {
609 for (i32u n2 = 0;n2 < 2;n2++)
610 {
611 atom * atmr = bvect_r[n1]->atmr[n2];
612
613 i32s atmi = NOT_DEFINED;
614 for (i32u n3 = 0;n3 < half;n3++)
615 {
616 if (atab1[n3] == atmr) // search reactant side!
617 {
618 atmi = n3;
619 break;
620 }
621 }
622
623 if (atmi == NOT_DEFINED)
624 {
625 assertion_failed(__FILE__, __LINE__, "search r1 failed!");
626 }
627
628 bool unique = true;
629 for (i32u n3 = 0;n3 < patoms.size();n3++)
630 {
631 if (patoms[n3] == (i32u) atmi) unique = false;
632 }
633
634 if (unique) patoms.push_back((i32u) atmi);
635 }
636
637 bond tmpb = bond(bvect_r[n1]->atmr[0], bvect_r[n1]->atmr[1], bondtype());
638 iter_bl itb = find(mdl->GetBondsBegin(), mdl->GetBondsEnd(), tmpb);
639
640 if (itb == mdl->GetBondsEnd())
641 {
642 assertion_failed(__FILE__, __LINE__, "search r2 failed!");
643 }
644
645 rbonds.push_back(& (* itb)); // we are on the reactant side -> the bond is ok!
646 }
647
648 for (i32u n1 = 0;n1 < bvect_p.size();n1++)
649 {
650 i32s stored_atmi[2] = { NOT_DEFINED, NOT_DEFINED };
651
652 for (i32u n2 = 0;n2 < 2;n2++)
653 {
654 atom * atmr = bvect_p[n1]->atmr[n2];
655
656 i32s atmi = NOT_DEFINED;
657 for (i32u n3 = 0;n3 < half;n3++)
658 {
659 if (atab2[n3] == atmr) // search product side!
660 {
661 atmi = n3;
662 break;
663 }
664 }
665
666 stored_atmi[n2] = atmi;
667 if (atmi == NOT_DEFINED)
668 {
669 assertion_failed(__FILE__, __LINE__, "search p1 failed!");
670 }
671
672 bool unique = true;
673 for (i32u n3 = 0;n3 < patoms.size();n3++)
674 {
675 if (patoms[n3] == (i32u) atmi) unique = false;
676 }
677
678 if (unique) patoms.push_back((i32u) atmi);
679 }
680
681 bond tmpb = bond(bvect_p[n1]->atmr[0], bvect_p[n1]->atmr[1], bondtype());
682 iter_bl itb = find(mdl->GetBondsBegin(), mdl->GetBondsEnd(), tmpb);
683
684 if (itb == mdl->GetBondsEnd())
685 {
686 assertion_failed(__FILE__, __LINE__, "search p2 failed!");
687 }
688
689 // at the product side, we need to duplicate the bond in reactant side!!!
690 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
691
692 tmpb = bond(atab1[stored_atmi[0]], atab1[stored_atmi[1]], (* itb).bt);
693 mdl->AddBond(tmpb); pbonds.push_back(& mdl->GetLastBond());
694
695 mdl->RemoveBond(itb); // this is in fact not needed ; would be deleted anyway below!!!
696 }
697
698 // remove the duplicate atoms...
699
700 for (i32u n1 = 0;n1 < half;n1++)
701 {
702 iter_al it1 = mdl->GetAtomsBegin();
703 while (it1 != mdl->GetAtomsEnd())
704 {
705 if (& (* it1) == atab2[n1]) break;
706 else it1++;
707 }
708
709 if (it1 == mdl->GetAtomsEnd())
710 {
711 assertion_failed(__FILE__, __LINE__, "atom search failed!");
712 }
713
714 mdl->RemoveAtom(it1);
715 }
716
717 delete[] atab1; atab1 = NULL;
718 delete[] atab2; atab2 = NULL;
719
720 // ...and then superimpose the initial reactant/product structures.
721
722 superimpose si(mdl, 0, 1);
723 for (i32s n1 = 0;n1 < 100;n1++)
724 {
725 si.TakeCGStep(conjugate_gradient::Newton2An);
726 cout << _("step = ") << n1 << _(" value = ") << si.optval << endl;
727 }
728
729 si.Transform();
730
731 // ok, the geometry setup is ready.
732
733 target[0] = new f64[mdl->GetAtomCount() * 3];
734 target[1] = new f64[mdl->GetAtomCount() * 3];
735
736 SetTarget(0, 1); // use products as target for reactants...
737 SetTarget(1, 0); // use reactants as target for products...
738
739 // if no patoms were found (that is, no change in bonds; for example amine inversion) set the atom that shifts most!
740
741 if (!patoms.size())
742 {
743 i32s index = NOT_DEFINED; fGL maxd = 0.0; i32s counter = 0;
744 for (iter_al it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
745 {
746 const fGL * crd1 = (* it1).GetCRD(0);
747 const fGL * crd2 = (* it1).GetCRD(1);
748
749 fGL dist = 0.0; fGL tmp1;
750 tmp1 = crd2[0] - crd1[0]; dist += tmp1 * tmp1;
751 tmp1 = crd2[1] - crd1[1]; dist += tmp1 * tmp1;
752 tmp1 = crd2[2] - crd1[2]; dist += tmp1 * tmp1;
753 dist = sqrt(dist);
754
755 if (dist > maxd)
756 {
757 maxd = dist;
758 index = counter;
759 }
760
761 counter++;
762 }
763
764 patoms.push_back(index);
765
766 ostringstream txts;
767 txts << _("no patoms found; using ") << index << _(" as a default.") << endl << ends;
768 mdl->PrintToLog(txts.str().c_str());
769
770 cout << txts.str().c_str();
771 }
772
773 // refine the initial structures + other initializations.
774
775 suQM->DiscardCurrentEngine(); // make sure to release mopac_lock!!!
776
777 // 20061025 TH ; the memory management inside libmopac7 seems to be
778 // more or less messed up ; so far we have created and destroyed the
779 // engine object many times here but since versions 1.9x of ghemical
780 // the transition state search has been VERY unstable. valgrind shows
781 // memory mess-ups inside libmopac7 so here we have a design where
782 // the engine object is created and destroyed only once...
783
784 engine * eng = suQM->CreateEngineByIndex(suQM->GetCurrEngIndex());
785 engQM = dynamic_cast<eng1_qm *>(eng); eng = NULL;
786
787 if (engQM == NULL)
788 {
789 assertion_failed(__FILE__, __LINE__, "engQM == NULL");
790 }
791
792 bool problems = false;
793 for (int ii = 0;ii < engQM->GetAtomCount();ii++)
794 {
795 if (engQM->l2g_qm[ii] != ii) problems = true;
796 }
797
798 if (problems)
799 {
800 assertion_failed(__FILE__, __LINE__, "l2g_qm is bad!");
801 }
802
803 for (i32s n1 = 0;n1 < 2;n1++)
804 {
805 CopyCRD(mdl, engQM, n1);
806 engQM->tss_ref_str = target[n1];
807 engQM->tss_force_const = fc[n1];
808
809 geomopt * opt = new geomopt(engQM, 20, 0.0125, 10.0); // optimal settings?!?!?
810
811 engQM->Compute(0);
812 f64 prev = engQM->energy;
813
814 i32s failed_steps = 0;
815 for (i32s n2 = 0;n2 < 100;n2++)
816 {
817 opt->TakeCGStep(conjugate_gradient::Newton2An);
818 cout << n2 << " " << opt->optval << " " << opt->optstp << endl;
819
820 f64 diff = prev - opt->optval;
821 if (diff > 0.001) failed_steps = 0;
822 else failed_steps++;
823
824 prev = opt->optval;
825
826 if (failed_steps >= 20) break;
827 }
828
829 delete opt; opt = NULL;
830
831 // store the calculated energy and the optimized structure...
832
833 engQM->Compute(0);
834 energy1[n1] = engQM->energy;
835 energy2[n1] = engQM->energy - engQM->tss_force_const * engQM->tss_delta_ene;
836 CopyCRD(engQM, mdl, n1);
837 }
838
839 SetTarget(0, 1); // update ; use products as target for reactants...
840 SetTarget(1, 0); // update ; use reactants as target for products...
841
842 t_ene[0] = energy1[0] + deltae;
843 t_ene[1] = energy1[1] + deltae;
844
845 p[0] = -10000.0; // set initial reaction coordinates...
846 p[1] = +10000.0; // ...that are used only to sort the results.
847
848 ready[0] = ready[1] = false;
849 last_de[0] = last_de[1] = 0.0;
850 }
851
~transition_state_search(void)852 transition_state_search::~transition_state_search(void)
853 {
854 if (target[0] != NULL)
855 {
856 delete[] target[0];
857 target[0] = NULL;
858 }
859
860 if (target[1] != NULL)
861 {
862 delete[] target[1];
863 target[1] = NULL;
864 }
865
866 if (engQM != NULL)
867 {
868 #ifdef ENABLE_MOPAC7
869 // 20061025 TH ; even the single deletion of the engine object
870 // seems to cause crashes. :( so try to just leave it there...
871 eng1_qm_mopac * zzz = dynamic_cast<eng1_qm_mopac *>(engQM);
872 if (zzz != NULL) zzz->FixMeLaterTSS(); else delete engQM;
873 #else // ENABLE_MOPAC7
874 delete engQM;
875 #endif // ENABLE_MOPAC7
876
877 engQM = NULL;
878 }
879 }
880
Run(i32s * prev_not_stored)881 void transition_state_search::Run(i32s * prev_not_stored)
882 {
883 if (init_failed)
884 {
885 assertion_failed(__FILE__, __LINE__, "tss init failed!");
886 }
887
888 for (i32s n1 = 0;n1 < 2;n1++)
889 {
890 if (prev_not_stored[n1] == 1) continue;
891
892 // perform constrained optimization...
893
894 CopyCRD(mdl, engQM, n1);
895 engQM->tss_ref_str = target[n1];
896 engQM->tss_force_const = fc[n1];
897
898 geomopt * opt = new geomopt(engQM, 10, 0.0125, 10.0); // optimal settings?!?!?
899
900 i32s n2 = 0;
901 while (true)
902 {
903 opt->TakeCGStep(conjugate_gradient::Newton2An);
904 cout << n2 << " " << opt->optval << " " << opt->optstp << endl;
905
906 if (n2 != 0 && !(n2 % 20))
907 {
908 if (engQM->tss_delta_ene < 1.0e-15) // flipped??? test using zero???
909 {
910 ////////////////////////////////////////////////////////////
911 //cout << "FLIPPED!!! n1= " << n1 << endl; int xx;cin>>xx;
912 ////////////////////////////////////////////////////////////
913 CopyCRD(mdl, engQM, n1);
914 ready[n1] = true;
915 break;
916 }
917
918 const f64 limit = deltae * 0.10;
919 if (opt->optval > (t_ene[n1] - limit) && opt->optval < (t_ene[n1] + limit)) break;
920
921 f64 diff = t_ene[n1] - opt->optval;
922 engQM->tss_force_const += diff / engQM->tss_delta_ene * 0.50; // optimal???
923 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
924 //cout << "n2 = " << n2 << " E = " << opt->optval << " t_ene = " << t_ene[n1] << " fc = " << engQM->tss_force_const << " de = " << engQM->tss_delta_ene << endl; int yy;cin>>yy;
925 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
926 }
927
928 n2++;
929 }
930
931 delete opt; opt = NULL;
932
933 if (n1 == 0) p[n1] += 1.0;
934 if (n1 == 1) p[n1] -= 1.0;
935
936 // copy the optimized structure and update energy...
937
938 engQM->Compute(0);
939 energy1[n1] = engQM->energy;
940 energy2[n1] = engQM->energy - engQM->tss_force_const * engQM->tss_delta_ene;
941 fc[n1] = engQM->tss_force_const; last_de[n1] = engQM->tss_delta_ene;
942 CopyCRD(engQM, mdl, n1);
943 }
944 }
945
UpdateTargets(bool * update)946 void transition_state_search::UpdateTargets(bool * update)
947 {
948 if (init_failed)
949 {
950 assertion_failed(__FILE__, __LINE__, "tss init failed!");
951 }
952
953 if (!update[0] && !update[1]) return;
954
955 bool check_fc[2] = { false, false };
956 for (i32s n1 = 0;n1 < 2;n1++)
957 {
958 if (!update[n1]) continue;
959
960 t_ene[n1] = energy1[n1] + deltae;
961
962 SetTarget(!n1, n1);
963 check_fc[!n1] = true;
964 }
965
966 for (i32s n1 = 0;n1 < 2;n1++)
967 {
968 if (last_de[n1] < 1.0e-15) continue; // flipped??? test using zero???
969 if (!check_fc[n1]) continue;
970
971 ///////////////////////////////////////////////////////
972 //cout << "CHECK FC ; " << n1 << endl; int xx;cin>>xx;
973 ///////////////////////////////////////////////////////
974
975 CopyCRD(mdl, engQM, n1);
976 engQM->tss_ref_str = target[n1];
977 engQM->tss_force_const = fc[n1];
978
979 engQM->Compute(0);
980
981 fc[n1] *= last_de[n1] / engQM->tss_delta_ene;
982 last_de[n1] = engQM->tss_delta_ene;
983 }
984 }
985
SetTarget(i32s index,i32s cset)986 void transition_state_search::SetTarget(i32s index, i32s cset)
987 {
988 if (init_failed)
989 {
990 assertion_failed(__FILE__, __LINE__, "tss init failed!");
991 }
992
993 iter_al it1; i32u counter = 0;
994 for (it1 = mdl->GetAtomsBegin();it1 != mdl->GetAtomsEnd();it1++)
995 {
996 const fGL * crd = (* it1).GetCRD(cset);
997
998 target[index][counter++] = crd[0];
999 target[index][counter++] = crd[1];
1000 target[index][counter++] = crd[2];
1001 }
1002 }
1003
1004 /*################################################################################################*/
1005
1006 // Szabo A, Ostlund N : "Modern Quantum Chemistry : introduction to advanced electronic
1007 // structure theory" McGraw-Hill, 1989
1008
stationary_state_search(engine * p1,i32s p2,f64 p3,f64 p4)1009 stationary_state_search::stationary_state_search(engine * p1, i32s p2, f64 p3, f64 p4) : conjugate_gradient(p2, p3, p4)
1010 {
1011 eng = p1;
1012 SetNGDelta(1.0e-4);
1013
1014 my_d1 = new f64[eng->GetAtomCount() * 3];
1015
1016 for (i32s n1 = 0;n1 < eng->GetAtomCount();n1++)
1017 {
1018 for (i32s n2 = 0;n2 < 3;n2++)
1019 {
1020 AddVar(& eng->crd[n1 * 3 + n2], & my_d1[n1 * 3 + n2]);
1021 }
1022 }
1023 }
1024
~stationary_state_search(void)1025 stationary_state_search::~stationary_state_search(void)
1026 {
1027 delete[] my_d1;
1028 }
1029
GetValue(void)1030 f64 stationary_state_search::GetValue(void)
1031 {
1032 eng->Compute(1); // request energy and gradient!!!
1033
1034 f64 sigma = 0.0;
1035 for (i32s n1 = 0;n1 < eng->GetAtomCount();n1++)
1036 {
1037 for (i32s n2 = 0;n2 < 3;n2++)
1038 {
1039 // the scaling factor 0.01 here is just for conveniency;
1040 // conjgrad/linesearch is easier to handle this way.
1041
1042 f64 tmp1 = eng->d1[n1 * 3 + n2] * 0.01;
1043 sigma += tmp1 * tmp1;
1044 }
1045 }
1046
1047 return sigma;
1048 }
1049
1050 /*################################################################################################*/
1051
1052 // eof
1053